1 /*  2 *******************************************************************************  3 \file pri.c  4 \brief Prime numbers  5 \project bee2 [cryptographic library]  6 \author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]  7 \created 2012.08.13  8 \version 2016.07.15  9 \license This program is released under the GNU General Public License  10 version 3. See Copyright Notices in bee2/info.h.  11 *******************************************************************************  12 */  13 14 #include "bee2/core/mem.h"  15 #include "bee2/core/prng.h"  16 #include "bee2/core/util.h"  17 #include "bee2/core/word.h"  18 #include "bee2/math/pri.h"  19 #include "bee2/math/ww.h"  20 #include "bee2/math/zm.h"  21 22 /*  23 *******************************************************************************  24 Факторная база: первые 1024 нечетных простых  25 26 Построены в Mathematica:  27 \code  28  n = 1024;  29  Table[Prime[i], {i, 2, n + 1}]  30 \endcode  31 *******************************************************************************  32 */  33 34 static const word _base[] =  35 {  36  3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,  37  73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,  38  157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,  39  239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,  40  331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,  41  421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,  42  509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,  43  613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,  44  709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,  45  821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,  46  919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,  47  1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,  48  1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,  49  1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,  50  1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361,  51  1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,  52  1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531,  53  1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,  54  1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699,  55  1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,  56  1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,  57  1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,  58  1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083,  59  2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,  60  2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273,  61  2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,  62  2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441,  63  2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,  64  2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663,  65  2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,  66  2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,  67  2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,  68  2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023,  69  3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,  70  3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251,  71  3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,  72  3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,  73  3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,  74  3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617,  75  3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,  76  3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,  77  3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,  78  3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013,  79  4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,  80  4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219,  81  4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,  82  4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423,  83  4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,  84  4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639,  85  4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,  86  4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,  87  4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,  88  4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023,  89  5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147,  90  5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261,  91  5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,  92  5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,  93  5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563,  94  5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659,  95  5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779,  96  5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,  97  5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981,  98  5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089,  99  6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199,  100  6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287,  101  6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,  102  6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491,  103  6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607,  104  6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709,  105  6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827,  106  6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,  107  6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013,  108  7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129,  109  7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243,  110  7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369,  111  7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,  112  7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577,  113  7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681,  114  7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,  115  7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901,  116  7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017,  117  8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123,  118  8147, 8161, 8167,  119 };  120 121 1 size_t priBaseSize()  122 {  123 1  return COUNT_OF(_base);  124 }  125 126 1 word priBasePrime(size_t i)  127 {  128 1  ASSERT(i < priBaseSize());  129 1  return _base[i];  130 }  131 132 /*  133 *******************************************************************************  134 Произведения последовательных простых из факторной базы: множители  135 накапливаются до тех пор, пока произведение помещается в машинное слово.  136 Тривиальные произведения из одного множителя не фиксируются.  137 138 Построены в Mathematica:  139 \code  140  n = 1024;  141  w = B_PER_W;  142  prods = {};  143  For[i = 2, i <= n + 1,  144  For[p = 1; l = 0,  145  i + l <= n + 1 && p * Prime[i + ++l] < 2^w,  146  p *= Prime[i + l]];  147  i += l;  148  If[l > 1, AppendTo[prods, {p, l}]  149  ]  150  prods  151 \endcode  152 *******************************************************************************  153 */  154 155 typedef struct pri_prod_t  156 {  157  word prod; /*< произведение простых из факторной базы */  158  size_t num; /*< количество множителей */  159 } pri_prod_t;  160 161 static const pri_prod_t _prods[] =  162 {  163 #if (B_PER_W == 16)  164  {15015u, 5}, {7429u, 3}, {33263u, 3}, {1763u, 2}, {2491u, 2},  165  {3599u, 2}, {4757u, 2}, {5767u, 2}, {7387u, 2}, {9797u, 2}, {11021u, 2},  166  {12317u, 2}, {16637u, 2}, {19043u, 2}, {22499u, 2}, {25591u, 2},  167  {28891u, 2}, {32399u, 2}, {36863u, 2}, {39203u, 2}, {47053u, 2},  168  {51983u, 2}, {55687u, 2}, {60491u, 2},  169 #elif (B_PER_W == 32)  170  {3234846615u, 9}, {95041567u, 5}, {907383479u, 5}, {4132280413u, 5},  171  {121330189u, 4}, {257557397u, 4}, {490995677u, 4}, {842952707u, 4},  172  {1314423991u, 4}, {2125525169u, 4}, {3073309843u, 4}, {16965341u, 3},  173  {20193023u, 3}, {23300239u, 3}, {29884301u, 3}, {35360399u, 3},  174  {42749359u, 3}, {49143869u, 3}, {56466073u, 3}, {65111573u, 3},  175  {76027969u, 3}, {84208541u, 3}, {94593973u, 3}, {103569859u, 3},  176  {119319383u, 3}, {133390067u, 3}, {154769821u, 3}, {178433279u, 3},  177  {193397129u, 3}, {213479407u, 3}, {229580147u, 3}, {250367549u, 3},  178  {271661713u, 3}, {293158127u, 3}, {319512181u, 3}, {357349471u, 3},  179  {393806449u, 3}, {422400701u, 3}, {452366557u, 3}, {507436351u, 3},  180  {547978913u, 3}, {575204137u, 3}, {627947039u, 3}, {666785731u, 3},  181  {710381447u, 3}, {777767161u, 3}, {834985999u, 3}, {894826021u, 3},  182  {951747481u, 3}, {1019050649u, 3}, {1072651369u, 3}, {1125878063u, 3},  183  {1185362993u, 3}, {1267745273u, 3}, {1322520163u, 3}, {1391119619u, 3},  184  {1498299287u, 3}, {1608372013u, 3}, {1700725291u, 3}, {1805418283u, 3},  185  {1871456063u, 3}, {2008071007u, 3}, {2115193573u, 3}, {2178429527u, 3},  186  {2246284699u, 3}, {2385788087u, 3}, {2591986471u, 3}, {2805004793u, 3},  187  {2922149239u, 3}, {3021320083u, 3}, {3118412617u, 3}, {3265932301u, 3},  188  {3332392423u, 3}, {3523218343u, 3}, {3711836171u, 3}, {3837879163u, 3},  189  {3991792529u, 3}, {4139646463u, 3}, {4233155587u, 3}, {2663399u, 2},  190  {2755591u, 2}, {2782223u, 2}, {2873021u, 2}, {2903591u, 2}, {2965283u, 2},  191  {3017153u, 2}, {3062491u, 2}, {3125743u, 2}, {3186221u, 2}, {3221989u, 2},  192  {3301453u, 2}, {3381857u, 2}, {3474487u, 2}, {3504383u, 2}, {3526883u, 2},  193  {3590989u, 2}, {3648091u, 2}, {3732623u, 2}, {3802499u, 2}, {3904567u, 2},  194  {3960091u, 2}, {3992003u, 2}, {4028033u, 2}, {4088459u, 2}, {4137131u, 2},  195  {4235339u, 2}, {4305589u, 2}, {4347221u, 2}, {4384811u, 2}, {4460543u, 2},  196  {4536899u, 2}, {4575317u, 2}, {4613879u, 2}, {4708819u, 2}, {4862021u, 2},  197  {4915073u, 2}, {5008643u, 2}, {5048993u, 2}, {5143823u, 2}, {5184713u, 2},  198  {5244091u, 2}, {5303773u, 2}, {5391563u, 2}, {5475599u, 2}, {5517797u, 2},  199  {5588447u, 2}, {5659637u, 2}, {5692987u, 2}, {5740807u, 2}, {5827387u, 2},  200  {5904851u, 2}, {5973127u, 2}, {6066353u, 2}, {6125621u, 2}, {6310063u, 2},  201  {6426209u, 2}, {6482107u, 2}, {6522907u, 2}, {6682189u, 2}, {6765137u, 2},  202  {6859157u, 2}, {6969551u, 2}, {7064963u, 2}, {7112873u, 2}, {7182391u, 2},  203  {7225343u, 2}, {7268407u, 2}, {7338677u, 2}, {7376647u, 2}, {7452899u, 2},  204  {7535009u, 2}, {7617551u, 2}, {7745053u, 2}, {7806427u, 2}, {7851203u, 2},  205  {7986227u, 2}, {8065591u, 2}, {8145307u, 2}, {8236819u, 2}, {8363639u, 2},  206  {8444827u, 2}, {8538059u, 2}, {8678867u, 2}, {8761591u, 2}, {8820899u, 2},  207  {8999999u, 2}, {9090209u, 2}, {9180851u, 2}, {9272009u, 2}, {9388087u, 2},  208  {9492557u, 2}, {9603701u, 2}, {9734399u, 2}, {9922331u, 2}, {10036223u, 2},  209  {10137847u, 2}, {10220773u, 2}, {10323353u, 2}, {10400609u, 2},  210  {10575503u, 2}, {10614563u, 2}, {10791029u, 2}, {10916407u, 2},  211  {10995847u, 2}, {11062267u, 2}, {11135533u, 2}, {11242573u, 2},  212  {11329931u, 2}, {11431097u, 2}, {11553137u, 2}, {11716829u, 2},  213  {11923193u, 2}, {11985443u, 2}, {12027023u, 2}, {12215009u, 2},  214  {12348187u, 2}, {12446783u, 2}, {12503287u, 2}, {12559927u, 2},  215  {12659363u, 2}, {12787751u, 2}, {12873719u, 2}, {13032091u, 2},  216  {13104391u, 2}, {13205947u, 2}, {13329737u, 2}, {13483583u, 2},  217  {13571807u, 2}, {13682597u, 2}, {13793771u, 2}, {13912891u, 2},  218  {14062379u, 2}, {14197823u, 2}, {14333747u, 2}, {14439991u, 2},  219  {14607683u, 2}, {14745551u, 2}, {14837903u, 2}, {14976851u, 2},  220  {15093209u, 2}, {15280277u, 2}, {15350723u, 2}, {15413467u, 2},  221  {15499933u, 2}, {15657749u, 2}, {15959989u, 2}, {16040021u, 2},  222  {16128247u, 2}, {16192567u, 2}, {16402499u, 2}, {16524161u, 2},  223  {16687189u, 2}, {16777207u, 2}, {16966097u, 2}, {17065157u, 2},  224  {17189267u, 2}, {17288963u, 2}, {17547577u, 2}, {17757787u, 2},  225  {17842151u, 2}, {17943671u, 2}, {18045479u, 2}, {18147599u, 2},  226  {18249983u, 2}, {18369787u, 2}, {18593119u, 2}, {18818243u, 2},  227  {18948593u, 2}, {19079399u, 2}, {19307227u, 2}, {19492189u, 2},  228  {19642543u, 2}, {19793597u, 2}, {19891591u, 2}, {20088323u, 2},  229  {20249951u, 2}, {20385221u, 2}, {20439437u, 2}, {20684303u, 2},  230  {20830087u, 2}, {21040553u, 2}, {21159991u, 2}, {21427577u, 2},  231  {21538877u, 2}, {21622499u, 2}, {21715591u, 2}, {21864967u, 2},  232  {22061773u, 2}, {22297283u, 2}, {22382357u, 2}, {22610009u, 2},  233  {22896221u, 2}, {22953677u, 2}, {23039999u, 2}, {23184221u, 2},  234  {23483491u, 2}, {23755867u, 2}, {23970767u, 2}, {24147371u, 2},  235  {24324623u, 2}, {24403591u, 2}, {24542107u, 2}, {24681023u, 2},  236  {24800351u, 2}, {24960007u, 2}, {25060027u, 2}, {25160231u, 2},  237  {25310897u, 2}, {25553009u, 2}, {25796237u, 2}, {25938613u, 2},  238  {26050807u, 2}, {26173447u, 2}, {26522491u, 2}, {26718557u, 2},  239  {26873831u, 2}, {27071173u, 2}, {27342437u, 2}, {27405221u, 2},  240  {27741253u, 2}, {27878399u, 2}, {28089991u, 2}, {28259807u, 2},  241  {28515551u, 2}, {28793731u, 2}, {29052091u, 2}, {29192393u, 2},  242  {29322221u, 2}, {29430589u, 2}, {29582717u, 2}, {29658907u, 2},  243  {29964667u, 2}, {30041357u, 2}, {30272003u, 2}, {30393133u, 2},  244  {30514567u, 2}, {30735767u, 2}, {30980347u, 2}, {31102913u, 2},  245  {31438193u, 2}, {31809599u, 2}, {31911197u, 2}, {31979021u, 2},  246  {32080871u, 2}, {32330587u, 2}, {32455793u, 2}, {32649787u, 2},  247  {32936117u, 2}, {33016507u, 2}, {33419957u, 2}, {33593591u, 2},  248  {33756091u, 2}, {33918967u, 2}, {34117277u, 2}, {34222499u, 2},  249  {34327877u, 2}, {34433423u, 2}, {34574399u, 2}, {34809991u, 2},  250  {35105621u, 2}, {35354867u, 2}, {35808247u, 2}, {36108077u, 2},  251  {36397073u, 2}, {36542021u, 2}, {36723551u, 2}, {36917767u, 2},  252  {37088099u, 2}, {37295413u, 2}, {37527851u, 2}, {37675019u, 2},  253  {37908613u, 2}, {38254081u, 2}, {38452397u, 2}, {38613787u, 2},  254  {38750609u, 2}, {39087479u, 2}, {39262747u, 2}, {39363067u, 2},  255  {39601813u, 2}, {39765611u, 2}, {39942391u, 2}, {40106873u, 2},  256  {40297079u, 2}, {40449599u, 2}, {40576891u, 2}, {40755431u, 2},  257  {41075137u, 2}, {41447723u, 2}, {41731519u, 2}, {41951513u, 2},  258  {42327811u, 2}, {42745363u, 2}, {42928703u, 2}, {43112347u, 2},  259  {43217467u, 2}, {43428019u, 2}, {43731733u, 2}, {44155961u, 2},  260  {44355599u, 2}, {44568967u, 2}, {44756099u, 2}, {44916803u, 2},  261  {45077771u, 2}, {45360221u, 2}, {45724643u, 2}, {45968399u, 2},  262  {46131263u, 2}, {46416869u, 2}, {46621583u, 2}, {46744553u, 2},  263  {47059591u, 2}, {47196899u, 2}, {47485817u, 2}, {47734277u, 2},  264  {48052399u, 2}, {48358091u, 2}, {48497287u, 2}, {48636667u, 2},  265  {48818153u, 2}, {48985997u, 2}, {49224247u, 2}, {49463053u, 2},  266  {49702451u, 2}, {50041451u, 2}, {50495227u, 2}, {50751367u, 2},  267  {50979479u, 2}, {51380143u, 2}, {51696091u, 2}, {51969677u, 2},  268  {52070647u, 2}, {52316273u, 2}, {52490021u, 2}, {52823599u, 2},  269  {53319179u, 2}, {53509189u, 2}, {53758223u, 2}, {54022499u, 2},  270  {54479017u, 2}, {54967387u, 2}, {55383283u, 2}, {55621763u, 2},  271  {55935437u, 2}, {56070143u, 2}, {56294993u, 2}, {56550391u, 2},  272  {56746073u, 2}, {56911927u, 2}, {57062891u, 2}, {57259453u, 2},  273  {57456391u, 2}, {57608099u, 2}, {57836021u, 2}, {58216819u, 2},  274  {58461307u, 2}, {58844237u, 2}, {59043847u, 2}, {59213009u, 2},  275  {59444051u, 2}, {59675621u, 2}, {60015973u, 2}, {60186563u, 2},  276  {60699677u, 2}, {61152391u, 2}, {61387189u, 2}, {61779551u, 2},  277  {62015621u, 2}, {62110157u, 2}, {62473207u, 2}, {62773913u, 2},  278  {62964221u, 2}, {63202499u, 2}, {63648259u, 2}, {64160099u, 2},  279  {64448663u, 2}, {64899127u, 2}, {65205589u, 2}, {65415743u, 2},  280  {65561393u, 2}, {65836987u, 2}, {66178081u, 2}, {66650887u, 2},  281 #else  282  {16294579238595022365u, 15}, {7145393598349078859u, 10},  283  {6408001374760705163u, 9}, {690862709424854779u, 8},  284  {4312024209383942993u, 8}, {71235931512604841u, 7},  285  {192878245514479103u, 7}, {542676746453092519u, 7},  286  {1230544604996048471u, 7}, {2618501576975440661u, 7},  287  {4771180125133726009u, 7}, {9247077179230889629u, 7},  288  {32156968791364271u, 6}, {46627620659631719u, 6}, {64265583549260393u, 6},  289  {88516552714582021u, 6}, {131585967012906751u, 6},  290  {182675399263485151u, 6}, {261171077386532413u, 6},  291  {346060227726080771u, 6}, {448604664249794309u, 6},  292  {621993868801161359u, 6}, {813835565706097817u, 6},  293  {1050677302683430441u, 6}, {1294398862104002783u, 6},  294  {1615816556891330179u, 6}, {1993926996710486603u, 6},  295  {2626074105497143999u, 6}, {3280430033433832817u, 6},  296  {4076110663011485663u, 6}, {4782075577404875363u, 6},  297  {5906302864496324923u, 6}, {7899206880638488339u, 6},  298  {9178333502078117453u, 6}, {10680076322389870367u, 6},  299  {12622882367374918799u, 6}, {14897925470078818423u, 6},  300  {17264316336968551717u, 6}, {11896905306684389u, 5},  301  {13580761294555417u, 5}, {15289931661301991u, 5},  302  {17067874133764579u, 5}, {19008757261780379u, 5},  303  {21984658219193689u, 5}, {23721541361298551u, 5},  304  {26539432378378657u, 5}, {30167221680049747u, 5},  305  {32433198277139683u, 5}, {35517402656173043u, 5},  306  {39100537712055041u, 5}, {42477532426853543u, 5},  307  {45618621452253523u, 5}, {52071972962579407u, 5},  308  {57329264013213233u, 5}, {61692083285823527u, 5},  309  {66885169838978461u, 5}, {72186879569637319u, 5},  310  {77103033998665567u, 5}, {82549234838454463u, 5},  311  {89609394623390063u, 5}, {100441814079170659u, 5},  312  {109045745121501371u, 5}, {120230527473437819u, 5},  313  {131125107904515419u, 5}, {138612182127286823u, 5},  314  {144712752835963307u, 5}, {152692680370726429u, 5},  315  {164664356404541573u, 5}, {175376065798883557u, 5},  316  {187958301132741257u, 5}, {203342285718459187u, 5},  317  {219115706321995421u, 5}, {235226887496676263u, 5},  318  {253789253193479219u, 5}, {271717583502831491u, 5},  319  {293266389497362763u, 5}, {321821627692439603u, 5},  320  {339856237957830049u, 5}, {362469273063260281u, 5},  321  {390268963330916339u, 5}, {408848490015359209u, 5},  322  {429644565036857699u, 5}, {458755816747679897u, 5},  323  {495450768525623033u, 5}, {523240424009891327u, 5},  324  {551070603968128061u, 5}, {574205321266688311u, 5},  325  {606829434176923693u, 5}, {637763212653336997u, 5},  326  {676538378976146257u, 5}, {710263471119657661u, 5},  327  {754496879875465343u, 5}, {800075738315885429u, 5},  328  {845197573085733239u, 5}, {894146362391888161u, 5},  329  {930105507041885771u, 5}, {985345849616172623u, 5},  330  {1040222328124784927u, 5}, {1091468150538871153u, 5},  331  {1150933747479716653u, 5}, {1210604027868555713u, 5},  332  {1277530693373553361u, 5}, {1350088100087645657u, 5},  333  {1398676120233167591u, 5}, {1459450139327525269u, 5},  334  {1555755169940697937u, 5}, {1645735334920325819u, 5},  335  {1732866357938791147u, 5}, {1815492864312158099u, 5},  336  {1894564116319543619u, 5}, {1993720023757886939u, 5},  337  {2103356633892712673u, 5}, {2180099035358103487u, 5},  338  {2277315690161244011u, 5}, {2390146379836558999u, 5},  339  {2522126262040806983u, 5}, {2613887383383648311u, 5},  340  {2795412145600606001u, 5}, {2919958494381348367u, 5},  341  {3012266980379247553u, 5}, {3119360766859522543u, 5},  342  {3216618305468232557u, 5}, {3385071039891962579u, 5},  343  {3509427475807939163u, 5}, {3700008514672760651u, 5},  344  {3873423910591589033u, 5}, {4050200067084600439u, 5},  345  {4233429923647833421u, 5}, {4472862244562412787u, 5},  346  {4638587045132438407u, 5}, {4765110805097342489u, 5},  347  {4951886102290887619u, 5}, {5103665412856065733u, 5},  348  {5306636943410213377u, 5}, {5581202660702121667u, 5},  349  {5774946339890457283u, 5}, {5948565823654343479u, 5},  350  {6175776426345604697u, 5}, {6454381412132929663u, 5},  351  {6685489970462824483u, 5}, {6864273057227912189u, 5},  352  {7020466399135670969u, 5}, {7326535375987923521u, 5},  353  {7795297680723533551u, 5}, {8101368883577379127u, 5},  354  {8353548973662446233u, 5}, {8642941459163335097u, 5},  355  {8989536585105548947u, 5}, {9281597047973093449u, 5},  356  {9623989560555822323u, 5}, {9884958654427267811u, 5},  357  {10161260209201654649u, 5}, {10427299413974952277u, 5},  358  {10759022378261015069u, 5}, {11290266633494854903u, 5},  359  {11852839900193264543u, 5}, {12209591760462366097u, 5},  360  {12604874462407914499u, 5}, {13152204116524238149u, 5},  361  {13487108616207513373u, 5}, {13935742926215786057u, 5},  362  {14426307807825845411u, 5}, {14869398407236523377u, 5},  363  {15287602532539390847u, 5}, {15824557271701228177u, 5},  364  {16348641675671844607u, 5}, {16684838939207505557u, 5},  365  {17148166373867218913u, 5}, {17832011695956938489u, 5},  366  {2587278248197793u, 4}, {2656152548121013u, 4}, {2706094705771019u, 4},  367  {2746082268411733u, 4}, {2816510930505221u, 4}, {2876558914811147u, 4},  368  {2943092641403483u, 4}, {3044274349991521u, 4}, {3111227620115431u, 4},  369  {3156468307693999u, 4}, {3209012615864543u, 4}, {3247559087000957u, 4},  370  {3289921520014123u, 4}, {3331823223534079u, 4}, {3403431328122433u, 4},  371  {3474390126259739u, 4}, {3519861126859459u, 4}, {3581490458694233u, 4},  372  {3653304933840151u, 4}, {3753973384118899u, 4}, {3831297220366171u, 4},  373  {3880220695063499u, 4}, {3952510531166773u, 4}, {4022729025799241u, 4},  374  {4135032598497637u, 4}, {4231785801620803u, 4}, {4288747235209999u, 4},  375  {4356965458481947u, 4}, {66650887u, 2},  376 #endif  377 };  378 379 1 void priBaseMod(word mods[], const word a[], size_t n, size_t count)  380 {  381  size_t i, j;  382  // pre  383 1  ASSERT(wwIsValid(a, n));  384 1  ASSERT(count <= priBaseSize());  385 1  ASSERT(wwIsValid(mods, count));  386  // пробегаем произведения простых из факторной базы  387 1  for (i = j = 0; i < count && j < COUNT_OF(_prods); ++j)  388  {  389 1  size_t num = _prods[j].num;  390  // t <- a mod _base[i] * ... * _base[i + num - 1]  391 1  word t = zzModW(a, n, _prods[j].prod);  392  // mods[i] <- t mod _base[i]  393 1  while (num-- && i < count)  394 1  mods[i] = t % _base[i], ++i;  395  }  396  // пробегаем оставшиеся простые из факторной базы  397 1  for (; i < count; ++i)  398 0  mods[i] = zzModW(a, n, _base[i]);  399 }  400 401 /*  402 *******************************************************************************  403 Использование факторной базы  404 405 \todo Алгоритм Бернштейна выделения гладкой части  406 [http://cr.yp.to/factorization/smoothparts-20040510.pdf]:  407  z <- p1 p2 .... pm \mod a  408  y <- z^{2^e} \mod a (e --- минимальное натуральное т.ч. 3^{2^e} > a)  409  return \gcd(x, y)  410 *******************************************************************************  411 */  412 413 1 bool_t priIsSieved(const word a[], size_t n, size_t base_count, void* stack)  414 {  415  // переменные в stack  416  word* mods;  417  // pre  418 1  ASSERT(base_count <= priBaseSize());  419  // четное?  420 1  n = wwWordSize(a, n);  421 1  if (zzIsEven(a, n))  422 0  return FALSE;  423  // малое a?  424 1  if (n == 1)  425  // при необходимости скоррректировать факторную базу  426 0  while (base_count > 0 && priBasePrime(base_count - 1) > a[0])  427 0  --base_count;  428  // раскладка stack  429 1  mods = (word*)stack;  430 1  stack = mods + base_count;  431  // найти остатки  432 1  priBaseMod(mods, a, n, base_count);  433  // есть нулевые остатки?  434 1  while (base_count--)  435 1  if (mods[base_count] == 0)  436 1  return FALSE;  437  // нет  438 1  return TRUE;  439 }  440 441 1 size_t priIsSieved_deep(size_t base_count)  442 {  443 1  return O_OF_W(base_count);  444 }  445 446 1 bool_t priIsSmooth(const word a[], size_t n, size_t base_count, void* stack)  447 {  448  register size_t i;  449  register word mod;  450  // переменные в stack  451 1  word* t = (word*)stack;  452 1  stack = t + n;  453  // pre  454 1  ASSERT(base_count <= priBaseSize());  455  // t <- a  456 1  wwCopy(t, a, n);  457  // разделить t на степень 2  458 1  i = wwLoZeroBits(t, n);  459 1  wwShLo(t, n, i);  460 1  n = wwWordSize(t, n);  461 1  if (wwIsZero(t, n))  462  {  463 0  i = 0;  464 0  return TRUE;  465  }  466  // цикл по простым из факторной базы  467 1  for (i = 0; i < base_count;)  468  {  469 1  mod = _base[i] < WORD_BIT_HALF ?  470 1  zzModW2(t, n, _base[i]) : zzModW(t, n, _base[i]);  471  // делится на простое?  472 1  if (mod == 0)  473  {  474 1  zzDivW(t, t, n, _base[i]);  475 1  n = wwWordSize(t, n);  476 1  if (wwIsZero(t, n))  477  {  478 0  i = 0, mod = 0;  479 0  return TRUE;  480  }  481  }  482  // .. не делится  483  else  484 1  ++i;  485  }  486 1  i = 0, mod = 0;  487 1  return FALSE;  488 }  489 490 1 size_t priIsSmooth_deep(size_t n)  491 {  492 1  return O_OF_W(n);  493 }  494 495 /*  496 *******************************************************************************  497 Проверка простоты малых чисел  498 499 Применяется тест Миллера --- Рабина со специально подобранными основаниями.  500 Успешное завершение теста на всех основаниях из списка _base16 гарантирует  501 простоту чисел вплоть до 1373653, из списка _base32 --- вплоть  502 до 4759123141, из списка _base64 --- для всех 64-разрядных чисел.  503 504 Список оснований: http://miller-rabin.appspot.com (список обновляется).  505 506 \todo Проверка на малые делители позволит уменьшить число оснований /  507 сократить время проверки?  508 *******************************************************************************  509 */  510 511 const word _bases16[] = {2, 3};  512 const word _bases32[] = {2, 7, 61};  513 #if (B_PER_W == 64)  514 const word _bases64[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};  515 #endif  516 517 1 bool_t priIsPrimeW(register word a, void* stack)  518 {  519  const word* bases;  520  register word r;  521  register size_t s;  522  register size_t iter;  523  register size_t i;  524  register word base;  525  register dword prod;  526  // маленькое или четное?  527 1  if (a <= 3 || a % 2 == 0)  528 1  return a == 2 || a == 3;  529  // a - 1 = 2^s r (r -- нечетное)  530 1  for (r = a - 1, s = 0; r % 2 == 0; r >>= 1, ++s);  531 1  ASSERT(s > 0 && WORD_BIT_POS(s) * r + 1 == a);  532  // выбираем базис  533 #if (B_PER_W == 16)  534  bases = _bases16, iter = COUNT_OF(_bases16);  535 #elif (B_PER_W == 32)  536  if (a < 1373653)  537  bases = _bases16, iter = COUNT_OF(_bases16);  538  else  539  bases = _bases32, iter = COUNT_OF(_bases32);  540 #elif (B_PER_W == 64)  541 1  if (a < 1373653)  542 1  bases = _bases16, iter = COUNT_OF(_bases16);  543 0  else if (a < 4759123141)  544 0  bases = _bases32, iter = COUNT_OF(_bases32);  545  else  546 0  bases = _bases64, iter = COUNT_OF(_bases64);  547 #endif  548  // итерации  549 1  while (iter--)  550  {  551  // _bases[iter]^r \equiv \pm 1 \mod a?  552 1  base = zzPowerModW(bases[iter], r, a, stack);  553 1  if (base == 1 || base == a - 1)  554 1  continue;  555  // base^{2^i} \equiv - 1\mod a?  556 1  for (i = s - 1; i--;)  557  {  558 1  prod = base;  559 1  prod *= base, base = prod % a;  560 1  if (base == a - 1)  561 1  break;  562 1  if (base == 1)  563  {  564 0  r = base = 0, s = iter = i = 0, prod = 0;  565 0  return FALSE;  566  }  567  }  568 1  if (i == SIZE_MAX)  569  {  570 1  r = base = 0, s = iter = i = 0, prod = 0;  571 1  return FALSE;  572  }  573  }  574 1  r = base = 0, s = iter = i = 0, prod = 0;  575 1  return TRUE;  576 }  577 578 1 size_t priIsPrimeW_deep()  579 {  580 1  return zzPowerModW_deep();  581 }  582 583 /*  584 *******************************************************************************  585 Тест Рабина -- Миллера  586 587 При операциях \mod a может использоваться умножение Монтгомери. При этом  588 случайное основание выбирается отличным от R \mod a, -R \mod a (R = B^n),  589 а получаемые степени сравниваются не с \pm 1 \mod a, а с \pm R \mod a.  590 591 В начале функции проверяется, что a >= 49. Поэтому основание base будет  592 совпадать с \pm R \mod a с вероятностью p < 2/(49 - 1) = 1/24 и число попыток  593 генерации можно ограничить величиной  594 B_PER_IMPOSSLIBLE / log_2(24) < B_PER_IMPOSSLIBLE / 4.5.  595 *******************************************************************************  596 */  597 598 1 bool_t priRMTest(const word a[], size_t n, size_t iter, void* stack)  599 {  600  register size_t s;  601  register size_t m;  602  register size_t i;  603  // переменные в stack  604 1  word* r = (word*)stack;  605 1  word* base = r + n;  606 1  qr_o* qr = (qr_o*)(base + n);  607 1  octet* combo_state = (octet*)qr + zmCreate_keep(O_OF_W(n));  608 1  stack = combo_state + prngCOMBO_keep();  609  // pre  610 1  ASSERT(wwIsValid(a, n));  611  // нормализация  612 1  n = wwWordSize(a, n);  613  // четное?  614 1  if (zzIsEven(a, n))  615 0  return wwCmpW(a, n, 2) == 0;  616  // маленькое?  617 1  if (n == 1 && a[0] < 49)  618 1  return a[0] != 1 && (a[0] == 3 || a[0] % 3) && (a[0] == 5 || a[0] % 5);  619  // подготовить генератор  620 1  prngCOMBOStart(combo_state, utilNonce32());  621  // создать кольцо  622 1  wwTo(base, O_OF_W(n), a);  623 1  zmCreate(qr, (octet*)base, memNonZeroSize(base, O_OF_W(n)), stack);  624  // a - 1 = r 2^s (r -- нечетное)  625 1  wwCopy(r, a, n);  626 1  zzSubW2(r, n, 1);  627 1  s = wwLoZeroBits(r, n);  628 1  wwShLo(r, n, s);  629 1  m = wwWordSize(r, n);  630  // итерации  631 1  while (iter--)  632  {  633  // base <-R {1, \ldots, a - 1} \ {\pm one}  634 1  i = 0;  635  do  636 1  if (i++ * 45 > B_PER_IMPOSSIBLE * 10 ||  637 1  !zzRandNZMod(base, a, n, prngCOMBOStepR, combo_state))  638  {  639 0  s = m = 0;  640 0  return FALSE;  641  }  642 1  while (wwEq(base, qr->unity, n) || zzIsSumEq(a, base, qr->unity, n));  643  // base <- base^r \mod a  644 1  qrPower(base, base, r, m, qr, stack);  645  // base == \pm one => тест пройден  646 1  if (wwEq(base, qr->unity, n) ||  647 1  zzIsSumEq(a, base, qr->unity, n))  648 1  continue;  649  // base^{2^i} \equiv -1 \mod a?  650 1  for (i = s; i--;)  651  {  652 1  qrSqr(base, base, qr, stack);  653 1  if (wwEq(base, qr->unity, n))  654  {  655 0  s = m = i = 0;  656 0  return FALSE;  657  }  658 1  if (zzIsSumEq(a, base, qr->unity, n))  659 1  break;  660  }  661 1  if (i == SIZE_MAX)  662  {  663 1  s = m = i = 0;  664 1  return FALSE;  665  }  666  }  667 1  s = m = i = 0;  668  // простое  669 1  return TRUE;  670 }  671 672 1 size_t priRMTest_deep(size_t n)  673 {  674 1  size_t qr_deep = zmCreate_deep(O_OF_W(n));  675 1  return O_OF_W(2 * n) + zmCreate_keep(O_OF_W(n)) + prngCOMBO_keep() +  676 1  utilMax(2,  677  qr_deep,  678  qrPower_deep(n, n, qr_deep));  679 }  680 681 1 bool_t priIsPrime(const word a[], size_t n, void* stack)  682 {  683 1  return priRMTest(a, n, (B_PER_IMPOSSIBLE + 1) / 2, stack);  684 }  685 686 1 size_t priIsPrime_deep(size_t n)  687 {  688 1  return priRMTest_deep(n);  689 }  690 691 /*  692 *******************************************************************************  693 Простое Софи Жермен  694 695 q -- нечетное простое => p = 2q + 1 -- простое <=> [теорема Демитко]  696 1) 2^2 \not\equiv 1 \mod p;  697 2) 2^2q \equiv 1 \mod p.  698 *******************************************************************************  699 */  700 701 1 bool_t priIsSGPrime(const word q[], size_t n, void* stack)  702 {  703 1  size_t no = O_OF_W(n + 1);  704  // переменные в stack  705  word* p;  706  qr_o* qr;  707  // pre  708 1  ASSERT(zzIsOdd(q, n) && wwCmpW(q, n, 1) > 0);  709  // раскладка стек  710 1  p = (word*)stack;  711 1  qr = (qr_o*)(p + n + 1);  712 1  stack = (octet*)qr + zmCreate_keep(no);  713  // p <- 2q + 1  714 1  wwCopy(p, q, n);  715 1  p[n] = 0;  716 1  wwShHi(p, n + 1, 1);  717 1  ++p[0];  718  // создать кольцо \mod p  719 1  no = wwOctetSize(p, n + 1);  720 1  wwTo(p, no, p);  721 1  zmCreate(qr, (octet*)p, no, stack);  722  // p <- 4^q (в кольце qr)  723 1  qrAdd(p, qr->unity, qr->unity, qr);  724 1  qrAdd(p, p, p, qr);  725 1  qrPower(p, p, q, n, qr, stack);  726  // 4^q == 1 (в кольце qr)?  727 1  return qrCmp(p, qr->unity, qr) == 0;  728 }  729 730 1 size_t priIsSGPrime_deep(size_t n)  731 {  732 1  const size_t no = O_OF_W(n + 1);  733 1  const size_t qr_deep = zmCreate_deep(no);  734 1  return no + zmCreate_keep(no) +  735 1  utilMax(2,  736  qr_deep,  737  qrPower_deep(n + 1, n, qr_deep));  738 }  739 740 /*  741 *******************************************************************************  742 Следующее простое  743 *******************************************************************************  744 */  745 746 1 bool_t priNextPrimeW(word p[1], register word a, void* stack)  747 {  748  register size_t l;  749  // p <- a, l <- битовая длина a  750 1  p[0] = a, l = wwBitSize(p, 1);  751  // 0-битовых и 1-битовых простых не существует  752 1  if (l <= 1)  753 0  return FALSE;  754  // сделать p нечетным  755 1  p[0] |= 1;  756  // поиск  757 1  while (!priIsPrimeW(p[0], stack))  758  {  759 1  p[0] += 2;  760 1  if (wwBitSize(p, 1) != l)  761  {  762 0  l = 0;  763 0  return FALSE;  764  }  765  }  766 1  l = 0;  767 1  return TRUE;  768 }  769 770 1 size_t priNextPrimeW_deep()  771 {  772 1  return priIsPrimeW_deep();  773 }  774 775 1 bool_t priNextPrime(word p[], const word a[], size_t n, size_t trials,  776  size_t base_count, size_t iter, void* stack)  777 {  778  size_t l;  779  size_t i;  780  bool_t base_success;  781  // переменные в stack  782  word* mods;  783  // pre  784 1  ASSERT(wwIsSameOrDisjoint(a, p, n));  785 1  ASSERT(base_count <= priBaseSize());  786  // раскладка stack  787 1  mods = (word*)stack;  788 1  stack = mods + base_count;  789  // l <- битовая длина a  790 1  l = wwBitSize(a, n);  791  // 0-битовых и 1-битовых простых не существует  792 1  if (l <= 1)  793 0  return FALSE;  794  // p <- минимальное нечетное >= a  795 1  wwCopy(p, a, n);  796 1  p[0] |= 1;  797  // малое p?  798 1  if (n == 1)  799  // при необходимости скоррректировать факторную базу  800 1  while (base_count > 0 && priBasePrime(base_count - 1) >= p[0])  801 0  --base_count;  802  // рассчитать остатки от деления на малые простые  803 1  priBaseMod(mods, p, n, base_count);  804 1  for (i = 0, base_success = TRUE; i < base_count; ++i)  805 1  if (mods[i] == 0)  806  {  807 1  base_success = FALSE;  808 1  break;  809  }  810  // попытки  811 1  while (trials == SIZE_MAX || trials--)  812  {  813  // проверка простоты  814 1  if (base_success && priRMTest(p, n, iter, stack))  815 1  return TRUE;  816  // к следующему кандидату  817 1  if (zzAddW2(p, n, 2) || wwBitSize(p, n) > l)  818 1  return FALSE;  819 1  for (i = 0, base_success = TRUE; i < base_count; ++i)  820  {  821 1  if (mods[i] < _base[i] - 2)  822 1  mods[i] += 2;  823 1  else if (mods[i] == _base[i] - 1)  824 1  mods[i] = 1;  825  else  826 1  mods[i] = 0, base_success = FALSE;  827  }  828  }  829 0  return FALSE;  830 }  831 832 1 size_t priNextPrime_deep(size_t n, size_t base_count)  833 {  834 1  return base_count * O_PER_W + priRMTest_deep(n);  835 }  836 837 /*  838 *******************************************************************************  839 Расширение простого  840 841 Теорема Демитко. Если q -- нечетное простое, p = 2qr + 1, где 2r < 4q + 1,  842 и выполнены условия:  843 1) 2^{2qr} \equiv 1 \mod p;  844 2) 2^{2r} \not\equiv 1 \mod p,  845 то p -- простое.  846 847 \remark Если l <= 2 * lq, где lq = bitlen(q), то условие 2r < 4q + 1 будет  848 выполнено:  849 2r = (p - 1) / q < 2^l / (2^{lq - 1}} = 2^{l - lq + 1} <= 2^{lq + 2} < 4q.  850 851 Построение p:  852 1) t <-R {2^{l - 2} + 1,..., 2^{l - 1} - 1};  853 2) t <- t + 2^{l - 1};  854 3) r <- ceil(t / q);  855 4) p <- 2qr + 1;  856 5) если bitlen(p) != l, то вернуться к шагу 1.  857 858 \remark Если t укладывается в m слов, q -- в n слов, то r на шаге 3)  859 укладывается в m - n + 1 слов. Действительно, максимальное r получается  860 при t = B^m - 1, q = B^{n - 1} и равняется B^{m - n + 1} - 1.  861 *******************************************************************************  862 */  863 864 1 bool_t priExtendPrime(word p[], size_t l, const word q[], size_t n,  865  size_t trials, size_t base_count, gen_i rng, void* rng_state, void* stack)  866 {  867 1  const size_t m = W_OF_B(l);  868 1  const size_t mo = O_OF_B(l);  869  size_t i;  870  // переменные в stack  871  word* r;  872  word* t;  873  word* four;  874  word* mods;  875  word* mods1;  876  qr_o* qr;  877  // pre  878 1  ASSERT(wwIsDisjoint2(q, n, p, m));  879 1  ASSERT(zzIsOdd(q, n) && wwCmpW(q, n, 3) >= 0);  880 1  ASSERT(wwBitSize(q, n) + 1 <= l && l <= 2 * wwBitSize(q, n));  881 1  ASSERT(base_count <= priBaseSize());  882 1  ASSERT(rng != 0);  883  // подкорректировать n  884 1  n = wwWordSize(q, n);  885  // раскладка stack  886 1  r = (word*)stack;  887 1  t = r + m - n + 1;  888 1  four = t + m + 1;  889 1  mods = four + m;  890 1  mods1 = mods + base_count;  891 1  qr = (qr_o*)(mods1 + base_count);  892 1  stack = (octet*)qr + zmCreate_keep(mo);  893  // малое p?  894 1  if (l < B_PER_W)  895  // при необходимости уменьшить факторную базу  896 1  while (base_count > 0 &&  897 1  priBasePrime(base_count - 1) > WORD_BIT_POS(l - 1))  898 0  --base_count;  899  // попытки  900 1  while (trials == SIZE_MAX || trials--)  901  {  902  // t <-R [2^{l - 2}, 2^{l - 1})  903 1  rng(t, mo, rng_state);  904 1  wwFrom(t, t, mo);  905 1  wwTrimHi(t, m, l - 2);  906 1  wwSetBit(t, l - 2, 1);  907  // r <- t \div q  908 1  zzDiv(r, t, t, m, q, n, stack);  909 1  if (!wwIsZero(t, m))  910 1  VERIFY(zzAddW2(r, m - n + 1, 1) == 0);  911  // t <- q * r  912 1  zzMul(t, q, n, r, m - n + 1, stack);  913 1  if (wwBitSize(t, m) > l - 1)  914 1  continue;  915  // p <- 2 * t + 1  916 1  wwCopy(p, t, m);  917 1  wwShHi(p, m, 1);  918 1  ++p[0];  919 1  ASSERT(wwBitSize(p, m) == l);  920  // рассчитать вычеты p, 2q по малым модулям  921 1  priBaseMod(mods, p, m, base_count);  922 1  priBaseMod(mods1, q, n, base_count);  923 1  for (i = 0; i < base_count; ++i)  924 1  if ((mods1[i] += mods1[i]) >= _base[i])  925 1  mods1[i] -= _base[i];  926  // проверка простоты  927  while (1)  928  {  929  // p делится на малые простые?  930 1  for (i = 0; i < base_count; ++i)  931 1  if (mods[i] == 0)  932 1  break;  933  // не делится: тест Демитко  934 1  if (i == base_count)  935  {  936  // создать кольцо вычетов \mod p  937 1  wwTo(t, mo, p);  938 1  zmCreate(qr, (octet*)t, mo, stack);  939  // four <- 4 [в кольце qr]  940 1  qrAdd(four, qr->unity, qr->unity, qr);  941 1  qrAdd(four, four, four, qr);  942  // 4^r \mod p != 1?  943 1  qrPower(t, four, r, m - n + 1, qr, stack);  944 1  if (qrCmp(t, qr->unity, qr) != 0)  945  {  946  // (4^r)^q \mod p == 1?  947 1  qrPower(t, t, q, n, qr, stack);  948 1  if (qrCmp(t, qr->unity, qr) == 0)  949 1  return TRUE;  950  }  951  }  952  // p <- p + 2q, переполнение?  953 1  if (zzAddW2(p + n, m - n, zzAdd2(p, q, n)) ||  954 1  zzAddW2(p + n, m - n, zzAdd2(p, q, n)) ||  955 1  wwBitSize(p, m) > l)  956  break;  957  // r <- r + 1, t <- t + q, пересчитать mods  958 1  zzAddW2(r, m - n + 1, 1);  959 1  zzAddW2(t + n, m - n, zzAdd2(t, q, n));  960 1  for (i = 0; i < base_count; ++i)  961 1  if ((mods[i] += mods1[i]) >= _base[i])  962 1  mods[i] -= _base[i];  963  // к следующей попытке  964 1  if (trials != SIZE_MAX && trials-- == 0)  965 0  return FALSE;  966  }  967  }  968 0  return FALSE;  969 }  970 971 1 size_t priExtendPrime_deep(size_t l, size_t n, size_t base_count)  972 {  973 1  const size_t m = W_OF_B(l);  974 1  const size_t mo = O_OF_B(l);  975 1  const size_t qr_deep = zmCreate_deep(mo);  976 1  ASSERT(m >= n);  977 1  return O_OF_W(m - n + 1 + m + 1 + m + 2 * base_count) + zmCreate_keep(mo) +  978 1  utilMax(4,  979  zzDiv_deep(m, n),  980 1  zzMul_deep(n, m - n + 1),  981  qr_deep,  982  qrPower_deep(m, m, qr_deep));  983 } 

Read our documentation on viewing source code .