1/* 2 * lucas_chk - test all primes of the form h*2^n-1, 1<=h<200 and n <= high_n 3 * 4 * Copyright (C) 1999,2021 Landon Curt Noll 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 1991/01/11 05:41:43 21 * File existed as early as: 1991 22 * 23 * chongo <was here> /\oo/\ http://www.isthe.com/chongo/ 24 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 25 */ 26 27/* 28 * primes of the form h*2^n-1 for 1<=h<200 and 1<=n<1000 29 * 30 * For all 0 <= i < prime_cnt, h_p[i]*2^n_p[i]-1 is prime. 31 * 32 * These values were taken from: 33 * 34 * "Prime numbers and Computer Methods for Factorization", by Hans Riesel, 35 * Birkhauser, 1985, pp 384-387. 36 * 37 * This routine assumes that the file "lucas.cal" has been loaded. 38 * 39 * NOTE: There are several errors in Riesel's table that have been corrected 40 * in this file: 41 * 42 * 193*2^87-1 is prime 43 * 193*2^97-1 is NOT prime 44 * 199*2^211-1 is prime 45 * 199*2^221-1 is NOT prime 46 */ 47 48 49static prime_cnt = 1145; /* number of primes in the list */ 50 51/* h = prime parameters */ 52static mat h_p[prime_cnt] = { 53 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, /* element 0 */ 54 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 55 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 56 3, 3, 3, 3, 3, 3, 3, 3, 3, 5, 57 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 58 5, 5, 5, 5, 5, 5, 7, 7, 7, 7, 59 7, 7, 7, 7, 9, 9, 9, 9, 9, 9, 60 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 61 9, 9, 9, 11, 11, 11, 11, 11, 11, 11, 62 11, 11, 11, 13, 13, 13, 13, 13, 13, 15, 63 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, /* 100 */ 64 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 65 15, 15, 17, 17, 17, 17, 17, 17, 17, 17, 66 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 67 17, 17, 19, 19, 19, 19, 19, 19, 19, 19, 68 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 69 19, 19, 21, 21, 21, 21, 21, 21, 21, 21, 70 21, 21, 21, 21, 21, 21, 21, 21, 23, 23, 71 23, 23, 23, 23, 23, 23, 23, 25, 25, 25, 72 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 73 25, 25, 25, 27, 27, 27, 27, 27, 27, 27, /* 200 */ 74 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 75 27, 27, 27, 27, 27, 27, 27, 29, 29, 29, 76 29, 29, 31, 31, 31, 31, 31, 31, 31, 31, 77 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 78 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 79 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 80 33, 33, 33, 33, 35, 35, 35, 35, 35, 35, 81 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 82 35, 37, 39, 39, 39, 39, 39, 39, 39, 39, 83 39, 41, 41, 41, 41, 41, 41, 41, 41, 41, /* 300 */ 84 41, 41, 41, 41, 43, 43, 43, 43, 43, 45, 85 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 86 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 87 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 88 45, 45, 45, 45, 45, 47, 47, 47, 47, 49, 89 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 90 49, 49, 49, 49, 49, 49, 51, 51, 51, 51, 91 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 92 51, 53, 53, 53, 53, 53, 53, 53, 53, 53, 93 53, 55, 55, 55, 55, 55, 55, 55, 55, 55, /* 400 */ 94 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 95 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 96 57, 57, 57, 57, 57, 57, 57, 57, 59, 59, 97 59, 59, 59, 59, 61, 61, 61, 61, 61, 61, 98 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 99 61, 63, 63, 63, 63, 63, 63, 63, 63, 63, 100 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 101 63, 63, 63, 63, 65, 65, 65, 65, 65, 65, 102 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 103 65, 65, 67, 67, 67, 67, 67, 67, 67, 67, /* 500 */ 104 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 105 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 106 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 107 69, 69, 71, 71, 71, 73, 73, 73, 73, 73, 108 73, 75, 75, 75, 75, 75, 75, 75, 75, 75, 109 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 110 75, 75, 75, 75, 75, 75, 75, 77, 77, 77, 111 77, 77, 77, 77, 77, 77, 77, 77, 77, 79, 112 79, 79, 79, 79, 79, 79, 79, 79, 79, 79, 113 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, /* 600 */ 114 81, 81, 81, 83, 83, 83, 83, 83, 83, 83, 115 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 116 83, 83, 83, 83, 83, 85, 85, 85, 85, 85, 117 85, 85, 85, 85, 87, 87, 87, 87, 87, 87, 118 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 119 87, 87, 87, 87, 87, 87, 89, 89, 89, 89, 120 89, 89, 89, 89, 89, 91, 91, 91, 91, 91, 121 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 122 91, 91, 91, 91, 91, 91, 91, 93, 93, 93, 123 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, /* 700 */ 124 93, 93, 93, 93, 93, 95, 95, 95, 95, 95, 125 95, 95, 95, 95, 95, 97, 97, 97, 97, 97, 126 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 127 99, 99, 99, 99, 99, 99, 101, 101, 101, 101, 128 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 129 103, 103, 103, 105, 105, 105, 105, 105, 105, 105, 130 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 131 105, 105, 107, 107, 107, 107, 107, 107, 107, 107, 132 107, 107, 107, 107, 107, 107, 109, 109, 109, 109, 133 109, 113, 113, 113, 113, 113, 113, 113, 113, 113, /* 800 */ 134 113, 115, 115, 115, 115, 115, 115, 115, 115, 115, 135 115, 115, 115, 115, 115, 115, 115, 119, 119, 119, 136 119, 119, 119, 119, 119, 121, 121, 121, 121, 121, 137 121, 121, 121, 121, 121, 121, 121, 125, 125, 125, 138 125, 125, 125, 127, 127, 131, 131, 131, 131, 131, 139 131, 131, 131, 131, 131, 133, 133, 133, 133, 133, 140 133, 133, 133, 133, 133, 133, 133, 133, 137, 137, 141 137, 137, 139, 139, 139, 139, 139, 139, 139, 139, 142 139, 139, 139, 139, 139, 139, 139, 139, 139, 139, 143 139, 139, 139, 139, 139, 139, 139, 139, 139, 143, /* 900 */ 144 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 145 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 146 143, 143, 143, 145, 145, 145, 145, 145, 145, 145, 147 145, 145, 145, 145, 149, 149, 149, 149, 149, 149, 148 149, 151, 151, 151, 155, 155, 155, 155, 155, 155, 149 155, 155, 155, 155, 155, 155, 157, 157, 157, 157, 150 157, 157, 157, 157, 157, 161, 161, 161, 161, 161, 151 161, 161, 161, 161, 161, 161, 161, 161, 161, 161, 152 163, 163, 163, 163, 167, 167, 167, 167, 167, 167, 153 167, 167, 167, 167, 167, 167, 169, 169, 169, 169, /* 1000 */ 154 169, 169, 169, 169, 169, 169, 169, 169, 173, 173, 155 173, 173, 173, 173, 173, 173, 173, 173, 173, 173, 156 173, 173, 173, 173, 175, 175, 175, 175, 175, 175, 157 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 158 179, 179, 179, 181, 181, 181, 181, 181, 181, 181, 159 181, 181, 181, 181, 181, 181, 181, 181, 181, 181, 160 181, 181, 181, 181, 181, 181, 181, 181, 185, 185, 161 185, 185, 185, 185, 185, 185, 185, 185, 187, 187, 162 187, 187, 187, 191, 193, 193, 193, 193, 193, 193, 163 193, 197, 197, 197, 197, 197, 197, 197, 197, 197, /* 1100 */ 164 197, 197, 197, 197, 197, 197, 197, 197, 197, 199, 165 199, 199, 199, 199, 199, 199, 199, 199, 199, 199, 166 199, 199, 199, 199, 199, 199, 199, 199, 199, 199, 167 199, 199, 199, 199, 199 168}; 169 170 171/* n (exponent) prime parameters */ 172static mat n_p[prime_cnt] = { 173 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, /* element 0 */ 174 107, 127, 521, 607, 1, 2, 3, 4, 6, 7, 175 11, 18, 34, 38, 43, 55, 64, 76, 94, 103, 176 143, 206, 216, 306, 324, 391, 458, 470, 827, 2, 177 4, 8, 10, 12, 14, 18, 32, 48, 54, 72, 178 148, 184, 248, 270, 274, 420, 1, 5, 9, 17, 179 21, 29, 45, 177, 1, 3, 7, 13, 15, 21, 180 43, 63, 99, 109, 159, 211, 309, 343, 415, 469, 181 781, 871, 939, 2, 26, 50, 54, 126, 134, 246, 182 354, 362, 950, 3, 7, 23, 287, 291, 795, 1, 183 2, 4, 5, 10, 14, 17, 31, 41, 73, 80, /* 100 */ 184 82, 116, 125, 145, 157, 172, 202, 224, 266, 289, 185 293, 463, 2, 4, 6, 16, 20, 36, 54, 60, 186 96, 124, 150, 252, 356, 460, 612, 654, 664, 698, 187 702, 972, 1, 3, 5, 21, 41, 49, 89, 133, 188 141, 165, 189, 293, 305, 395, 651, 665, 771, 801, 189 923, 953, 1, 2, 3, 7, 10, 13, 18, 27, 190 37, 51, 74, 157, 271, 458, 530, 891, 4, 6, 191 12, 46, 72, 244, 264, 544, 888, 3, 9, 11, 192 17, 23, 35, 39, 75, 105, 107, 155, 215, 335, 193 635, 651, 687, 1, 2, 4, 5, 8, 10, 14, /* 200 */ 194 28, 37, 38, 70, 121, 122, 160, 170, 253, 329, 195 362, 454, 485, 500, 574, 892, 962, 4, 16, 76, 196 148, 184, 1, 5, 7, 11, 13, 23, 33, 35, 197 37, 47, 115, 205, 235, 271, 409, 739, 837, 887, 198 2, 3, 6, 8, 10, 22, 35, 42, 43, 46, 199 56, 91, 102, 106, 142, 190, 208, 266, 330, 360, 200 382, 462, 503, 815, 2, 6, 10, 20, 44, 114, 201 146, 156, 174, 260, 306, 380, 654, 686, 702, 814, 202 906, 1, 3, 24, 105, 153, 188, 605, 795, 813, 203 839, 2, 10, 14, 18, 50, 114, 122, 294, 362, /* 300 */ 204 554, 582, 638, 758, 7, 31, 67, 251, 767, 1, 205 2, 3, 4, 5, 6, 8, 9, 14, 15, 16, 206 22, 28, 29, 36, 37, 54, 59, 85, 93, 117, 207 119, 161, 189, 193, 256, 308, 322, 327, 411, 466, 208 577, 591, 902, 928, 946, 4, 14, 70, 78, 1, 209 5, 7, 9, 13, 15, 29, 33, 39, 55, 81, 210 95, 205, 279, 581, 807, 813, 1, 9, 10, 19, 211 22, 57, 69, 97, 141, 169, 171, 195, 238, 735, 212 885, 2, 6, 8, 42, 50, 62, 362, 488, 642, 213 846, 1, 3, 5, 7, 15, 33, 41, 57, 69, /* 400 */ 214 75, 77, 131, 133, 153, 247, 305, 351, 409, 471, 215 1, 2, 4, 5, 8, 10, 20, 22, 25, 26, 216 32, 44, 62, 77, 158, 317, 500, 713, 12, 16, 217 72, 160, 256, 916, 3, 5, 9, 13, 17, 19, 218 25, 39, 63, 67, 75, 119, 147, 225, 419, 715, 219 895, 2, 3, 8, 11, 14, 16, 28, 32, 39, 220 66, 68, 91, 98, 116, 126, 164, 191, 298, 323, 221 443, 714, 758, 759, 4, 6, 12, 22, 28, 52, 222 78, 94, 124, 162, 174, 192, 204, 304, 376, 808, 223 930, 972, 5, 9, 21, 45, 65, 77, 273, 677, /* 500 */ 224 1, 4, 5, 7, 9, 11, 13, 17, 19, 23, 225 29, 37, 49, 61, 79, 99, 121, 133, 141, 164, 226 173, 181, 185, 193, 233, 299, 313, 351, 377, 540, 227 569, 909, 2, 14, 410, 7, 11, 19, 71, 79, 228 131, 1, 3, 5, 6, 18, 19, 20, 22, 28, 229 29, 39, 43, 49, 75, 85, 92, 111, 126, 136, 230 159, 162, 237, 349, 381, 767, 969, 2, 4, 14, 231 26, 58, 60, 64, 100, 122, 212, 566, 638, 1, 232 3, 7, 15, 43, 57, 61, 75, 145, 217, 247, 233 3, 5, 11, 17, 21, 27, 81, 101, 107, 327, /* 600 */ 234 383, 387, 941, 2, 4, 8, 10, 14, 18, 22, 235 24, 26, 28, 36, 42, 58, 64, 78, 158, 198, 236 206, 424, 550, 676, 904, 5, 11, 71, 113, 115, 237 355, 473, 563, 883, 1, 2, 8, 9, 10, 12, 238 22, 29, 32, 50, 57, 69, 81, 122, 138, 200, 239 296, 514, 656, 682, 778, 881, 4, 8, 12, 24, 240 48, 52, 64, 84, 96, 1, 3, 9, 13, 15, 241 17, 19, 23, 47, 57, 67, 73, 77, 81, 83, 242 191, 301, 321, 435, 867, 869, 917, 3, 4, 7, 243 10, 15, 18, 19, 24, 27, 39, 60, 84, 111, /* 700 */ 244 171, 192, 222, 639, 954, 2, 6, 26, 32, 66, 245 128, 170, 288, 320, 470, 1, 9, 45, 177, 585, 246 1, 4, 5, 7, 8, 11, 19, 25, 28, 35, 247 65, 79, 212, 271, 361, 461, 10, 18, 54, 70, 248 3, 7, 11, 19, 63, 75, 95, 127, 155, 163, 249 171, 283, 563, 2, 3, 5, 6, 8, 9, 25, 250 32, 65, 113, 119, 155, 177, 299, 335, 426, 462, 251 617, 896, 10, 12, 18, 24, 28, 40, 90, 132, 252 214, 238, 322, 532, 858, 940, 9, 149, 177, 419, 253 617, 8, 14, 74, 80, 274, 334, 590, 608, 614, /* 800 */ 254 650, 1, 3, 11, 13, 19, 21, 31, 49, 59, 255 69, 73, 115, 129, 397, 623, 769, 12, 16, 52, 256 160, 192, 216, 376, 436, 1, 3, 21, 27, 37, 257 43, 91, 117, 141, 163, 373, 421, 2, 4, 44, 258 182, 496, 904, 25, 113, 2, 14, 34, 38, 42, 259 78, 90, 178, 778, 974, 3, 11, 15, 19, 31, 260 59, 75, 103, 163, 235, 375, 615, 767, 2, 18, 261 38, 62, 1, 5, 7, 9, 15, 19, 21, 35, 262 37, 39, 41, 49, 69, 111, 115, 141, 159, 181, 263 201, 217, 487, 567, 677, 765, 811, 841, 917, 2, /* 900 */ 264 4, 6, 8, 12, 18, 26, 32, 34, 36, 42, 265 60, 78, 82, 84, 88, 154, 174, 208, 256, 366, 266 448, 478, 746, 5, 13, 15, 31, 77, 151, 181, 267 245, 445, 447, 883, 4, 16, 48, 60, 240, 256, 268 304, 5, 221, 641, 2, 8, 14, 16, 44, 46, 269 82, 172, 196, 254, 556, 806, 1, 5, 33, 121, 270 125, 305, 445, 473, 513, 2, 6, 18, 22, 34, 271 54, 98, 122, 146, 222, 306, 422, 654, 682, 862, 272 3, 31, 63, 303, 4, 6, 8, 10, 16, 32, 273 38, 42, 52, 456, 576, 668, 1, 5, 11, 17, /* 1000 */ 274 67, 137, 157, 203, 209, 227, 263, 917, 2, 4, 275 6, 16, 32, 50, 76, 80, 96, 104, 162, 212, 276 230, 260, 480, 612, 1, 3, 9, 21, 23, 41, 277 47, 57, 69, 83, 193, 249, 291, 421, 433, 997, 278 8, 68, 108, 3, 5, 7, 9, 11, 17, 23, 279 31, 35, 43, 47, 83, 85, 99, 101, 195, 267, 280 281, 363, 391, 519, 623, 653, 673, 701, 2, 6, 281 10, 18, 26, 40, 46, 78, 230, 542, 1, 17, 282 21, 53, 253, 226, 3, 15, 27, 63, 87, 135, 283 543, 2, 16, 20, 22, 40, 82, 112, 178, 230, /* 1100 */ 284 302, 304, 328, 374, 442, 472, 500, 580, 694, 1, 285 5, 7, 15, 19, 23, 25, 27, 43, 65, 99, 286 125, 141, 165, 201, 211, 331, 369, 389, 445, 461, 287 463, 467, 513, 583, 835 288}; 289 290 291/* obtain our required libs */ 292read -once "lucas.cal"; 293 294 295/* 296 * lucas_chk - check the lucas function on known primes 297 * 298 * This function tests entries in the above h_p, n_p table 299 * when n_p is below a given limit. 300 * 301 * input: 302 * high_n skip tests on n_p[i] > high_n 303 * [quiet] if given and != 0, then do not print individual test results 304 * 305 * returns: 306 * 1 all is OK 307 * 0 something went wrong 308 */ 309define 310lucas_chk(high_n, quiet) 311{ 312 local i; /* index */ 313 local result; /* 0 => non-prime, 1 => prime, -1 => bad test */ 314 local error; /* number of errors and bad tests found */ 315 316 /* 317 * firewall 318 */ 319 if (!isint(high_n)) { 320 ldebug("test_lucas", "high_n is non-int"); 321 quit "FATAL: bad args: high_n must be an integer"; 322 } 323 if (param(0) == 1) { 324 quiet = 0; 325 } 326 327 /* 328 * scan thru the above prime table 329 */ 330 error = 0; 331 for (i=0; i < prime_cnt; ++i) { 332 333 /* skip primes where h>=2^n */ 334 if (highbit(h_p[i]) >= n_p[i]) { 335 if (config("resource_debug") & 8) { 336 print "h>=2^n skip:", h_p[i]:"*2^":n_p[i]:"-1"; 337 } 338 continue; 339 } 340 341 /* test the prime if it is small enough */ 342 if (n_p[i] <= high_n) { 343 344 /* test the table value */ 345 result = lucas(h_p[i], n_p[i]); 346 347 /* report the test */ 348 if (result == 0) { 349 print "ERROR, bad primality test of",\ 350 h_p[i]:"*2^":n_p[i]:"-1"; 351 ++error; 352 } else if (result == 1) { 353 if (quiet == 0) { 354 print h_p[i]:"*2^":n_p[i]:"-1 is prime"; 355 } 356 } else if (result == -1) { 357 print "ERROR, failed to compute v(1) for",\ 358 h_p[i]:"*2^":n_p[i]:"-1"; 359 ++error; 360 } else { 361 print "ERROR, bogus return value:", result; 362 ++error; 363 } 364 } 365 } 366 367 /* return the full status */ 368 if (error == 0) { 369 if (quiet == 0) { 370 print "lucas_chk(":high_n:") passed"; 371 } 372 return 1; 373 } else if (error == 1) { 374 print "lucas_chk(":high_n:") failed", error, "test"; 375 return 0; 376 } else { 377 print "lucas_chk(":high_n:") failed", error, "tests"; 378 return 0; 379 } 380} 381