1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
2516
2517
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2567
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589
2590
2591
2592
2593
2594
2595
2596
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609
2610
2611
2612
2613
2614
2615
2616
2617
2618
2619
2620
2621
2622
2623
2624
2625
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645
2646
2647
2648
2649
2650
2651
2652
2653
2654
2655
2656
2657
2658
2659
2660
2661
2662
2663
2664
2665
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698
2699
2700
2701
2702
2703
2704
2705
2706
2707
2708
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778
2779
2780
2781
2782
2783
2784
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803
2804
2805
2806
2807
2808
2809
2810
|
/*
* Multiprecision integer arithmetic, implementing mpint.h.
*/
#include <assert.h>
#include <limits.h>
#include <stdio.h>
#include "defs.h"
#include "misc.h"
#include "puttymem.h"
#include "mpint.h"
#include "mpint_i.h"
#define SIZE_T_BITS (CHAR_BIT * sizeof(size_t))
/*
* Inline helpers to take min and max of size_t values, used
* throughout this code.
*/
static inline size_t size_t_min(size_t a, size_t b)
{
return a < b ? a : b;
}
static inline size_t size_t_max(size_t a, size_t b)
{
return a > b ? a : b;
}
/*
* Helper to fetch a word of data from x with array overflow checking.
* If x is too short to have that word, 0 is returned.
*/
static inline BignumInt mp_word(mp_int *x, size_t i)
{
return i < x->nw ? x->w[i] : 0;
}
/*
* Shift an ordinary C integer by BIGNUM_INT_BITS, in a way that
* avoids writing a shift operator whose RHS is greater or equal to
* the size of the type, because that's undefined behaviour in C.
*
* In fact we must avoid even writing it in a definitely-untaken
* branch of an if, because compilers will sometimes warn about
* that. So you can't just write 'shift too big ? 0 : n >> shift',
* because even if 'shift too big' is a constant-expression
* evaluating to false, you can still get complaints about the
* else clause of the ?:.
*
* So we have to re-check _inside_ that clause, so that the shift
* count is reset to something nonsensical but safe in the case
* where the clause wasn't going to be taken anyway.
*/
static uintmax_t shift_right_by_one_word(uintmax_t n)
{
bool shift_too_big = BIGNUM_INT_BYTES >= sizeof(n);
return shift_too_big ? 0 :
n >> (shift_too_big ? 0 : BIGNUM_INT_BITS);
}
static uintmax_t shift_left_by_one_word(uintmax_t n)
{
bool shift_too_big = BIGNUM_INT_BYTES >= sizeof(n);
return shift_too_big ? 0 :
n << (shift_too_big ? 0 : BIGNUM_INT_BITS);
}
mp_int *mp_make_sized(size_t nw)
{
mp_int *x = snew_plus(mp_int, nw * sizeof(BignumInt));
assert(nw); /* we outlaw the zero-word mp_int */
x->nw = nw;
x->w = snew_plus_get_aux(x);
mp_clear(x);
return x;
}
mp_int *mp_new(size_t maxbits)
{
size_t words = (maxbits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
return mp_make_sized(words);
}
mp_int *mp_resize(mp_int *mp, size_t newmaxbits)
{
mp_int *copy = mp_new(newmaxbits);
mp_copy_into(copy, mp);
mp_free(mp);
return copy;
}
mp_int *mp_from_integer(uintmax_t n)
{
mp_int *x = mp_make_sized(
(sizeof(n) + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES);
for (size_t i = 0; i < x->nw; i++)
x->w[i] = n >> (i * BIGNUM_INT_BITS);
return x;
}
size_t mp_max_bytes(mp_int *x)
{
return x->nw * BIGNUM_INT_BYTES;
}
size_t mp_max_bits(mp_int *x)
{
return x->nw * BIGNUM_INT_BITS;
}
void mp_free(mp_int *x)
{
mp_clear(x);
smemclr(x, sizeof(*x));
sfree(x);
}
void mp_dump(FILE *fp, const char *prefix, mp_int *x, const char *suffix)
{
fprintf(fp, "%s0x", prefix);
for (size_t i = mp_max_bytes(x); i-- > 0 ;)
fprintf(fp, "%02X", mp_get_byte(x, i));
fputs(suffix, fp);
}
void mp_copy_into(mp_int *dest, mp_int *src)
{
size_t copy_nw = size_t_min(dest->nw, src->nw);
memmove(dest->w, src->w, copy_nw * sizeof(BignumInt));
smemclr(dest->w + copy_nw, (dest->nw - copy_nw) * sizeof(BignumInt));
}
void mp_copy_integer_into(mp_int *r, uintmax_t n)
{
for (size_t i = 0; i < r->nw; i++) {
r->w[i] = n;
n = shift_right_by_one_word(n);
}
}
/*
* Conditional selection is done by negating 'which', to give a mask
* word which is all 1s if which==1 and all 0s if which==0. Then you
* can select between two inputs a,b without data-dependent control
* flow by XORing them to get their difference; ANDing with the mask
* word to replace that difference with 0 if which==0; and XORing that
* into a, which will either turn it into b or leave it alone.
*
* This trick will be used throughout this code and taken as read the
* rest of the time (or else I'd be here all week typing comments),
* but I felt I ought to explain it in words _once_.
*/
void mp_select_into(mp_int *dest, mp_int *src0, mp_int *src1,
unsigned which)
{
BignumInt mask = -(BignumInt)(1 & which);
for (size_t i = 0; i < dest->nw; i++) {
BignumInt srcword0 = mp_word(src0, i), srcword1 = mp_word(src1, i);
dest->w[i] = srcword0 ^ ((srcword1 ^ srcword0) & mask);
}
}
void mp_cond_swap(mp_int *x0, mp_int *x1, unsigned swap)
{
assert(x0->nw == x1->nw);
volatile BignumInt mask = -(BignumInt)(1 & swap);
for (size_t i = 0; i < x0->nw; i++) {
BignumInt diff = (x0->w[i] ^ x1->w[i]) & mask;
x0->w[i] ^= diff;
x1->w[i] ^= diff;
}
}
void mp_clear(mp_int *x)
{
smemclr(x->w, x->nw * sizeof(BignumInt));
}
void mp_cond_clear(mp_int *x, unsigned clear)
{
BignumInt mask = ~-(BignumInt)(1 & clear);
for (size_t i = 0; i < x->nw; i++)
x->w[i] &= mask;
}
/*
* Common code between mp_from_bytes_{le,be} which reads bytes in an
* arbitrary arithmetic progression.
*/
static mp_int *mp_from_bytes_int(ptrlen bytes, size_t m, size_t c)
{
size_t nw = (bytes.len + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES;
nw = size_t_max(nw, 1);
mp_int *n = mp_make_sized(nw);
for (size_t i = 0; i < bytes.len; i++)
n->w[i / BIGNUM_INT_BYTES] |=
(BignumInt)(((const unsigned char *)bytes.ptr)[m*i+c]) <<
(8 * (i % BIGNUM_INT_BYTES));
return n;
}
mp_int *mp_from_bytes_le(ptrlen bytes)
{
return mp_from_bytes_int(bytes, 1, 0);
}
mp_int *mp_from_bytes_be(ptrlen bytes)
{
return mp_from_bytes_int(bytes, -1, bytes.len - 1);
}
static mp_int *mp_from_words(size_t nw, const BignumInt *w)
{
mp_int *x = mp_make_sized(nw);
memcpy(x->w, w, x->nw * sizeof(BignumInt));
return x;
}
/*
* Decimal-to-binary conversion: just go through the input string
* adding on the decimal value of each digit, and then multiplying the
* number so far by 10.
*/
mp_int *mp_from_decimal_pl(ptrlen decimal)
{
/* 196/59 is an upper bound (and also a continued-fraction
* convergent) for log2(10), so this conservatively estimates the
* number of bits that will be needed to store any number that can
* be written in this many decimal digits. */
assert(decimal.len < (~(size_t)0) / 196);
size_t bits = 196 * decimal.len / 59;
/* Now round that up to words. */
size_t words = bits / BIGNUM_INT_BITS + 1;
mp_int *x = mp_make_sized(words);
for (size_t i = 0; i < decimal.len; i++) {
mp_add_integer_into(x, x, ((const char *)decimal.ptr)[i] - '0');
if (i+1 == decimal.len)
break;
mp_mul_integer_into(x, x, 10);
}
return x;
}
mp_int *mp_from_decimal(const char *decimal)
{
return mp_from_decimal_pl(ptrlen_from_asciz(decimal));
}
/*
* Hex-to-binary conversion: _algorithmically_ simpler than decimal
* (none of those multiplications by 10), but there's some fiddly
* bit-twiddling needed to process each hex digit without diverging
* control flow depending on whether it's a letter or a number.
*/
mp_int *mp_from_hex_pl(ptrlen hex)
{
assert(hex.len <= (~(size_t)0) / 4);
size_t bits = hex.len * 4;
size_t words = (bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
words = size_t_max(words, 1);
mp_int *x = mp_make_sized(words);
for (size_t nibble = 0; nibble < hex.len; nibble++) {
BignumInt digit = ((const char *)hex.ptr)[hex.len-1 - nibble];
BignumInt lmask = ~-((BignumInt)((digit-'a')|('f'-digit))
>> (BIGNUM_INT_BITS-1));
BignumInt umask = ~-((BignumInt)((digit-'A')|('F'-digit))
>> (BIGNUM_INT_BITS-1));
BignumInt digitval = digit - '0';
digitval ^= (digitval ^ (digit - 'a' + 10)) & lmask;
digitval ^= (digitval ^ (digit - 'A' + 10)) & umask;
digitval &= 0xF; /* at least be _slightly_ nice about weird input */
size_t word_idx = nibble / (BIGNUM_INT_BYTES*2);
size_t nibble_within_word = nibble % (BIGNUM_INT_BYTES*2);
x->w[word_idx] |= digitval << (nibble_within_word * 4);
}
return x;
}
mp_int *mp_from_hex(const char *hex)
{
return mp_from_hex_pl(ptrlen_from_asciz(hex));
}
mp_int *mp_copy(mp_int *x)
{
return mp_from_words(x->nw, x->w);
}
uint8_t mp_get_byte(mp_int *x, size_t byte)
{
return 0xFF & (mp_word(x, byte / BIGNUM_INT_BYTES) >>
(8 * (byte % BIGNUM_INT_BYTES)));
}
unsigned mp_get_bit(mp_int *x, size_t bit)
{
return 1 & (mp_word(x, bit / BIGNUM_INT_BITS) >>
(bit % BIGNUM_INT_BITS));
}
uintmax_t mp_get_integer(mp_int *x)
{
uintmax_t toret = 0;
for (size_t i = x->nw; i-- > 0 ;)
toret = shift_left_by_one_word(toret) | x->w[i];
return toret;
}
void mp_set_bit(mp_int *x, size_t bit, unsigned val)
{
size_t word = bit / BIGNUM_INT_BITS;
assert(word < x->nw);
unsigned shift = (bit % BIGNUM_INT_BITS);
x->w[word] &= ~((BignumInt)1 << shift);
x->w[word] |= (BignumInt)(val & 1) << shift;
}
/*
* Helper function used here and there to normalise any nonzero input
* value to 1.
*/
static inline unsigned normalise_to_1(BignumInt n)
{
n = (n >> 1) | (n & 1); /* ensure top bit is clear */
n = (BignumInt)(-n) >> (BIGNUM_INT_BITS - 1); /* normalise to 0 or 1 */
return n;
}
static inline unsigned normalise_to_1_u64(uint64_t n)
{
n = (n >> 1) | (n & 1); /* ensure top bit is clear */
n = (-n) >> 63; /* normalise to 0 or 1 */
return n;
}
/*
* Find the highest nonzero word in a number. Returns the index of the
* word in x->w, and also a pair of output uint64_t in which that word
* appears in the high one shifted left by 'shift_wanted' bits, the
* words immediately below it occupy the space to the right, and the
* words below _that_ fill up the low one.
*
* If there is no nonzero word at all, the passed-by-reference output
* variables retain their original values.
*/
static inline void mp_find_highest_nonzero_word_pair(
mp_int *x, size_t shift_wanted, size_t *index,
uint64_t *hi, uint64_t *lo)
{
uint64_t curr_hi = 0, curr_lo = 0;
for (size_t curr_index = 0; curr_index < x->nw; curr_index++) {
BignumInt curr_word = x->w[curr_index];
unsigned indicator = normalise_to_1(curr_word);
curr_lo = (BIGNUM_INT_BITS < 64 ? (curr_lo >> BIGNUM_INT_BITS) : 0) |
(curr_hi << (64 - BIGNUM_INT_BITS));
curr_hi = (BIGNUM_INT_BITS < 64 ? (curr_hi >> BIGNUM_INT_BITS) : 0) |
((uint64_t)curr_word << shift_wanted);
if (hi) *hi ^= (curr_hi ^ *hi ) & -(uint64_t)indicator;
if (lo) *lo ^= (curr_lo ^ *lo ) & -(uint64_t)indicator;
if (index) *index ^= (curr_index ^ *index) & -(size_t) indicator;
}
}
size_t mp_get_nbits(mp_int *x)
{
/* Sentinel values in case there are no bits set at all: we
* imagine that there's a word at position -1 (i.e. the topmost
* fraction word) which is all 1s, because that way, we handle a
* zero input by considering its highest set bit to be the top one
* of that word, i.e. just below the units digit, i.e. at bit
* index -1, i.e. so we'll return 0 on output. */
size_t hiword_index = -(size_t)1;
uint64_t hiword64 = ~(BignumInt)0;
/*
* Find the highest nonzero word and its index.
*/
mp_find_highest_nonzero_word_pair(x, 0, &hiword_index, &hiword64, NULL);
BignumInt hiword = hiword64; /* in case BignumInt is a narrower type */
/*
* Find the index of the highest set bit within hiword.
*/
BignumInt hibit_index = 0;
for (size_t i = (1 << (BIGNUM_INT_BITS_BITS-1)); i != 0; i >>= 1) {
BignumInt shifted_word = hiword >> i;
BignumInt indicator =
(BignumInt)(-shifted_word) >> (BIGNUM_INT_BITS-1);
hiword ^= (shifted_word ^ hiword ) & -indicator;
hibit_index += i & -(size_t)indicator;
}
/*
* Put together the result.
*/
return (hiword_index << BIGNUM_INT_BITS_BITS) + hibit_index + 1;
}
/*
* Shared code between the hex and decimal output functions to get rid
* of leading zeroes on the output string. The idea is that we wrote
* out a fixed number of digits and a trailing \0 byte into 'buf', and
* now we want to shift it all left so that the first nonzero digit
* moves to buf[0] (or, if there are no nonzero digits at all, we move
* up by 'maxtrim', so that we return 0 as "0" instead of "").
*/
static void trim_leading_zeroes(char *buf, size_t bufsize, size_t maxtrim)
{
size_t trim = maxtrim;
/*
* Look for the first character not equal to '0', to find the
* shift count.
*/
if (trim > 0) {
for (size_t pos = trim; pos-- > 0 ;) {
uint8_t diff = buf[pos] ^ '0';
size_t mask = -((((size_t)diff) - 1) >> (SIZE_T_BITS - 1));
trim ^= (trim ^ pos) & ~mask;
}
}
/*
* Now do the shift, in log n passes each of which does a
* conditional shift by 2^i bytes if bit i is set in the shift
* count.
*/
uint8_t *ubuf = (uint8_t *)buf;
for (size_t logd = 0; bufsize >> logd; logd++) {
uint8_t mask = -(uint8_t)((trim >> logd) & 1);
size_t d = (size_t)1 << logd;
for (size_t i = 0; i+d < bufsize; i++) {
uint8_t diff = mask & (ubuf[i] ^ ubuf[i+d]);
ubuf[i] ^= diff;
ubuf[i+d] ^= diff;
}
}
}
/*
* Binary to decimal conversion. Our strategy here is to extract each
* decimal digit by finding the input number's residue mod 10, then
* subtract that off to give an exact multiple of 10, which then means
* you can safely divide by 10 by means of shifting right one bit and
* then multiplying by the inverse of 5 mod 2^n.
*/
char *mp_get_decimal(mp_int *x_orig)
{
mp_int *x = mp_copy(x_orig), *y = mp_make_sized(x->nw);
/*
* The inverse of 5 mod 2^lots is 0xccccccccccccccccccccd, for an
* appropriate number of 'c's. Manually construct an integer the
* right size.
*/
mp_int *inv5 = mp_make_sized(x->nw);
assert(BIGNUM_INT_BITS % 8 == 0);
for (size_t i = 0; i < inv5->nw; i++)
inv5->w[i] = BIGNUM_INT_MASK / 5 * 4;
inv5->w[0]++;
/*
* 146/485 is an upper bound (and also a continued-fraction
* convergent) of log10(2), so this is a conservative estimate of
* the number of decimal digits needed to store a value that fits
* in this many binary bits.
*/
assert(x->nw < (~(size_t)1) / (146 * BIGNUM_INT_BITS));
size_t bufsize = size_t_max(x->nw * (146 * BIGNUM_INT_BITS) / 485, 1) + 2;
char *outbuf = snewn(bufsize, char);
outbuf[bufsize - 1] = '\0';
/*
* Loop over the number generating digits from the least
* significant upwards, so that we write to outbuf in reverse
* order.
*/
for (size_t pos = bufsize - 1; pos-- > 0 ;) {
/*
* Find the current residue mod 10. We do this by first
* summing the bytes of the number, with all but the lowest
* one multiplied by 6 (because 256^i == 6 mod 10 for all
* i>0). That gives us a single word congruent mod 10 to the
* input number, and then we reduce it further by manual
* multiplication and shifting, just in case the compiler
* target implements the C division operator in a way that has
* input-dependent timing.
*/
uint32_t low_digit = 0, maxval = 0, mult = 1;
for (size_t i = 0; i < x->nw; i++) {
for (unsigned j = 0; j < BIGNUM_INT_BYTES; j++) {
low_digit += mult * (0xFF & (x->w[i] >> (8*j)));
maxval += mult * 0xFF;
mult = 6;
}
/*
* For _really_ big numbers, prevent overflow of t by
* periodically folding the top half of the accumulator
* into the bottom half, using the same rule 'multiply by
* 6 when shifting down by one or more whole bytes'.
*/
if (maxval > UINT32_MAX - (6 * 0xFF * BIGNUM_INT_BYTES)) {
low_digit = (low_digit & 0xFFFF) + 6 * (low_digit >> 16);
maxval = (maxval & 0xFFFF) + 6 * (maxval >> 16);
}
}
/*
* Final reduction of low_digit. We multiply by 2^32 / 10
* (that's the constant 0x19999999) to get a 64-bit value
* whose top 32 bits are the approximate quotient
* low_digit/10; then we subtract off 10 times that; and
* finally we do one last trial subtraction of 10 by adding 6
* (which sets bit 4 if the number was just over 10) and then
* testing bit 4.
*/
low_digit -= 10 * ((0x19999999ULL * low_digit) >> 32);
low_digit -= 10 * ((low_digit + 6) >> 4);
assert(low_digit < 10); /* make sure we did reduce fully */
outbuf[pos] = '0' + low_digit;
/*
* Now subtract off that digit, divide by 2 (using a right
* shift) and by 5 (using the modular inverse), to get the
* next output digit into the units position.
*/
mp_sub_integer_into(x, x, low_digit);
mp_rshift_fixed_into(y, x, 1);
mp_mul_into(x, y, inv5);
}
mp_free(x);
mp_free(y);
mp_free(inv5);
trim_leading_zeroes(outbuf, bufsize, bufsize - 2);
return outbuf;
}
/*
* Binary to hex conversion. Reasonably simple (only a spot of bit
* twiddling to choose whether to output a digit or a letter for each
* nibble).
*/
static char *mp_get_hex_internal(mp_int *x, uint8_t letter_offset)
{
size_t nibbles = x->nw * BIGNUM_INT_BYTES * 2;
size_t bufsize = nibbles + 1;
char *outbuf = snewn(bufsize, char);
outbuf[nibbles] = '\0';
for (size_t nibble = 0; nibble < nibbles; nibble++) {
size_t word_idx = nibble / (BIGNUM_INT_BYTES*2);
size_t nibble_within_word = nibble % (BIGNUM_INT_BYTES*2);
uint8_t digitval = 0xF & (x->w[word_idx] >> (nibble_within_word * 4));
uint8_t mask = -((digitval + 6) >> 4);
char digit = digitval + '0' + (letter_offset & mask);
outbuf[nibbles-1 - nibble] = digit;
}
trim_leading_zeroes(outbuf, bufsize, nibbles - 1);
return outbuf;
}
char *mp_get_hex(mp_int *x)
{
return mp_get_hex_internal(x, 'a' - ('0'+10));
}
char *mp_get_hex_uppercase(mp_int *x)
{
return mp_get_hex_internal(x, 'A' - ('0'+10));
}
/*
* Routines for reading and writing the SSH-1 and SSH-2 wire formats
* for multiprecision integers, declared in marshal.h.
*
* These can't avoid having control flow dependent on the true bit
* size of the number, because the wire format requires the number of
* output bytes to depend on that.
*/
void BinarySink_put_mp_ssh1(BinarySink *bs, mp_int *x)
{
size_t bits = mp_get_nbits(x);
size_t bytes = (bits + 7) / 8;
assert(bits < 0x10000);
put_uint16(bs, bits);
for (size_t i = bytes; i-- > 0 ;)
put_byte(bs, mp_get_byte(x, i));
}
void BinarySink_put_mp_ssh2(BinarySink *bs, mp_int *x)
{
size_t bytes = (mp_get_nbits(x) + 8) / 8;
put_uint32(bs, bytes);
for (size_t i = bytes; i-- > 0 ;)
put_byte(bs, mp_get_byte(x, i));
}
mp_int *BinarySource_get_mp_ssh1(BinarySource *src)
{
unsigned bitc = get_uint16(src);
ptrlen bytes = get_data(src, (bitc + 7) / 8);
if (get_err(src)) {
return mp_from_integer(0);
} else {
mp_int *toret = mp_from_bytes_be(bytes);
/* SSH-1.5 spec says that it's OK for the prefix uint16 to be
* _greater_ than the actual number of bits */
if (mp_get_nbits(toret) > bitc) {
src->err = BSE_INVALID;
mp_free(toret);
toret = mp_from_integer(0);
}
return toret;
}
}
mp_int *BinarySource_get_mp_ssh2(BinarySource *src)
{
ptrlen bytes = get_string(src);
if (get_err(src)) {
return mp_from_integer(0);
} else {
const unsigned char *p = bytes.ptr;
if ((bytes.len > 0 &&
((p[0] & 0x80) ||
(p[0] == 0 && (bytes.len <= 1 || !(p[1] & 0x80)))))) {
src->err = BSE_INVALID;
return mp_from_integer(0);
}
return mp_from_bytes_be(bytes);
}
}
/*
* Make an mp_int structure whose words array aliases a subinterval of
* some other mp_int. This makes it easy to read or write just the low
* or high words of a number, e.g. to add a number starting from a
* high bit position, or to reduce mod 2^{n*BIGNUM_INT_BITS}.
*
* The convention throughout this code is that when we store an mp_int
* directly by value, we always expect it to be an alias of some kind,
* so its words array won't ever need freeing. Whereas an 'mp_int *'
* has an owner, who knows whether it needs freeing or whether it was
* created by address-taking an alias.
*/
static mp_int mp_make_alias(mp_int *in, size_t offset, size_t len)
{
/*
* Bounds-check the offset and length so that we always return
* something valid, even if it's not necessarily the length the
* caller asked for.
*/
if (offset > in->nw)
offset = in->nw;
if (len > in->nw - offset)
len = in->nw - offset;
mp_int toret;
toret.nw = len;
toret.w = in->w + offset;
return toret;
}
/*
* A special case of mp_make_alias: in some cases we preallocate a
* large mp_int to use as scratch space (to avoid pointless
* malloc/free churn in recursive or iterative work).
*
* mp_alloc_from_scratch creates an alias of size 'len' to part of
* 'pool', and adjusts 'pool' itself so that further allocations won't
* overwrite that space.
*
* There's no free function to go with this. Typically you just copy
* the pool mp_int by value, allocate from the copy, and when you're
* done with those allocations, throw the copy away and go back to the
* original value of pool. (A mark/release system.)
*/
static mp_int mp_alloc_from_scratch(mp_int *pool, size_t len)
{
assert(len <= pool->nw);
mp_int toret = mp_make_alias(pool, 0, len);
*pool = mp_make_alias(pool, len, pool->nw);
return toret;
}
/*
* Internal component common to lots of assorted add/subtract code.
* Reads words from a,b; writes into w_out (which might be NULL if the
* output isn't even needed). Takes an input carry flag in 'carry',
* and returns the output carry. Each word read from b is ANDed with
* b_and and then XORed with b_xor.
*
* So you can implement addition by setting b_and to all 1s and b_xor
* to 0; you can subtract by making b_xor all 1s too (effectively
* bit-flipping b) and also passing 1 as the input carry (to turn
* one's complement into two's complement). And you can do conditional
* add/subtract by choosing b_and to be all 1s or all 0s based on a
* condition, because the value of b will be totally ignored if b_and
* == 0.
*/
static BignumCarry mp_add_masked_into(
BignumInt *w_out, size_t rw, mp_int *a, mp_int *b,
BignumInt b_and, BignumInt b_xor, BignumCarry carry)
{
for (size_t i = 0; i < rw; i++) {
BignumInt aword = mp_word(a, i), bword = mp_word(b, i), out;
bword = (bword & b_and) ^ b_xor;
BignumADC(out, carry, aword, bword, carry);
if (w_out)
w_out[i] = out;
}
return carry;
}
/*
* Like the public mp_add_into except that it returns the output carry.
*/
static inline BignumCarry mp_add_into_internal(mp_int *r, mp_int *a, mp_int *b)
{
return mp_add_masked_into(r->w, r->nw, a, b, ~(BignumInt)0, 0, 0);
}
void mp_add_into(mp_int *r, mp_int *a, mp_int *b)
{
mp_add_into_internal(r, a, b);
}
void mp_sub_into(mp_int *r, mp_int *a, mp_int *b)
{
mp_add_masked_into(r->w, r->nw, a, b, ~(BignumInt)0, ~(BignumInt)0, 1);
}
void mp_and_into(mp_int *r, mp_int *a, mp_int *b)
{
for (size_t i = 0; i < r->nw; i++) {
BignumInt aword = mp_word(a, i), bword = mp_word(b, i);
r->w[i] = aword & bword;
}
}
void mp_or_into(mp_int *r, mp_int *a, mp_int *b)
{
for (size_t i = 0; i < r->nw; i++) {
BignumInt aword = mp_word(a, i), bword = mp_word(b, i);
r->w[i] = aword | bword;
}
}
void mp_xor_into(mp_int *r, mp_int *a, mp_int *b)
{
for (size_t i = 0; i < r->nw; i++) {
BignumInt aword = mp_word(a, i), bword = mp_word(b, i);
r->w[i] = aword ^ bword;
}
}
void mp_bic_into(mp_int *r, mp_int *a, mp_int *b)
{
for (size_t i = 0; i < r->nw; i++) {
BignumInt aword = mp_word(a, i), bword = mp_word(b, i);
r->w[i] = aword & ~bword;
}
}
static void mp_cond_negate(mp_int *r, mp_int *x, unsigned yes)
{
BignumCarry carry = yes;
BignumInt flip = -(BignumInt)yes;
for (size_t i = 0; i < r->nw; i++) {
BignumInt xword = mp_word(x, i);
xword ^= flip;
BignumADC(r->w[i], carry, 0, xword, carry);
}
}
/*
* Similar to mp_add_masked_into, but takes a C integer instead of an
* mp_int as the masked operand.
*/
static BignumCarry mp_add_masked_integer_into(
BignumInt *w_out, size_t rw, mp_int *a, uintmax_t b,
BignumInt b_and, BignumInt b_xor, BignumCarry carry)
{
for (size_t i = 0; i < rw; i++) {
BignumInt aword = mp_word(a, i);
BignumInt bword = b;
b = shift_right_by_one_word(b);
BignumInt out;
bword = (bword ^ b_xor) & b_and;
BignumADC(out, carry, aword, bword, carry);
if (w_out)
w_out[i] = out;
}
return carry;
}
void mp_add_integer_into(mp_int *r, mp_int *a, uintmax_t n)
{
mp_add_masked_integer_into(r->w, r->nw, a, n, ~(BignumInt)0, 0, 0);
}
void mp_sub_integer_into(mp_int *r, mp_int *a, uintmax_t n)
{
mp_add_masked_integer_into(r->w, r->nw, a, n,
~(BignumInt)0, ~(BignumInt)0, 1);
}
/*
* Sets r to a + n << (word_index * BIGNUM_INT_BITS), treating
* word_index as secret data.
*/
static void mp_add_integer_into_shifted_by_words(
mp_int *r, mp_int *a, uintmax_t n, size_t word_index)
{
unsigned indicator = 0;
BignumCarry carry = 0;
for (size_t i = 0; i < r->nw; i++) {
/* indicator becomes 1 when we reach the index that the least
* significant bits of n want to be placed at, and it stays 1
* thereafter. */
indicator |= 1 ^ normalise_to_1(i ^ word_index);
/* If indicator is 1, we add the low bits of n into r, and
* shift n down. If it's 0, we add zero bits into r, and
* leave n alone. */
BignumInt bword = n & -(BignumInt)indicator;
uintmax_t new_n = shift_right_by_one_word(n);
n ^= (n ^ new_n) & -(uintmax_t)indicator;
BignumInt aword = mp_word(a, i);
BignumInt out;
BignumADC(out, carry, aword, bword, carry);
r->w[i] = out;
}
}
void mp_mul_integer_into(mp_int *r, mp_int *a, uint16_t n)
{
BignumInt carry = 0, mult = n;
for (size_t i = 0; i < r->nw; i++) {
BignumInt aword = mp_word(a, i);
BignumMULADD(carry, r->w[i], aword, mult, carry);
}
assert(!carry);
}
void mp_cond_add_into(mp_int *r, mp_int *a, mp_int *b, unsigned yes)
{
BignumInt mask = -(BignumInt)(yes & 1);
mp_add_masked_into(r->w, r->nw, a, b, mask, 0, 0);
}
void mp_cond_sub_into(mp_int *r, mp_int *a, mp_int *b, unsigned yes)
{
BignumInt mask = -(BignumInt)(yes & 1);
mp_add_masked_into(r->w, r->nw, a, b, mask, mask, 1 & mask);
}
/*
* Ordered comparison between unsigned numbers is done by subtracting
* one from the other and looking at the output carry.
*/
unsigned mp_cmp_hs(mp_int *a, mp_int *b)
{
size_t rw = size_t_max(a->nw, b->nw);
return mp_add_masked_into(NULL, rw, a, b, ~(BignumInt)0, ~(BignumInt)0, 1);
}
unsigned mp_hs_integer(mp_int *x, uintmax_t n)
{
BignumInt carry = 1;
size_t nwords = sizeof(n)/BIGNUM_INT_BYTES;
for (size_t i = 0, e = size_t_max(x->nw, nwords); i < e; i++) {
BignumInt nword = n;
n = shift_right_by_one_word(n);
BignumInt dummy_out;
BignumADC(dummy_out, carry, mp_word(x, i), ~nword, carry);
(void)dummy_out;
}
return carry;
}
/*
* Equality comparison is done by bitwise XOR of the input numbers,
* ORing together all the output words, and normalising the result
* using our careful normalise_to_1 helper function.
*/
unsigned mp_cmp_eq(mp_int *a, mp_int *b)
{
BignumInt diff = 0;
for (size_t i = 0, limit = size_t_max(a->nw, b->nw); i < limit; i++)
diff |= mp_word(a, i) ^ mp_word(b, i);
return 1 ^ normalise_to_1(diff); /* return 1 if diff _is_ zero */
}
unsigned mp_eq_integer(mp_int *x, uintmax_t n)
{
BignumInt diff = 0;
size_t nwords = sizeof(n)/BIGNUM_INT_BYTES;
for (size_t i = 0, e = size_t_max(x->nw, nwords); i < e; i++) {
BignumInt nword = n;
n = shift_right_by_one_word(n);
diff |= mp_word(x, i) ^ nword;
}
return 1 ^ normalise_to_1(diff); /* return 1 if diff _is_ zero */
}
static void mp_neg_into(mp_int *r, mp_int *a)
{
mp_int zero;
zero.nw = 0;
mp_sub_into(r, &zero, a);
}
mp_int *mp_add(mp_int *x, mp_int *y)
{
mp_int *r = mp_make_sized(size_t_max(x->nw, y->nw) + 1);
mp_add_into(r, x, y);
return r;
}
mp_int *mp_sub(mp_int *x, mp_int *y)
{
mp_int *r = mp_make_sized(size_t_max(x->nw, y->nw));
mp_sub_into(r, x, y);
return r;
}
/*
* Internal routine: multiply and accumulate in the trivial O(N^2)
* way. Sets r <- r + a*b.
*/
static void mp_mul_add_simple(mp_int *r, mp_int *a, mp_int *b)
{
BignumInt *aend = a->w + a->nw, *bend = b->w + b->nw, *rend = r->w + r->nw;
for (BignumInt *ap = a->w, *rp = r->w;
ap < aend && rp < rend; ap++, rp++) {
BignumInt adata = *ap, carry = 0, *rq = rp;
for (BignumInt *bp = b->w; bp < bend && rq < rend; bp++, rq++) {
BignumInt bdata = bp < bend ? *bp : 0;
BignumMULADD2(carry, *rq, adata, bdata, *rq, carry);
}
for (; rq < rend; rq++)
BignumADC(*rq, carry, carry, *rq, 0);
}
}
#ifndef KARATSUBA_THRESHOLD /* allow redefinition via -D for testing */
#define KARATSUBA_THRESHOLD 24
#endif
static inline size_t mp_mul_scratchspace_unary(size_t n)
{
/*
* Simplistic and overcautious bound on the amount of scratch
* space that the recursive multiply function will need.
*
* The rationale is: on the main Karatsuba branch of
* mp_mul_internal, which is the most space-intensive one, we
* allocate space for (a0+a1) and (b0+b1) (each just over half the
* input length n) and their product (the sum of those sizes, i.e.
* just over n itself). Then in order to actually compute the
* product, we do a recursive multiplication of size just over n.
*
* If all those 'just over' weren't there, and everything was
* _exactly_ half the length, you'd get the amount of space for a
* size-n multiply defined by the recurrence M(n) = 2n + M(n/2),
* which is satisfied by M(n) = 4n. But instead it's (2n plus a
* word or two) and M(n/2 plus a word or two). On the assumption
* that there's still some constant k such that M(n) <= kn, this
* gives us kn = 2n + w + k(n/2 + w), where w is a small constant
* (one or two words). That simplifies to kn/2 = 2n + (k+1)w, and
* since we don't even _start_ needing scratch space until n is at
* least 50, we can bound 2n + (k+1)w above by 3n, giving k=6.
*
* So I claim that 6n words of scratch space will suffice, and I
* check that by assertion at every stage of the recursion.
*/
return n * 6;
}
static size_t mp_mul_scratchspace(size_t rw, size_t aw, size_t bw)
{
size_t inlen = size_t_min(rw, size_t_max(aw, bw));
return mp_mul_scratchspace_unary(inlen);
}
static void mp_mul_internal(mp_int *r, mp_int *a, mp_int *b, mp_int scratch)
{
size_t inlen = size_t_min(r->nw, size_t_max(a->nw, b->nw));
assert(scratch.nw >= mp_mul_scratchspace_unary(inlen));
mp_clear(r);
if (inlen < KARATSUBA_THRESHOLD || a->nw == 0 || b->nw == 0) {
/*
* The input numbers are too small to bother optimising. Go
* straight to the simple primitive approach.
*/
mp_mul_add_simple(r, a, b);
return;
}
/*
* Karatsuba divide-and-conquer algorithm. We cut each input in
* half, so that it's expressed as two big 'digits' in a giant
* base D:
*
* a = a_1 D + a_0
* b = b_1 D + b_0
*
* Then the product is of course
*
* ab = a_1 b_1 D^2 + (a_1 b_0 + a_0 b_1) D + a_0 b_0
*
* and we compute the three coefficients by recursively calling
* ourself to do half-length multiplications.
*
* The clever bit that makes this worth doing is that we only need
* _one_ half-length multiplication for the central coefficient
* rather than the two that it obviouly looks like, because we can
* use a single multiplication to compute
*
* (a_1 + a_0) (b_1 + b_0) = a_1 b_1 + a_1 b_0 + a_0 b_1 + a_0 b_0
*
* and then we subtract the other two coefficients (a_1 b_1 and
* a_0 b_0) which we were computing anyway.
*
* Hence we get to multiply two numbers of length N in about three
* times as much work as it takes to multiply numbers of length
* N/2, which is obviously better than the four times as much work
* it would take if we just did a long conventional multiply.
*/
/* Break up the input as botlen + toplen, with botlen >= toplen.
* The 'base' D is equal to 2^{botlen * BIGNUM_INT_BITS}. */
size_t toplen = inlen / 2;
size_t botlen = inlen - toplen;
/* Alias bignums that address the two halves of a,b, and useful
* pieces of r. */
mp_int a0 = mp_make_alias(a, 0, botlen);
mp_int b0 = mp_make_alias(b, 0, botlen);
mp_int a1 = mp_make_alias(a, botlen, toplen);
mp_int b1 = mp_make_alias(b, botlen, toplen);
mp_int r0 = mp_make_alias(r, 0, botlen*2);
mp_int r1 = mp_make_alias(r, botlen, r->nw);
mp_int r2 = mp_make_alias(r, botlen*2, r->nw);
/* Recurse to compute a0*b0 and a1*b1, in their correct positions
* in the output bignum. They can't overlap. */
mp_mul_internal(&r0, &a0, &b0, scratch);
mp_mul_internal(&r2, &a1, &b1, scratch);
if (r->nw < inlen*2) {
/*
* The output buffer isn't large enough to require the whole
* product, so some of a1*b1 won't have been stored. In that
* case we won't try to do the full Karatsuba optimisation;
* we'll just recurse again to compute a0*b1 and a1*b0 - or at
* least as much of them as the output buffer size requires -
* and add each one in.
*/
mp_int s = mp_alloc_from_scratch(
&scratch, size_t_min(botlen+toplen, r1.nw));
mp_mul_internal(&s, &a0, &b1, scratch);
mp_add_into(&r1, &r1, &s);
mp_mul_internal(&s, &a1, &b0, scratch);
mp_add_into(&r1, &r1, &s);
return;
}
/* a0+a1 and b0+b1 */
mp_int asum = mp_alloc_from_scratch(&scratch, botlen+1);
mp_int bsum = mp_alloc_from_scratch(&scratch, botlen+1);
mp_add_into(&asum, &a0, &a1);
mp_add_into(&bsum, &b0, &b1);
/* Their product */
mp_int product = mp_alloc_from_scratch(&scratch, botlen*2+1);
mp_mul_internal(&product, &asum, &bsum, scratch);
/* Subtract off the outer terms we already have */
mp_sub_into(&product, &product, &r0);
mp_sub_into(&product, &product, &r2);
/* And add it in with the right offset. */
mp_add_into(&r1, &r1, &product);
}
void mp_mul_into(mp_int *r, mp_int *a, mp_int *b)
{
mp_int *scratch = mp_make_sized(mp_mul_scratchspace(r->nw, a->nw, b->nw));
mp_mul_internal(r, a, b, *scratch);
mp_free(scratch);
}
mp_int *mp_mul(mp_int *x, mp_int *y)
{
mp_int *r = mp_make_sized(x->nw + y->nw);
mp_mul_into(r, x, y);
return r;
}
void mp_lshift_fixed_into(mp_int *r, mp_int *a, size_t bits)
{
size_t words = bits / BIGNUM_INT_BITS;
size_t bitoff = bits % BIGNUM_INT_BITS;
for (size_t i = r->nw; i-- > 0 ;) {
if (i < words) {
r->w[i] = 0;
} else {
r->w[i] = mp_word(a, i - words);
if (bitoff != 0) {
r->w[i] <<= bitoff;
if (i > words)
r->w[i] |= mp_word(a, i - words - 1) >>
(BIGNUM_INT_BITS - bitoff);
}
}
}
}
void mp_rshift_fixed_into(mp_int *r, mp_int *a, size_t bits)
{
size_t words = bits / BIGNUM_INT_BITS;
size_t bitoff = bits % BIGNUM_INT_BITS;
for (size_t i = 0; i < r->nw; i++) {
r->w[i] = mp_word(a, i + words);
if (bitoff != 0) {
r->w[i] >>= bitoff;
r->w[i] |= mp_word(a, i + words + 1) << (BIGNUM_INT_BITS - bitoff);
}
}
}
mp_int *mp_lshift_fixed(mp_int *x, size_t bits)
{
size_t words = (bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
mp_int *r = mp_make_sized(x->nw + words);
mp_lshift_fixed_into(r, x, bits);
return r;
}
mp_int *mp_rshift_fixed(mp_int *x, size_t bits)
{
size_t words = bits / BIGNUM_INT_BITS;
size_t nw = x->nw - size_t_min(x->nw, words);
mp_int *r = mp_make_sized(size_t_max(nw, 1));
mp_rshift_fixed_into(r, x, bits);
return r;
}
/*
* Safe right shift is done using the same technique as
* trim_leading_zeroes above: you make an n-word left shift by
* composing an appropriate subset of power-of-2-sized shifts, so it
* takes log_2(n) loop iterations each of which does a different shift
* by a power of 2 words, using the usual bit twiddling to make the
* whole shift conditional on the appropriate bit of n.
*/
static void mp_rshift_safe_in_place(mp_int *r, size_t bits)
{
size_t wordshift = bits / BIGNUM_INT_BITS;
size_t bitshift = bits % BIGNUM_INT_BITS;
unsigned clear = (r->nw - wordshift) >> (CHAR_BIT * sizeof(size_t) - 1);
mp_cond_clear(r, clear);
for (unsigned bit = 0; r->nw >> bit; bit++) {
size_t word_offset = (size_t)1 << bit;
BignumInt mask = -(BignumInt)((wordshift >> bit) & 1);
for (size_t i = 0; i < r->nw; i++) {
BignumInt w = mp_word(r, i + word_offset);
r->w[i] ^= (r->w[i] ^ w) & mask;
}
}
/*
* That's done the shifting by words; now we do the shifting by
* bits.
*/
for (unsigned bit = 0; bit < BIGNUM_INT_BITS_BITS; bit++) {
unsigned shift = 1 << bit, upshift = BIGNUM_INT_BITS - shift;
BignumInt mask = -(BignumInt)((bitshift >> bit) & 1);
for (size_t i = 0; i < r->nw; i++) {
BignumInt w = ((r->w[i] >> shift) | (mp_word(r, i+1) << upshift));
r->w[i] ^= (r->w[i] ^ w) & mask;
}
}
}
mp_int *mp_rshift_safe(mp_int *x, size_t bits)
{
mp_int *r = mp_copy(x);
mp_rshift_safe_in_place(r, bits);
return r;
}
void mp_rshift_safe_into(mp_int *r, mp_int *x, size_t bits)
{
mp_copy_into(r, x);
mp_rshift_safe_in_place(r, bits);
}
static void mp_lshift_safe_in_place(mp_int *r, size_t bits)
{
size_t wordshift = bits / BIGNUM_INT_BITS;
size_t bitshift = bits % BIGNUM_INT_BITS;
/*
* Same strategy as mp_rshift_safe_in_place, but of course the
* other way up.
*/
unsigned clear = (r->nw - wordshift) >> (CHAR_BIT * sizeof(size_t) - 1);
mp_cond_clear(r, clear);
for (unsigned bit = 0; r->nw >> bit; bit++) {
size_t word_offset = (size_t)1 << bit;
BignumInt mask = -(BignumInt)((wordshift >> bit) & 1);
for (size_t i = r->nw; i-- > 0 ;) {
BignumInt w = mp_word(r, i - word_offset);
r->w[i] ^= (r->w[i] ^ w) & mask;
}
}
size_t downshift = BIGNUM_INT_BITS - bitshift;
size_t no_shift = (downshift >> BIGNUM_INT_BITS_BITS);
downshift &= ~-(size_t)no_shift;
BignumInt downshifted_mask = ~-(BignumInt)no_shift;
for (size_t i = r->nw; i-- > 0 ;) {
r->w[i] = (r->w[i] << bitshift) |
((mp_word(r, i-1) >> downshift) & downshifted_mask);
}
}
void mp_lshift_safe_into(mp_int *r, mp_int *x, size_t bits)
{
mp_copy_into(r, x);
mp_lshift_safe_in_place(r, bits);
}
void mp_reduce_mod_2to(mp_int *x, size_t p)
{
size_t word = p / BIGNUM_INT_BITS;
size_t mask = ((size_t)1 << (p % BIGNUM_INT_BITS)) - 1;
for (; word < x->nw; word++) {
x->w[word] &= mask;
mask = 0;
}
}
/*
* Inverse mod 2^n is computed by an iterative technique which doubles
* the number of bits at each step.
*/
mp_int *mp_invert_mod_2to(mp_int *x, size_t p)
{
/* Input checks: x must be coprime to the modulus, i.e. odd, and p
* can't be zero */
assert(x->nw > 0);
assert(x->w[0] & 1);
assert(p > 0);
size_t rw = (p + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
rw = size_t_max(rw, 1);
mp_int *r = mp_make_sized(rw);
size_t mul_scratchsize = mp_mul_scratchspace(2*rw, rw, rw);
mp_int *scratch_orig = mp_make_sized(6 * rw + mul_scratchsize);
mp_int scratch_per_iter = *scratch_orig;
mp_int mul_scratch = mp_alloc_from_scratch(
&scratch_per_iter, mul_scratchsize);
r->w[0] = 1;
for (size_t b = 1; b < p; b <<= 1) {
/*
* In each step of this iteration, we have the inverse of x
* mod 2^b, and we want the inverse of x mod 2^{2b}.
*
* Write B = 2^b for convenience, so we want x^{-1} mod B^2.
* Let x = x_0 + B x_1 + k B^2, with 0 <= x_0,x_1 < B.
*
* We want to find r_0 and r_1 such that
* (r_1 B + r_0) (x_1 B + x_0) == 1 (mod B^2)
*
* To begin with, we know r_0 must be the inverse mod B of
* x_0, i.e. of x, i.e. it is the inverse we computed in the
* previous iteration. So now all we need is r_1.
*
* Multiplying out, neglecting multiples of B^2, and writing
* x_0 r_0 = K B + 1, we have
*
* r_1 x_0 B + r_0 x_1 B + K B == 0 (mod B^2)
* => r_1 x_0 B == - r_0 x_1 B - K B (mod B^2)
* => r_1 x_0 == - r_0 x_1 - K (mod B)
* => r_1 == r_0 (- r_0 x_1 - K) (mod B)
*
* (the last step because we multiply through by the inverse
* of x_0, which we already know is r_0).
*/
mp_int scratch_this_iter = scratch_per_iter;
size_t Bw = (b + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
size_t B2w = (2*b + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
/* Start by finding K: multiply x_0 by r_0, and shift down. */
mp_int x0 = mp_alloc_from_scratch(&scratch_this_iter, Bw);
mp_copy_into(&x0, x);
mp_reduce_mod_2to(&x0, b);
mp_int r0 = mp_make_alias(r, 0, Bw);
mp_int Kshift = mp_alloc_from_scratch(&scratch_this_iter, B2w);
mp_mul_internal(&Kshift, &x0, &r0, mul_scratch);
mp_int K = mp_alloc_from_scratch(&scratch_this_iter, Bw);
mp_rshift_fixed_into(&K, &Kshift, b);
/* Now compute the product r_0 x_1, reusing the space of Kshift. */
mp_int x1 = mp_alloc_from_scratch(&scratch_this_iter, Bw);
mp_rshift_fixed_into(&x1, x, b);
mp_reduce_mod_2to(&x1, b);
mp_int r0x1 = mp_make_alias(&Kshift, 0, Bw);
mp_mul_internal(&r0x1, &r0, &x1, mul_scratch);
/* Add K to that. */
mp_add_into(&r0x1, &r0x1, &K);
/* Negate it. */
mp_neg_into(&r0x1, &r0x1);
/* Multiply by r_0. */
mp_int r1 = mp_alloc_from_scratch(&scratch_this_iter, Bw);
mp_mul_internal(&r1, &r0, &r0x1, mul_scratch);
mp_reduce_mod_2to(&r1, b);
/* That's our r_1, so add it on to r_0 to get the full inverse
* output from this iteration. */
mp_lshift_fixed_into(&K, &r1, (b % BIGNUM_INT_BITS));
size_t Bpos = b / BIGNUM_INT_BITS;
mp_int r1_position = mp_make_alias(r, Bpos, B2w-Bpos);
mp_add_into(&r1_position, &r1_position, &K);
}
/* Finally, reduce mod the precise desired number of bits. */
mp_reduce_mod_2to(r, p);
mp_free(scratch_orig);
return r;
}
static size_t monty_scratch_size(MontyContext *mc)
{
return 3*mc->rw + mc->pw + mp_mul_scratchspace(mc->pw, mc->rw, mc->rw);
}
MontyContext *monty_new(mp_int *modulus)
{
MontyContext *mc = snew(MontyContext);
mc->rw = modulus->nw;
mc->rbits = mc->rw * BIGNUM_INT_BITS;
mc->pw = mc->rw * 2 + 1;
mc->m = mp_copy(modulus);
mc->minus_minv_mod_r = mp_invert_mod_2to(mc->m, mc->rbits);
mp_neg_into(mc->minus_minv_mod_r, mc->minus_minv_mod_r);
mp_int *r = mp_make_sized(mc->rw + 1);
r->w[mc->rw] = 1;
mc->powers_of_r_mod_m[0] = mp_mod(r, mc->m);
mp_free(r);
for (size_t j = 1; j < lenof(mc->powers_of_r_mod_m); j++)
mc->powers_of_r_mod_m[j] = mp_modmul(
mc->powers_of_r_mod_m[0], mc->powers_of_r_mod_m[j-1], mc->m);
mc->scratch = mp_make_sized(monty_scratch_size(mc));
return mc;
}
void monty_free(MontyContext *mc)
{
mp_free(mc->m);
for (size_t j = 0; j < 3; j++)
mp_free(mc->powers_of_r_mod_m[j]);
mp_free(mc->minus_minv_mod_r);
mp_free(mc->scratch);
smemclr(mc, sizeof(*mc));
sfree(mc);
}
/*
* The main Montgomery reduction step.
*/
static mp_int monty_reduce_internal(MontyContext *mc, mp_int *x, mp_int scratch)
{
/*
* The trick with Montgomery reduction is that on the one hand we
* want to reduce the size of the input by a factor of about r,
* and on the other hand, the two numbers we just multiplied were
* both stored with an extra factor of r multiplied in. So we
* computed ar*br = ab r^2, but we want to return abr, so we need
* to divide by r - and if we can do that by _actually dividing_
* by r then this also reduces the size of the number.
*
* But we can only do that if the number we're dividing by r is a
* multiple of r. So first we must add an adjustment to it which
* clears its bottom 'rbits' bits. That adjustment must be a
* multiple of m in order to leave the residue mod n unchanged, so
* the question is, what multiple of m can we add to x to make it
* congruent to 0 mod r? And the answer is, x * (-m)^{-1} mod r.
*/
/* x mod r */
mp_int x_lo = mp_make_alias(x, 0, mc->rbits);
/* x * (-m)^{-1}, i.e. the number we want to multiply by m */
mp_int k = mp_alloc_from_scratch(&scratch, mc->rw);
mp_mul_internal(&k, &x_lo, mc->minus_minv_mod_r, scratch);
/* m times that, i.e. the number we want to add to x */
mp_int mk = mp_alloc_from_scratch(&scratch, mc->pw);
mp_mul_internal(&mk, mc->m, &k, scratch);
/* Add it to x */
mp_add_into(&mk, x, &mk);
/* Reduce mod r, by simply making an alias to the upper words of x */
mp_int toret = mp_make_alias(&mk, mc->rw, mk.nw - mc->rw);
/*
* We'll generally be doing this after a multiplication of two
* fully reduced values. So our input could be anything up to m^2,
* and then we added up to rm to it. Hence, the maximum value is
* rm+m^2, and after dividing by r, that becomes r + m(m/r) < 2r.
* So a single trial-subtraction will finish reducing to the
* interval [0,m).
*/
mp_cond_sub_into(&toret, &toret, mc->m, mp_cmp_hs(&toret, mc->m));
return toret;
}
void monty_mul_into(MontyContext *mc, mp_int *r, mp_int *x, mp_int *y)
{
assert(x->nw <= mc->rw);
assert(y->nw <= mc->rw);
mp_int scratch = *mc->scratch;
mp_int tmp = mp_alloc_from_scratch(&scratch, 2*mc->rw);
mp_mul_into(&tmp, x, y);
mp_int reduced = monty_reduce_internal(mc, &tmp, scratch);
mp_copy_into(r, &reduced);
mp_clear(mc->scratch);
}
mp_int *monty_mul(MontyContext *mc, mp_int *x, mp_int *y)
{
mp_int *toret = mp_make_sized(mc->rw);
monty_mul_into(mc, toret, x, y);
return toret;
}
mp_int *monty_modulus(MontyContext *mc)
{
return mc->m;
}
mp_int *monty_identity(MontyContext *mc)
{
return mc->powers_of_r_mod_m[0];
}
mp_int *monty_invert(MontyContext *mc, mp_int *x)
{
/* Given xr, we want to return x^{-1}r = (xr)^{-1} r^2 =
* monty_reduce((xr)^{-1} r^3) */
mp_int *tmp = mp_invert(x, mc->m);
mp_int *toret = monty_mul(mc, tmp, mc->powers_of_r_mod_m[2]);
mp_free(tmp);
return toret;
}
/*
* Importing a number into Montgomery representation involves
* multiplying it by r and reducing mod m. We use the general-purpose
* mp_modmul for this, in case the input number is out of range.
*/
mp_int *monty_import(MontyContext *mc, mp_int *x)
{
return mp_modmul(x, mc->powers_of_r_mod_m[0], mc->m);
}
void monty_import_into(MontyContext *mc, mp_int *r, mp_int *x)
{
mp_int *imported = monty_import(mc, x);
mp_copy_into(r, imported);
mp_free(imported);
}
/*
* Exporting a number means multiplying it by r^{-1}, which is exactly
* what monty_reduce does anyway, so we just do that.
*/
void monty_export_into(MontyContext *mc, mp_int *r, mp_int *x)
{
assert(x->nw <= 2*mc->rw);
mp_int reduced = monty_reduce_internal(mc, x, *mc->scratch);
mp_copy_into(r, &reduced);
mp_clear(mc->scratch);
}
mp_int *monty_export(MontyContext *mc, mp_int *x)
{
mp_int *toret = mp_make_sized(mc->rw);
monty_export_into(mc, toret, x);
return toret;
}
#define MODPOW_LOG2_WINDOW_SIZE 5
#define MODPOW_WINDOW_SIZE (1 << MODPOW_LOG2_WINDOW_SIZE)
mp_int *monty_pow(MontyContext *mc, mp_int *base, mp_int *exponent)
{
/*
* Modular exponentiation is done from the top down, using a
* fixed-window technique.
*
* We have a table storing every power of the base from base^0 up
* to base^{w-1}, where w is a small power of 2, say 2^k. (k is
* defined above as MODPOW_LOG2_WINDOW_SIZE, and w = 2^k is
* defined as MODPOW_WINDOW_SIZE.)
*
* We break the exponent up into k-bit chunks, from the bottom up,
* that is
*
* exponent = c_0 + 2^k c_1 + 2^{2k} c_2 + ... + 2^{nk} c_n
*
* and we compute base^exponent by computing in turn
*
* base^{c_n}
* base^{2^k c_n + c_{n-1}}
* base^{2^{2k} c_n + 2^k c_{n-1} + c_{n-2}}
* ...
*
* where each line is obtained by raising the previous line to the
* power 2^k (i.e. squaring it k times) and then multiplying in
* a value base^{c_i}, which we can look up in our table.
*
* Side-channel considerations: the exponent is secret, so
* actually doing a single table lookup by using a chunk of
* exponent bits as an array index would be an obvious leak of
* secret information into the cache. So instead, in each
* iteration, we read _all_ the table entries, and do a sequence
* of mp_select operations to leave just the one we wanted in the
* variable that will go into the multiplication. In other
* contexts (like software AES) that technique is so prohibitively
* slow that it makes you choose a strategy that doesn't use table
* lookups at all (we do bitslicing in preference); but here, this
* iteration through 2^k table elements is replacing k-1 bignum
* _multiplications_ that you'd have to use instead if you did
* simple square-and-multiply, and that makes it still a win.
*/
/* Table that holds base^0, ..., base^{w-1} */
mp_int *table[MODPOW_WINDOW_SIZE];
table[0] = mp_copy(monty_identity(mc));
for (size_t i = 1; i < MODPOW_WINDOW_SIZE; i++)
table[i] = monty_mul(mc, table[i-1], base);
/* out accumulates the output value */
mp_int *out = mp_make_sized(mc->rw);
mp_copy_into(out, monty_identity(mc));
/* table_entry will hold each value we get out of the table */
mp_int *table_entry = mp_make_sized(mc->rw);
/* Bit index of the chunk of bits we're working on. Start with the
* highest multiple of k strictly less than the size of our
* bignum, i.e. the highest-index chunk of bits that might
* conceivably contain any nonzero bit. */
size_t i = (exponent->nw * BIGNUM_INT_BITS) - 1;
i -= i % MODPOW_LOG2_WINDOW_SIZE;
bool first_iteration = true;
while (true) {
/* Construct the table index */
unsigned table_index = 0;
for (size_t j = 0; j < MODPOW_LOG2_WINDOW_SIZE; j++)
table_index |= mp_get_bit(exponent, i+j) << j;
/* Iterate through the table to do a side-channel-safe lookup,
* ending up with table_entry = table[table_index] */
mp_copy_into(table_entry, table[0]);
for (size_t j = 1; j < MODPOW_WINDOW_SIZE; j++) {
unsigned not_this_one =
((table_index ^ j) + MODPOW_WINDOW_SIZE - 1)
>> MODPOW_LOG2_WINDOW_SIZE;
mp_select_into(table_entry, table[j], table_entry, not_this_one);
}
if (!first_iteration) {
/* Multiply into the output */
monty_mul_into(mc, out, out, table_entry);
} else {
/* On the first iteration, we can save one multiplication
* by just copying */
mp_copy_into(out, table_entry);
first_iteration = false;
}
/* If that was the bottommost chunk of bits, we're done */
if (i == 0)
break;
/* Otherwise, square k times and go round again. */
for (size_t j = 0; j < MODPOW_LOG2_WINDOW_SIZE; j++)
monty_mul_into(mc, out, out, out);
i-= MODPOW_LOG2_WINDOW_SIZE;
}
for (size_t i = 0; i < MODPOW_WINDOW_SIZE; i++)
mp_free(table[i]);
mp_free(table_entry);
mp_clear(mc->scratch);
return out;
}
mp_int *mp_modpow(mp_int *base, mp_int *exponent, mp_int *modulus)
{
assert(modulus->nw > 0);
assert(modulus->w[0] & 1);
MontyContext *mc = monty_new(modulus);
mp_int *m_base = monty_import(mc, base);
mp_int *m_out = monty_pow(mc, m_base, exponent);
mp_int *out = monty_export(mc, m_out);
mp_free(m_base);
mp_free(m_out);
monty_free(mc);
return out;
}
/*
* Given two input integers a,b which are not both even, computes d =
* gcd(a,b) and also two integers A,B such that A*a - B*b = d. A,B
* will be the minimal non-negative pair satisfying that criterion,
* which is equivalent to saying that 0 <= A < b/d and 0 <= B < a/d.
*
* This algorithm is an adapted form of Stein's algorithm, which
* computes gcd(a,b) using only addition and bit shifts (i.e. without
* needing general division), using the following rules:
*
* - if both of a,b are even, divide off a common factor of 2
* - if one of a,b (WLOG a) is even, then gcd(a,b) = gcd(a/2,b), so
* just divide a by 2
* - if both of a,b are odd, then WLOG a>b, and gcd(a,b) =
* gcd(b,(a-b)/2).
*
* Sometimes this function is used for modular inversion, in which
* case we already know we expect the two inputs to be coprime, so to
* save time the 'both even' initial case is assumed not to arise (or
* to have been handled already by the caller). So this function just
* performs a sequence of reductions in the following form:
*
* - if a,b are both odd, sort them so that a > b, and replace a with
* b-a; otherwise sort them so that a is the even one
* - either way, now a is even and b is odd, so divide a by 2.
*
* The big change to Stein's algorithm is that we need the Bezout
* coefficients as output, not just the gcd. So we need to know how to
* generate those in each case, based on the coefficients from the
* reduced pair of numbers:
*
* - If a is even, and u,v are such that u*(a/2) + v*b = d:
* + if u is also even, then this is just (u/2)*a + v*b = d
* + otherwise, (u+b)*(a/2) + (v-a/2)*b is also equal to d, and
* since u and b are both odd, (u+b)/2 is an integer, so we have
* ((u+b)/2)*a + (v-a/2)*b = d.
*
* - If a,b are both odd, and u,v are such that u*b + v*(a-b) = d,
* then v*a + (u-v)*b = d.
*
* In the case where we passed from (a,b) to (b,(a-b)/2), we regard it
* as having first subtracted b from a and then halved a, so both of
* these transformations must be done in sequence.
*
* The code below transforms this from a recursive to an iterative
* algorithm. We first reduce a,b to 0,1, recording at each stage
* whether we did the initial subtraction, and whether we had to swap
* the two values; then we iterate backwards over that record of what
* we did, applying the above rules for building up the Bezout
* coefficients as we go. Of course, all the case analysis is done by
* the usual bit-twiddling conditionalisation to avoid data-dependent
* control flow.
*
* Also, since these mp_ints are generally treated as unsigned, we
* store the coefficients by absolute value, with the semantics that
* they always have opposite sign, and in the unwinding loop we keep a
* bit indicating whether Aa-Bb is currently expected to be +d or -d,
* so that we can do one final conditional adjustment if it's -d.
*
* Once the reduction rules have managed to reduce the input numbers
* to (0,d), then they are stable (the next reduction will always
* divide the even one by 2, which maps 0 to 0). So it doesn't matter
* if we do more steps of the algorithm than necessary; hence, for
* constant time, we just need to find the maximum number we could
* _possibly_ require, and do that many.
*
* If a,b < 2^n, at most 2n iterations are required. Proof: consider
* the quantity Q = log_2(a) + log_2(b). Every step halves one of the
* numbers (and may also reduce one of them further by doing a
* subtraction beforehand, but in the worst case, not by much or not
* at all). So Q reduces by at least 1 per iteration, and it starts
* off with a value at most 2n.
*
* The worst case inputs (I think) are where x=2^{n-1} and y=2^n-1
* (i.e. x is a power of 2 and y is all 1s). In that situation, the
* first n-1 steps repeatedly halve x until it's 1, and then there are
* n further steps each of which subtracts 1 from y and halves it.
*/
static void mp_bezout_into(mp_int *a_coeff_out, mp_int *b_coeff_out,
mp_int *gcd_out, mp_int *a_in, mp_int *b_in)
{
size_t nw = size_t_max(1, size_t_max(a_in->nw, b_in->nw));
/* Make mutable copies of the input numbers */
mp_int *a = mp_make_sized(nw), *b = mp_make_sized(nw);
mp_copy_into(a, a_in);
mp_copy_into(b, b_in);
/* Space to build up the output coefficients, with an extra word
* so that intermediate values can overflow off the top and still
* right-shift back down to the correct value */
mp_int *ac = mp_make_sized(nw + 1), *bc = mp_make_sized(nw + 1);
/* And a general-purpose temp register */
mp_int *tmp = mp_make_sized(nw);
/* Space to record the sequence of reduction steps to unwind. We
* make it a BignumInt for no particular reason except that (a)
* mp_make_sized conveniently zeroes the allocation and mp_free
* wipes it, and (b) this way I can use mp_dump() if I have to
* debug this code. */
size_t steps = 2 * nw * BIGNUM_INT_BITS;
mp_int *record = mp_make_sized(
(steps*2 + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
for (size_t step = 0; step < steps; step++) {
/*
* If a and b are both odd, we want to sort them so that a is
* larger. But if one is even, we want to sort them so that a
* is the even one.
*/
unsigned swap_if_both_odd = mp_cmp_hs(b, a);
unsigned swap_if_one_even = a->w[0] & 1;
unsigned both_odd = a->w[0] & b->w[0] & 1;
unsigned swap = swap_if_one_even ^ (
(swap_if_both_odd ^ swap_if_one_even) & both_odd);
mp_cond_swap(a, b, swap);
/*
* If a,b are both odd, then a is the larger number, so
* subtract the smaller one from it.
*/
mp_cond_sub_into(a, a, b, both_odd);
/*
* Now a is even, so divide it by two.
*/
mp_rshift_fixed_into(a, a, 1);
/*
* Record the two 1-bit values both_odd and swap.
*/
mp_set_bit(record, step*2, both_odd);
mp_set_bit(record, step*2+1, swap);
}
/*
* Now we expect to have reduced the two numbers to 0 and d,
* although we don't know which way round. (But we avoid checking
* this by assertion; sometimes we'll need to do this computation
* without giving away that we already know the inputs were bogus.
* So we'd prefer to just press on and return nonsense.)
*/
if (gcd_out) {
/*
* At this point we can return the actual gcd. Since one of
* a,b is it and the other is zero, the easiest way to get it
* is to add them together.
*/
mp_add_into(gcd_out, a, b);
}
/*
* If the caller _only_ wanted the gcd, and neither Bezout
* coefficient is even required, we can skip the entire unwind
* stage.
*/
if (a_coeff_out || b_coeff_out) {
/*
* The Bezout coefficients of a,b at this point are simply 0
* for whichever of a,b is zero, and 1 for whichever is
* nonzero. The nonzero number equals gcd(a,b), which by
* assumption is odd, so we can do this by just taking the low
* bit of each one.
*/
ac->w[0] = mp_get_bit(a, 0);
bc->w[0] = mp_get_bit(b, 0);
/*
* Overwrite a,b themselves with those same numbers. This has
* the effect of dividing both of them by d, which will
* arrange that during the unwind stage we generate the
* minimal coefficients instead of a larger pair.
*/
mp_copy_into(a, ac);
mp_copy_into(b, bc);
/*
* We'll maintain the invariant as we unwind that ac * a - bc
* * b is either +d or -d (or rather, +1/-1 after scaling by
* d), and we'll remember which. (We _could_ keep it at +d the
* whole time, but it would cost more work every time round
* the loop, so it's cheaper to fix that up once at the end.)
*
* Initially, the result is +d if a was the nonzero value after
* reduction, and -d if b was.
*/
unsigned minus_d = b->w[0];
for (size_t step = steps; step-- > 0 ;) {
/*
* Recover the data from the step we're unwinding.
*/
unsigned both_odd = mp_get_bit(record, step*2);
unsigned swap = mp_get_bit(record, step*2+1);
/*
* Unwind the division: if our coefficient of a is odd, we
* adjust the coefficients by +b and +a respectively.
*/
unsigned adjust = ac->w[0] & 1;
mp_cond_add_into(ac, ac, b, adjust);
mp_cond_add_into(bc, bc, a, adjust);
/*
* Now ac is definitely even, so we divide it by two.
*/
mp_rshift_fixed_into(ac, ac, 1);
/*
* Now unwind the subtraction, if there was one, by adding
* ac to bc.
*/
mp_cond_add_into(bc, bc, ac, both_odd);
/*
* Undo the transformation of the input numbers, by
* multiplying a by 2 and then adding b to a (the latter
* only if both_odd).
*/
mp_lshift_fixed_into(a, a, 1);
mp_cond_add_into(a, a, b, both_odd);
/*
* Finally, undo the swap. If we do swap, this also
* reverses the sign of the current result ac*a+bc*b.
*/
mp_cond_swap(a, b, swap);
mp_cond_swap(ac, bc, swap);
minus_d ^= swap;
}
/*
* Now we expect to have recovered the input a,b (or rather,
* the versions of them divided by d). But we might find that
* our current result is -d instead of +d, that is, we have
* A',B' such that A'a - B'b = -d.
*
* In that situation, we set A = b-A' and B = a-B', giving us
* Aa-Bb = ab - A'a - ab + B'b = +1.
*/
mp_sub_into(tmp, b, ac);
mp_select_into(ac, ac, tmp, minus_d);
mp_sub_into(tmp, a, bc);
mp_select_into(bc, bc, tmp, minus_d);
/*
* Now we really are done. Return the outputs.
*/
if (a_coeff_out)
mp_copy_into(a_coeff_out, ac);
if (b_coeff_out)
mp_copy_into(b_coeff_out, bc);
}
mp_free(a);
mp_free(b);
mp_free(ac);
mp_free(bc);
mp_free(tmp);
mp_free(record);
}
mp_int *mp_invert(mp_int *x, mp_int *m)
{
mp_int *result = mp_make_sized(m->nw);
mp_bezout_into(result, NULL, NULL, x, m);
return result;
}
void mp_gcd_into(mp_int *a, mp_int *b, mp_int *gcd, mp_int *A, mp_int *B)
{
/*
* Identify shared factors of 2. To do this we OR the two numbers
* to get something whose lowest set bit is in the right place,
* remove all higher bits by ANDing it with its own negation, and
* use mp_get_nbits to find the location of the single remaining
* set bit.
*/
mp_int *tmp = mp_make_sized(size_t_max(a->nw, b->nw));
for (size_t i = 0; i < tmp->nw; i++)
tmp->w[i] = mp_word(a, i) | mp_word(b, i);
BignumCarry carry = 1;
for (size_t i = 0; i < tmp->nw; i++) {
BignumInt negw;
BignumADC(negw, carry, 0, ~tmp->w[i], carry);
tmp->w[i] &= negw;
}
size_t shift = mp_get_nbits(tmp) - 1;
mp_free(tmp);
/*
* Make copies of a,b with those shared factors of 2 divided off,
* so that at least one is odd (which is the precondition for
* mp_bezout_into). Compute the gcd of those.
*/
mp_int *as = mp_rshift_safe(a, shift);
mp_int *bs = mp_rshift_safe(b, shift);
mp_bezout_into(A, B, gcd, as, bs);
mp_free(as);
mp_free(bs);
/*
* And finally shift the gcd back up (unless the caller didn't
* even ask for it), to put the shared factors of 2 back in.
*/
if (gcd)
mp_lshift_safe_in_place(gcd, shift);
}
mp_int *mp_gcd(mp_int *a, mp_int *b)
{
mp_int *gcd = mp_make_sized(size_t_min(a->nw, b->nw));
mp_gcd_into(a, b, gcd, NULL, NULL);
return gcd;
}
unsigned mp_coprime(mp_int *a, mp_int *b)
{
mp_int *gcd = mp_gcd(a, b);
unsigned toret = mp_eq_integer(gcd, 1);
mp_free(gcd);
return toret;
}
static uint32_t recip_approx_32(uint32_t x)
{
/*
* Given an input x in [2^31,2^32), i.e. a uint32_t with its high
* bit set, this function returns an approximation to 2^63/x,
* computed using only multiplications and bit shifts just in case
* the C divide operator has non-constant time (either because the
* underlying machine instruction does, or because the operator
* expands to a library function on a CPU without hardware
* division).
*
* The coefficients are derived from those of the degree-9
* polynomial which is the minimax-optimal approximation to that
* function on the given interval (generated using the Remez
* algorithm), converted into integer arithmetic with shifts used
* to maximise the number of significant bits at every state. (A
* sort of 'static floating point' - the exponent is statically
* known at every point in the code, so it never needs to be
* stored at run time or to influence runtime decisions.)
*
* Exhaustive iteration over the whole input space shows the
* largest possible error to be 1686.54. (The input value
* attaining that bound is 4226800006 == 0xfbefd986, whose true
* reciprocal is 2182116973.540... == 0x8210766d.8a6..., whereas
* this function returns 2182115287 == 0x82106fd7.)
*/
uint64_t r = 0x92db03d6ULL;
r = 0xf63e71eaULL - ((r*x) >> 34);
r = 0xb63721e8ULL - ((r*x) >> 34);
r = 0x9c2da00eULL - ((r*x) >> 33);
r = 0xaada0bb8ULL - ((r*x) >> 32);
r = 0xf75cd403ULL - ((r*x) >> 31);
r = 0xecf97a41ULL - ((r*x) >> 31);
r = 0x90d876cdULL - ((r*x) >> 31);
r = 0x6682799a0ULL - ((r*x) >> 26);
return r;
}
void mp_divmod_into(mp_int *n, mp_int *d, mp_int *q_out, mp_int *r_out)
{
assert(!mp_eq_integer(d, 0));
/*
* We do division by using Newton-Raphson iteration to converge to
* the reciprocal of d (or rather, R/d for R a sufficiently large
* power of 2); then we multiply that reciprocal by n; and we
* finish up with conditional subtraction.
*
* But we have to do it in a fixed number of N-R iterations, so we
* need some error analysis to know how many we might need.
*
* The iteration is derived by defining f(r) = d - R/r.
* Differentiating gives f'(r) = R/r^2, and the Newton-Raphson
* formula applied to those functions gives
*
* r_{i+1} = r_i - f(r_i) / f'(r_i)
* = r_i - (d - R/r_i) r_i^2 / R
* = r_i (2 R - d r_i) / R
*
* Now let e_i be the error in a given iteration, in the sense
* that
*
* d r_i = R + e_i
* i.e. e_i/R = (r_i - r_true) / r_true
*
* so e_i is the _relative_ error in r_i.
*
* We must also introduce a rounding-error term, because the
* division by R always gives an integer. This might make the
* output off by up to 1 (in the negative direction, because
* right-shifting gives floor of the true quotient). So when we
* divide by R, we must imagine adding some f in [0,1). Then we
* have
*
* d r_{i+1} = d r_i (2 R - d r_i) / R - d f
* = (R + e_i) (R - e_i) / R - d f
* = (R^2 - e_i^2) / R - d f
* = R - (e_i^2 / R + d f)
* => e_{i+1} = - (e_i^2 / R + d f)
*
* The sum of two positive quantities is bounded above by twice
* their max, and max |f| = 1, so we can bound this as follows:
*
* |e_{i+1}| <= 2 max (e_i^2/R, d)
* |e_{i+1}/R| <= 2 max ((e_i/R)^2, d/R)
* log2 |R/e_{i+1}| <= min (2 log2 |R/e_i|, log2 |R/d|) - 1
*
* which tells us that the number of 'good' bits - i.e.
* log2(R/e_i) - very nearly doubles at every iteration (apart
* from that subtraction of 1), until it gets to the same size as
* log2(R/d). In other words, the size of R in bits has to be the
* size of denominator we're putting in, _plus_ the amount of
* precision we want to get back out.
*
* So when we multiply n (the input numerator) by our final
* reciprocal approximation r, but actually r differs from R/d by
* up to 2, then it follows that
*
* n/d - nr/R = n/d - [ n (R/d + e) ] / R
* = n/d - [ (n/d) R + n e ] / R
* = -ne/R
* => 0 <= n/d - nr/R < 2n/R
*
* so our computed quotient can differ from the true n/d by up to
* 2n/R. Hence, as long as we also choose R large enough that 2n/R
* is bounded above by a constant, we can guarantee a bounded
* number of final conditional-subtraction steps.
*/
/*
* Get at least 32 of the most significant bits of the input
* number.
*/
size_t hiword_index = 0;
uint64_t hibits = 0, lobits = 0;
mp_find_highest_nonzero_word_pair(d, 64 - BIGNUM_INT_BITS,
&hiword_index, &hibits, &lobits);
/*
* Make a shifted combination of those two words which puts the
* topmost bit of the number at bit 63.
*/
size_t shift_up = 0;
for (size_t i = BIGNUM_INT_BITS_BITS; i-- > 0;) {
size_t sl = (size_t)1 << i; /* left shift count */
size_t sr = 64 - sl; /* complementary right-shift count */
/* Should we shift up? */
unsigned indicator = 1 ^ normalise_to_1_u64(hibits >> sr);
/* If we do, what will we get? */
uint64_t new_hibits = (hibits << sl) | (lobits >> sr);
uint64_t new_lobits = lobits << sl;
size_t new_shift_up = shift_up + sl;
/* Conditionally swap those values in. */
hibits ^= (hibits ^ new_hibits ) & -(uint64_t)indicator;
lobits ^= (lobits ^ new_lobits ) & -(uint64_t)indicator;
shift_up ^= (shift_up ^ new_shift_up ) & -(size_t) indicator;
}
/*
* So now we know the most significant 32 bits of d are at the top
* of hibits. Approximate the reciprocal of those bits.
*/
lobits = (uint64_t)recip_approx_32(hibits >> 32) << 32;
hibits = 0;
/*
* And shift that up by as many bits as the input was shifted up
* just now, so that the product of this approximation and the
* actual input will be close to a fixed power of two regardless
* of where the MSB was.
*
* I do this in another log n individual passes, partly in case
* the CPU's register-controlled shift operation isn't
* time-constant, and also in case the compiler code-generates
* uint64_t shifts out of a variable number of smaller-word shift
* instructions, e.g. by splitting up into cases.
*/
for (size_t i = BIGNUM_INT_BITS_BITS; i-- > 0;) {
size_t sl = (size_t)1 << i; /* left shift count */
size_t sr = 64 - sl; /* complementary right-shift count */
/* Should we shift up? */
unsigned indicator = 1 & (shift_up >> i);
/* If we do, what will we get? */
uint64_t new_hibits = (hibits << sl) | (lobits >> sr);
uint64_t new_lobits = lobits << sl;
/* Conditionally swap those values in. */
hibits ^= (hibits ^ new_hibits ) & -(uint64_t)indicator;
lobits ^= (lobits ^ new_lobits ) & -(uint64_t)indicator;
}
/*
* The product of the 128-bit value now in hibits:lobits with the
* 128-bit value we originally retrieved in the same variables
* will be in the vicinity of 2^191. So we'll take log2(R) to be
* 191, plus a multiple of BIGNUM_INT_BITS large enough to allow R
* to hold the combined sizes of n and d.
*/
size_t log2_R;
{
size_t max_log2_n = (n->nw + d->nw) * BIGNUM_INT_BITS;
log2_R = max_log2_n + 3;
log2_R -= size_t_min(191, log2_R);
log2_R = (log2_R + BIGNUM_INT_BITS - 1) & ~(BIGNUM_INT_BITS - 1);
log2_R += 191;
}
/* Number of words in a bignum capable of holding numbers the size
* of twice R. */
size_t rw = ((log2_R+2) + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS;
/*
* Now construct our full-sized starting reciprocal approximation.
*/
mp_int *r_approx = mp_make_sized(rw);
size_t output_bit_index;
{
/* Where in the input number did the input 128-bit value come from? */
size_t input_bit_index =
(hiword_index * BIGNUM_INT_BITS) - (128 - BIGNUM_INT_BITS);
/* So how far do we need to shift our 64-bit output, if the
* product of those two fixed-size values is 2^191 and we want
* to make it 2^log2_R instead? */
output_bit_index = log2_R - 191 - input_bit_index;
/* If we've done all that right, it should be a whole number
* of words. */
assert(output_bit_index % BIGNUM_INT_BITS == 0);
size_t output_word_index = output_bit_index / BIGNUM_INT_BITS;
mp_add_integer_into_shifted_by_words(
r_approx, r_approx, lobits, output_word_index);
mp_add_integer_into_shifted_by_words(
r_approx, r_approx, hibits,
output_word_index + 64 / BIGNUM_INT_BITS);
}
/*
* Make the constant 2*R, which we'll need in the iteration.
*/
mp_int *two_R = mp_make_sized(rw);
BignumInt top_word = (BignumInt)1 << ((log2_R+1) % BIGNUM_INT_BITS);
mp_add_integer_into_shifted_by_words(
two_R, two_R, top_word, (log2_R+1) / BIGNUM_INT_BITS);
/*
* Scratch space.
*/
mp_int *dr = mp_make_sized(rw + d->nw);
mp_int *diff = mp_make_sized(size_t_max(rw, dr->nw));
mp_int *product = mp_make_sized(rw + diff->nw);
size_t scratchsize = size_t_max(
mp_mul_scratchspace(dr->nw, r_approx->nw, d->nw),
mp_mul_scratchspace(product->nw, r_approx->nw, diff->nw));
mp_int *scratch = mp_make_sized(scratchsize);
mp_int product_shifted = mp_make_alias(
product, log2_R / BIGNUM_INT_BITS, product->nw);
/*
* Initial error estimate: the 32-bit output of recip_approx_32
* differs by less than 2048 (== 2^11) from the true top 32 bits
* of the reciprocal, so the relative error is at most 2^11
* divided by the 32-bit reciprocal, which at worst is 2^11/2^31 =
* 2^-20. So even in the worst case, we have 20 good bits of
* reciprocal to start with.
*/
size_t good_bits = 31 - 11;
size_t good_bits_needed = BIGNUM_INT_BITS * n->nw + 4; /* add a few */
/*
* Now do Newton-Raphson iterations until we have reason to think
* they're not converging any more.
*/
while (good_bits < good_bits_needed) {
/*
* Compute the next iterate.
*/
mp_mul_internal(dr, r_approx, d, *scratch);
mp_sub_into(diff, two_R, dr);
mp_mul_internal(product, r_approx, diff, *scratch);
mp_rshift_fixed_into(r_approx, &product_shifted,
log2_R % BIGNUM_INT_BITS);
/*
* Adjust the error estimate.
*/
good_bits = good_bits * 2 - 1;
}
mp_free(dr);
mp_free(diff);
mp_free(product);
mp_free(scratch);
/*
* Now we've got our reciprocal, we can compute the quotient, by
* multiplying in n and then shifting down by log2_R bits.
*/
mp_int *quotient_full = mp_mul(r_approx, n);
mp_int quotient_alias = mp_make_alias(
quotient_full, log2_R / BIGNUM_INT_BITS, quotient_full->nw);
mp_int *quotient = mp_make_sized(n->nw);
mp_rshift_fixed_into(quotient, "ient_alias, log2_R % BIGNUM_INT_BITS);
/*
* Next, compute the remainder.
*/
mp_int *remainder = mp_make_sized(d->nw);
mp_mul_into(remainder, quotient, d);
mp_sub_into(remainder, n, remainder);
/*
* Finally, two conditional subtractions to fix up any remaining
* rounding error. (I _think_ one should be enough, but this
* routine isn't time-critical enough to take chances.)
*/
unsigned q_correction = 0;
for (unsigned iter = 0; iter < 2; iter++) {
unsigned need_correction = mp_cmp_hs(remainder, d);
mp_cond_sub_into(remainder, remainder, d, need_correction);
q_correction += need_correction;
}
mp_add_integer_into(quotient, quotient, q_correction);
/*
* Now we should have a perfect answer, i.e. 0 <= r < d.
*/
assert(!mp_cmp_hs(remainder, d));
if (q_out)
mp_copy_into(q_out, quotient);
if (r_out)
mp_copy_into(r_out, remainder);
mp_free(r_approx);
mp_free(two_R);
mp_free(quotient_full);
mp_free(quotient);
mp_free(remainder);
}
mp_int *mp_div(mp_int *n, mp_int *d)
{
mp_int *q = mp_make_sized(n->nw);
mp_divmod_into(n, d, q, NULL);
return q;
}
mp_int *mp_mod(mp_int *n, mp_int *d)
{
mp_int *r = mp_make_sized(d->nw);
mp_divmod_into(n, d, NULL, r);
return r;
}
uint32_t mp_mod_known_integer(mp_int *x, uint32_t m)
{
uint64_t reciprocal = ((uint64_t)1 << 48) / m;
uint64_t accumulator = 0;
for (size_t i = mp_max_bytes(x); i-- > 0 ;) {
accumulator = 0x100 * accumulator + mp_get_byte(x, i);
/*
* Let A be the value in 'accumulator' at this point, and let
* R be the value it will have after we subtract quot*m below.
*
* Lemma 1: if A < 2^48, then R < 2m.
*
* Proof:
*
* By construction, we have 2^48/m - 1 < reciprocal <= 2^48/m.
* Multiplying that by the accumulator gives
*
* A/m * 2^48 - A < unshifted_quot <= A/m * 2^48
* i.e. 0 <= (A/m * 2^48) - unshifted_quot < A
* i.e. 0 <= A/m - unshifted_quot/2^48 < A/2^48
*
* So when we shift this quotient right by 48 bits, i.e. take
* the floor of (unshifted_quot/2^48), the value we take the
* floor of is at most A/2^48 less than the true rational
* value A/m that we _wanted_ to take the floor of.
*
* Provided A < 2^48, this is less than 1. So the quotient
* 'quot' that we've just produced is either the true quotient
* floor(A/m), or one less than it. Hence, the output value R
* is less than 2m. []
*
* Lemma 2: if A < 2^16 m, then the multiplication of
* accumulator*reciprocal does not overflow.
*
* Proof: as above, we have reciprocal <= 2^48/m. Multiplying
* by A gives unshifted_quot <= 2^48 * A / m < 2^48 * 2^16 =
* 2^64. []
*/
uint64_t unshifted_quot = accumulator * reciprocal;
uint64_t quot = unshifted_quot >> 48;
accumulator -= quot * m;
}
/*
* Theorem 1: accumulator < 2m at the end of every iteration of
* this loop.
*
* Proof: induction on the above loop.
*
* Base case: at the start of the first loop iteration, the
* accumulator is 0, which is certainly < 2m.
*
* Inductive step: in each loop iteration, we take a value at most
* 2m-1, multiply it by 2^8, and add another byte less than 2^8 to
* generate the input value A to the reduction process above. So
* we have A < 2m * 2^8 - 1. We know m < 2^32 (because it was
* passed in as a uint32_t), so A < 2^41, which is enough to allow
* us to apply Lemma 1, showing that the value of 'accumulator' at
* the end of the loop is still < 2m. []
*
* Corollary: we need at most one final subtraction of m to
* produce the canonical residue of x mod m, i.e. in the range
* [0,m).
*
* Theorem 2: no multiplication in the inner loop overflows.
*
* Proof: in Theorem 1 we established A < 2m * 2^8 - 1 in every
* iteration. That is less than m * 2^16, so Lemma 2 applies.
*
* The other multiplication, of quot * m, cannot overflow because
* quot is at most A/m, so quot*m <= A < 2^64. []
*/
uint32_t result = accumulator;
uint32_t reduced = result - m;
uint32_t select = -(reduced >> 31);
result = reduced ^ ((result ^ reduced) & select);
assert(result < m);
return result;
}
mp_int *mp_nthroot(mp_int *y, unsigned n, mp_int *remainder_out)
{
/*
* Allocate scratch space.
*/
mp_int **alloc, **powers, **newpowers, *scratch;
size_t nalloc = 2*(n+1)+1;
alloc = snewn(nalloc, mp_int *);
for (size_t i = 0; i < nalloc; i++)
alloc[i] = mp_make_sized(y->nw + 1);
powers = alloc;
newpowers = alloc + (n+1);
scratch = alloc[2*n+2];
/*
* We're computing the rounded-down nth root of y, i.e. the
* maximal x such that x^n <= y. We try to add 2^i to it for each
* possible value of i, starting from the largest one that might
* fit (i.e. such that 2^{n*i} fits in the size of y) downwards to
* i=0.
*
* We track all the smaller powers of x in the array 'powers'. In
* each iteration, if we update x, we update all of those values
* to match.
*/
mp_copy_integer_into(powers[0], 1);
for (size_t s = mp_max_bits(y) / n + 1; s-- > 0 ;) {
/*
* Let b = 2^s. We need to compute the powers (x+b)^i for each
* i, starting from our recorded values of x^i.
*/
for (size_t i = 0; i < n+1; i++) {
/*
* (x+b)^i = x^i
* + (i choose 1) x^{i-1} b
* + (i choose 2) x^{i-2} b^2
* + ...
* + b^i
*/
uint16_t binom = 1; /* coefficient of b^i */
mp_copy_into(newpowers[i], powers[i]);
for (size_t j = 0; j < i; j++) {
/* newpowers[i] += binom * powers[j] * 2^{(i-j)*s} */
mp_mul_integer_into(scratch, powers[j], binom);
mp_lshift_fixed_into(scratch, scratch, (i-j) * s);
mp_add_into(newpowers[i], newpowers[i], scratch);
uint32_t binom_mul = binom;
binom_mul *= (i-j);
binom_mul /= (j+1);
assert(binom_mul < 0x10000);
binom = binom_mul;
}
}
/*
* Now, is the new value of x^n still <= y? If so, update.
*/
unsigned newbit = mp_cmp_hs(y, newpowers[n]);
for (size_t i = 0; i < n+1; i++)
mp_select_into(powers[i], powers[i], newpowers[i], newbit);
}
if (remainder_out)
mp_sub_into(remainder_out, y, powers[n]);
mp_int *root = mp_new(mp_max_bits(y) / n);
mp_copy_into(root, powers[1]);
for (size_t i = 0; i < nalloc; i++)
mp_free(alloc[i]);
sfree(alloc);
return root;
}
mp_int *mp_modmul(mp_int *x, mp_int *y, mp_int *modulus)
{
mp_int *product = mp_mul(x, y);
mp_int *reduced = mp_mod(product, modulus);
mp_free(product);
return reduced;
}
mp_int *mp_modadd(mp_int *x, mp_int *y, mp_int *modulus)
{
mp_int *sum = mp_add(x, y);
mp_int *reduced = mp_mod(sum, modulus);
mp_free(sum);
return reduced;
}
mp_int *mp_modsub(mp_int *x, mp_int *y, mp_int *modulus)
{
mp_int *diff = mp_make_sized(size_t_max(x->nw, y->nw));
mp_sub_into(diff, x, y);
unsigned negate = mp_cmp_hs(y, x);
mp_cond_negate(diff, diff, negate);
mp_int *residue = mp_mod(diff, modulus);
mp_cond_negate(residue, residue, negate);
/* If we've just negated the residue, then it will be < 0 and need
* the modulus adding to it to make it positive - *except* if the
* residue was zero when we negated it. */
unsigned make_positive = negate & ~mp_eq_integer(residue, 0);
mp_cond_add_into(residue, residue, modulus, make_positive);
mp_free(diff);
return residue;
}
static mp_int *mp_modadd_in_range(mp_int *x, mp_int *y, mp_int *modulus)
{
mp_int *sum = mp_make_sized(modulus->nw);
unsigned carry = mp_add_into_internal(sum, x, y);
mp_cond_sub_into(sum, sum, modulus, carry | mp_cmp_hs(sum, modulus));
return sum;
}
static mp_int *mp_modsub_in_range(mp_int *x, mp_int *y, mp_int *modulus)
{
mp_int *diff = mp_make_sized(modulus->nw);
mp_sub_into(diff, x, y);
mp_cond_add_into(diff, diff, modulus, 1 ^ mp_cmp_hs(x, y));
return diff;
}
mp_int *monty_add(MontyContext *mc, mp_int *x, mp_int *y)
{
return mp_modadd_in_range(x, y, mc->m);
}
mp_int *monty_sub(MontyContext *mc, mp_int *x, mp_int *y)
{
return mp_modsub_in_range(x, y, mc->m);
}
void mp_min_into(mp_int *r, mp_int *x, mp_int *y)
{
mp_select_into(r, x, y, mp_cmp_hs(x, y));
}
void mp_max_into(mp_int *r, mp_int *x, mp_int *y)
{
mp_select_into(r, y, x, mp_cmp_hs(x, y));
}
mp_int *mp_min(mp_int *x, mp_int *y)
{
mp_int *r = mp_make_sized(size_t_min(x->nw, y->nw));
mp_min_into(r, x, y);
return r;
}
mp_int *mp_max(mp_int *x, mp_int *y)
{
mp_int *r = mp_make_sized(size_t_max(x->nw, y->nw));
mp_max_into(r, x, y);
return r;
}
mp_int *mp_power_2(size_t power)
{
mp_int *x = mp_new(power + 1);
mp_set_bit(x, power, 1);
return x;
}
struct ModsqrtContext {
mp_int *p; /* the prime */
MontyContext *mc; /* for doing arithmetic mod p */
/* Decompose p-1 as 2^e k, for positive integer e and odd k */
size_t e;
mp_int *k;
mp_int *km1o2; /* (k-1)/2 */
/* The user-provided value z which is not a quadratic residue mod
* p, and its kth power. Both in Montgomery form. */
mp_int *z, *zk;
};
ModsqrtContext *modsqrt_new(mp_int *p, mp_int *any_nonsquare_mod_p)
{
ModsqrtContext *sc = snew(ModsqrtContext);
memset(sc, 0, sizeof(ModsqrtContext));
sc->p = mp_copy(p);
sc->mc = monty_new(sc->p);
sc->z = monty_import(sc->mc, any_nonsquare_mod_p);
/* Find the lowest set bit in p-1. Since this routine expects p to
* be non-secret (typically a well-known standard elliptic curve
* parameter), for once we don't need clever bit tricks. */
for (sc->e = 1; sc->e < BIGNUM_INT_BITS * p->nw; sc->e++)
if (mp_get_bit(p, sc->e))
break;
sc->k = mp_rshift_fixed(p, sc->e);
sc->km1o2 = mp_rshift_fixed(sc->k, 1);
/* Leave zk to be filled in lazily, since it's more expensive to
* compute. If this context turns out never to be needed, we can
* save the bulk of the setup time this way. */
return sc;
}
static void modsqrt_lazy_setup(ModsqrtContext *sc)
{
if (!sc->zk)
sc->zk = monty_pow(sc->mc, sc->z, sc->k);
}
void modsqrt_free(ModsqrtContext *sc)
{
monty_free(sc->mc);
mp_free(sc->p);
mp_free(sc->z);
mp_free(sc->k);
mp_free(sc->km1o2);
if (sc->zk)
mp_free(sc->zk);
sfree(sc);
}
mp_int *mp_modsqrt(ModsqrtContext *sc, mp_int *x, unsigned *success)
{
mp_int *mx = monty_import(sc->mc, x);
mp_int *mroot = monty_modsqrt(sc, mx, success);
mp_free(mx);
mp_int *root = monty_export(sc->mc, mroot);
mp_free(mroot);
return root;
}
/*
* Modular square root, using an algorithm more or less similar to
* Tonelli-Shanks but adapted for constant time.
*
* The basic idea is to write p-1 = k 2^e, where k is odd and e > 0.
* Then the multiplicative group mod p (call it G) has a sequence of
* e+1 nested subgroups G = G_0 > G_1 > G_2 > ... > G_e, where each
* G_i is exactly half the size of G_{i-1} and consists of all the
* squares of elements in G_{i-1}. So the innermost group G_e has
* order k, which is odd, and hence within that group you can take a
* square root by raising to the power (k+1)/2.
*
* Our strategy is to iterate over these groups one by one and make
* sure the number x we're trying to take the square root of is inside
* each one, by adjusting it if it isn't.
*
* Suppose g is a primitive root of p, i.e. a generator of G_0. (We
* don't actually need to know what g _is_; we just imagine it for the
* sake of understanding.) Then G_i consists of precisely the (2^i)th
* powers of g, and hence, you can tell if a number is in G_i if
* raising it to the power k 2^{e-i} gives 1. So the conceptual
* algorithm goes: for each i, test whether x is in G_i by that
* method. If it isn't, then the previous iteration ensured it's in
* G_{i-1}, so it will be an odd power of g^{2^{i-1}}, and hence
* multiplying by any other odd power of g^{2^{i-1}} will give x' in
* G_i. And we have one of those, because our non-square z is an odd
* power of g, so z^{2^{i-1}} is an odd power of g^{2^{i-1}}.
*
* (There's a special case in the very first iteration, where we don't
* have a G_{i-1}. If it turns out that x is not even in G_1, that
* means it's not a square, so we set *success to 0. We still run the
* rest of the algorithm anyway, for the sake of constant time, but we
* don't give a hoot what it returns.)
*
* When we get to the end and have x in G_e, then we can take its
* square root by raising to (k+1)/2. But of course that's not the
* square root of the original input - it's only the square root of
* the adjusted version we produced during the algorithm. To get the
* true output answer we also have to multiply by a power of z,
* namely, z to the power of _half_ whatever we've been multiplying in
* as we go along. (The power of z we multiplied in must have been
* even, because the case in which we would have multiplied in an odd
* power of z is the i=0 case, in which we instead set the failure
* flag.)
*
* The code below is an optimised version of that basic idea, in which
* we _start_ by computing x^k so as to be able to test membership in
* G_i by only a few squarings rather than a full from-scratch modpow
* every time; we also start by computing our candidate output value
* x^{(k+1)/2}. So when the above description says 'adjust x by z^i'
* for some i, we have to adjust our running values of x^k and
* x^{(k+1)/2} by z^{ik} and z^{ik/2} respectively (the latter is safe
* because, as above, i is always even). And it turns out that we
* don't actually have to store the adjusted version of x itself at
* all - we _only_ keep those two powers of it.
*/
mp_int *monty_modsqrt(ModsqrtContext *sc, mp_int *x, unsigned *success)
{
modsqrt_lazy_setup(sc);
mp_int *scratch_to_free = mp_make_sized(3 * sc->mc->rw);
mp_int scratch = *scratch_to_free;
/*
* Compute toret = x^{(k+1)/2}, our starting point for the output
* square root, and also xk = x^k which we'll use as we go along
* for knowing when to apply correction factors. We do this by
* first computing x^{(k-1)/2}, then multiplying it by x, then
* multiplying the two together.
*/
mp_int *toret = monty_pow(sc->mc, x, sc->km1o2);
mp_int xk = mp_alloc_from_scratch(&scratch, sc->mc->rw);
mp_copy_into(&xk, toret);
monty_mul_into(sc->mc, toret, toret, x);
monty_mul_into(sc->mc, &xk, toret, &xk);
mp_int tmp = mp_alloc_from_scratch(&scratch, sc->mc->rw);
mp_int power_of_zk = mp_alloc_from_scratch(&scratch, sc->mc->rw);
mp_copy_into(&power_of_zk, sc->zk);
for (size_t i = 0; i < sc->e; i++) {
mp_copy_into(&tmp, &xk);
for (size_t j = i+1; j < sc->e; j++)
monty_mul_into(sc->mc, &tmp, &tmp, &tmp);
unsigned eq1 = mp_cmp_eq(&tmp, monty_identity(sc->mc));
if (i == 0) {
/* One special case: if x=0, then no power of x will ever
* equal 1, but we should still report success on the
* grounds that 0 does have a square root mod p. */
*success = eq1 | mp_eq_integer(x, 0);
} else {
monty_mul_into(sc->mc, &tmp, toret, &power_of_zk);
mp_select_into(toret, &tmp, toret, eq1);
monty_mul_into(sc->mc, &power_of_zk,
&power_of_zk, &power_of_zk);
monty_mul_into(sc->mc, &tmp, &xk, &power_of_zk);
mp_select_into(&xk, &tmp, &xk, eq1);
}
}
mp_free(scratch_to_free);
return toret;
}
mp_int *mp_random_bits_fn(size_t bits, random_read_fn_t random_read)
{
size_t bytes = (bits + 7) / 8;
uint8_t *randbuf = snewn(bytes, uint8_t);
random_read(randbuf, bytes);
if (bytes)
randbuf[0] &= (2 << ((bits-1) & 7)) - 1;
mp_int *toret = mp_from_bytes_be(make_ptrlen(randbuf, bytes));
smemclr(randbuf, bytes);
sfree(randbuf);
return toret;
}
mp_int *mp_random_upto_fn(mp_int *limit, random_read_fn_t rf)
{
/*
* It would be nice to generate our random numbers in such a way
* as to make every possible outcome literally equiprobable. But
* we can't do that in constant time, so we have to go for a very
* close approximation instead. I'm going to take the view that a
* factor of (1+2^-128) between the probabilities of two outcomes
* is acceptable on the grounds that you'd have to examine so many
* outputs to even detect it.
*/
mp_int *unreduced = mp_random_bits_fn(mp_max_bits(limit) + 128, rf);
mp_int *reduced = mp_mod(unreduced, limit);
mp_free(unreduced);
return reduced;
}
mp_int *mp_random_in_range_fn(mp_int *lo, mp_int *hi, random_read_fn_t rf)
{
mp_int *n_outcomes = mp_sub(hi, lo);
mp_int *addend = mp_random_upto_fn(n_outcomes, rf);
mp_int *result = mp_make_sized(hi->nw);
mp_add_into(result, addend, lo);
mp_free(addend);
mp_free(n_outcomes);
return result;
}
|