* Copyright (c) 1994 Landon Curt Noll
* Permission is granted to use, distribute, or modify this source,
* provided that this copyright notice remains intact.
* chongo@toad.com -or- ...!{pyramid,sun,uunet}!hoptoad!chongo
* primes of the form h*2^n-1 for 1<=h<200 and 1<=n<1000
* For all 0 <= i < prime_cnt, h_p[i]*2^n_p[i]-1 is prime.
* These values were taken from:
* "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
* Birkhauser, 1985, pp 384-387.
* This routine assumes that the file "lucas.cal" has been loaded.
* NOTE: There are several errors in Riesel's table that have been corrected
* 193*2^97-1 is NOT prime
* 199*2^221-1 is NOT prime
static prime_cnt = 1145; /* number of primes in the list */
/* h = prime parameters */
static mat h_p[prime_cnt] = {
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, /* element 0 */
1, 1, 1, 1, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 7, 7, 7, 7,
7, 7, 7, 7, 9, 9, 9, 9, 9, 9,
9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
9, 9, 9, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 13, 13, 13, 13, 13, 13, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, /* 100 */
15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 17, 17, 17, 17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
17, 17, 19, 19, 19, 19, 19, 19, 19, 19,
19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
19, 19, 21, 21, 21, 21, 21, 21, 21, 21,
21, 21, 21, 21, 21, 21, 21, 21, 23, 23,
23, 23, 23, 23, 23, 23, 23, 25, 25, 25,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
25, 25, 25, 27, 27, 27, 27, 27, 27, 27, /* 200 */
27, 27, 27, 27, 27, 27, 27, 27, 27, 27,
27, 27, 27, 27, 27, 27, 27, 29, 29, 29,
29, 29, 31, 31, 31, 31, 31, 31, 31, 31,
31, 31, 31, 31, 31, 31, 31, 31, 31, 31,
33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
33, 33, 33, 33, 35, 35, 35, 35, 35, 35,
35, 35, 35, 35, 35, 35, 35, 35, 35, 35,
35, 37, 39, 39, 39, 39, 39, 39, 39, 39,
39, 41, 41, 41, 41, 41, 41, 41, 41, 41, /* 300 */
41, 41, 41, 41, 43, 43, 43, 43, 43, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 47, 47, 47, 47, 49,
49, 49, 49, 49, 49, 49, 49, 49, 49, 49,
49, 49, 49, 49, 49, 49, 51, 51, 51, 51,
51, 51, 51, 51, 51, 51, 51, 51, 51, 51,
51, 53, 53, 53, 53, 53, 53, 53, 53, 53,
53, 55, 55, 55, 55, 55, 55, 55, 55, 55, /* 400 */
55, 55, 55, 55, 55, 55, 55, 55, 55, 55,
57, 57, 57, 57, 57, 57, 57, 57, 57, 57,
57, 57, 57, 57, 57, 57, 57, 57, 59, 59,
59, 59, 59, 59, 61, 61, 61, 61, 61, 61,
61, 61, 61, 61, 61, 61, 61, 61, 61, 61,
61, 63, 63, 63, 63, 63, 63, 63, 63, 63,
63, 63, 63, 63, 63, 63, 63, 63, 63, 63,
63, 63, 63, 63, 65, 65, 65, 65, 65, 65,
65, 65, 65, 65, 65, 65, 65, 65, 65, 65,
65, 65, 67, 67, 67, 67, 67, 67, 67, 67, /* 500 */
69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
69, 69, 71, 71, 71, 73, 73, 73, 73, 73,
73, 75, 75, 75, 75, 75, 75, 75, 75, 75,
75, 75, 75, 75, 75, 75, 75, 75, 75, 75,
75, 75, 75, 75, 75, 75, 75, 77, 77, 77,
77, 77, 77, 77, 77, 77, 77, 77, 77, 79,
79, 79, 79, 79, 79, 79, 79, 79, 79, 79,
81, 81, 81, 81, 81, 81, 81, 81, 81, 81, /* 600 */
81, 81, 81, 83, 83, 83, 83, 83, 83, 83,
83, 83, 83, 83, 83, 83, 83, 83, 83, 83,
83, 83, 83, 83, 83, 85, 85, 85, 85, 85,
85, 85, 85, 85, 87, 87, 87, 87, 87, 87,
87, 87, 87, 87, 87, 87, 87, 87, 87, 87,
87, 87, 87, 87, 87, 87, 89, 89, 89, 89,
89, 89, 89, 89, 89, 91, 91, 91, 91, 91,
91, 91, 91, 91, 91, 91, 91, 91, 91, 91,
91, 91, 91, 91, 91, 91, 91, 93, 93, 93,
93, 93, 93, 93, 93, 93, 93, 93, 93, 93, /* 700 */
93, 93, 93, 93, 93, 95, 95, 95, 95, 95,
95, 95, 95, 95, 95, 97, 97, 97, 97, 97,
99, 99, 99, 99, 99, 99, 99, 99, 99, 99,
99, 99, 99, 99, 99, 99, 101, 101, 101, 101,
103, 103, 103, 103, 103, 103, 103, 103, 103, 103,
103, 103, 103, 105, 105, 105, 105, 105, 105, 105,
105, 105, 105, 105, 105, 105, 105, 105, 105, 105,
105, 105, 107, 107, 107, 107, 107, 107, 107, 107,
107, 107, 107, 107, 107, 107, 109, 109, 109, 109,
109, 113, 113, 113, 113, 113, 113, 113, 113, 113, /* 800 */
113, 115, 115, 115, 115, 115, 115, 115, 115, 115,
115, 115, 115, 115, 115, 115, 115, 119, 119, 119,
119, 119, 119, 119, 119, 121, 121, 121, 121, 121,
121, 121, 121, 121, 121, 121, 121, 125, 125, 125,
125, 125, 125, 127, 127, 131, 131, 131, 131, 131,
131, 131, 131, 131, 131, 133, 133, 133, 133, 133,
133, 133, 133, 133, 133, 133, 133, 133, 137, 137,
137, 137, 139, 139, 139, 139, 139, 139, 139, 139,
139, 139, 139, 139, 139, 139, 139, 139, 139, 139,
139, 139, 139, 139, 139, 139, 139, 139, 139, 143, /* 900 */
143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
143, 143, 143, 145, 145, 145, 145, 145, 145, 145,
145, 145, 145, 145, 149, 149, 149, 149, 149, 149,
149, 151, 151, 151, 155, 155, 155, 155, 155, 155,
155, 155, 155, 155, 155, 155, 157, 157, 157, 157,
157, 157, 157, 157, 157, 161, 161, 161, 161, 161,
161, 161, 161, 161, 161, 161, 161, 161, 161, 161,
163, 163, 163, 163, 167, 167, 167, 167, 167, 167,
167, 167, 167, 167, 167, 167, 169, 169, 169, 169, /* 1000 */
169, 169, 169, 169, 169, 169, 169, 169, 173, 173,
173, 173, 173, 173, 173, 173, 173, 173, 173, 173,
173, 173, 173, 173, 175, 175, 175, 175, 175, 175,
175, 175, 175, 175, 175, 175, 175, 175, 175, 175,
179, 179, 179, 181, 181, 181, 181, 181, 181, 181,
181, 181, 181, 181, 181, 181, 181, 181, 181, 181,
181, 181, 181, 181, 181, 181, 181, 181, 185, 185,
185, 185, 185, 185, 185, 185, 185, 185, 187, 187,
187, 187, 187, 191, 193, 193, 193, 193, 193, 193,
193, 197, 197, 197, 197, 197, 197, 197, 197, 197, /* 1100 */
197, 197, 197, 197, 197, 197, 197, 197, 197, 199,
199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
/* n (exponent) prime parameters */
static mat n_p[prime_cnt] = {
2, 3, 5, 7, 13, 17, 19, 31, 61, 89, /* element 0 */
107, 127, 521, 607, 1, 2, 3, 4, 6, 7,
11, 18, 34, 38, 43, 55, 64, 76, 94, 103,
143, 206, 216, 306, 324, 391, 458, 470, 827, 2,
4, 8, 10, 12, 14, 18, 32, 48, 54, 72,
148, 184, 248, 270, 274, 420, 1, 5, 9, 17,
21, 29, 45, 177, 1, 3, 7, 13, 15, 21,
43, 63, 99, 109, 159, 211, 309, 343, 415, 469,
781, 871, 939, 2, 26, 50, 54, 126, 134, 246,
354, 362, 950, 3, 7, 23, 287, 291, 795, 1,
2, 4, 5, 10, 14, 17, 31, 41, 73, 80, /* 100 */
82, 116, 125, 145, 157, 172, 202, 224, 266, 289,
293, 463, 2, 4, 6, 16, 20, 36, 54, 60,
96, 124, 150, 252, 356, 460, 612, 654, 664, 698,
702, 972, 1, 3, 5, 21, 41, 49, 89, 133,
141, 165, 189, 293, 305, 395, 651, 665, 771, 801,
923, 953, 1, 2, 3, 7, 10, 13, 18, 27,
37, 51, 74, 157, 271, 458, 530, 891, 4, 6,
12, 46, 72, 244, 264, 544, 888, 3, 9, 11,
17, 23, 35, 39, 75, 105, 107, 155, 215, 335,
635, 651, 687, 1, 2, 4, 5, 8, 10, 14, /* 200 */
28, 37, 38, 70, 121, 122, 160, 170, 253, 329,
362, 454, 485, 500, 574, 892, 962, 4, 16, 76,
148, 184, 1, 5, 7, 11, 13, 23, 33, 35,
37, 47, 115, 205, 235, 271, 409, 739, 837, 887,
2, 3, 6, 8, 10, 22, 35, 42, 43, 46,
56, 91, 102, 106, 142, 190, 208, 266, 330, 360,
382, 462, 503, 815, 2, 6, 10, 20, 44, 114,
146, 156, 174, 260, 306, 380, 654, 686, 702, 814,
906, 1, 3, 24, 105, 153, 188, 605, 795, 813,
839, 2, 10, 14, 18, 50, 114, 122, 294, 362, /* 300 */
554, 582, 638, 758, 7, 31, 67, 251, 767, 1,
2, 3, 4, 5, 6, 8, 9, 14, 15, 16,
22, 28, 29, 36, 37, 54, 59, 85, 93, 117,
119, 161, 189, 193, 256, 308, 322, 327, 411, 466,
577, 591, 902, 928, 946, 4, 14, 70, 78, 1,
5, 7, 9, 13, 15, 29, 33, 39, 55, 81,
95, 205, 279, 581, 807, 813, 1, 9, 10, 19,
22, 57, 69, 97, 141, 169, 171, 195, 238, 735,
885, 2, 6, 8, 42, 50, 62, 362, 488, 642,
846, 1, 3, 5, 7, 15, 33, 41, 57, 69, /* 400 */
75, 77, 131, 133, 153, 247, 305, 351, 409, 471,
1, 2, 4, 5, 8, 10, 20, 22, 25, 26,
32, 44, 62, 77, 158, 317, 500, 713, 12, 16,
72, 160, 256, 916, 3, 5, 9, 13, 17, 19,
25, 39, 63, 67, 75, 119, 147, 225, 419, 715,
895, 2, 3, 8, 11, 14, 16, 28, 32, 39,
66, 68, 91, 98, 116, 126, 164, 191, 298, 323,
443, 714, 758, 759, 4, 6, 12, 22, 28, 52,
78, 94, 124, 162, 174, 192, 204, 304, 376, 808,
930, 972, 5, 9, 21, 45, 65, 77, 273, 677, /* 500 */
1, 4, 5, 7, 9, 11, 13, 17, 19, 23,
29, 37, 49, 61, 79, 99, 121, 133, 141, 164,
173, 181, 185, 193, 233, 299, 313, 351, 377, 540,
569, 909, 2, 14, 410, 7, 11, 19, 71, 79,
131, 1, 3, 5, 6, 18, 19, 20, 22, 28,
29, 39, 43, 49, 75, 85, 92, 111, 126, 136,
159, 162, 237, 349, 381, 767, 969, 2, 4, 14,
26, 58, 60, 64, 100, 122, 212, 566, 638, 1,
3, 7, 15, 43, 57, 61, 75, 145, 217, 247,
3, 5, 11, 17, 21, 27, 81, 101, 107, 327, /* 600 */
383, 387, 941, 2, 4, 8, 10, 14, 18, 22,
24, 26, 28, 36, 42, 58, 64, 78, 158, 198,
206, 424, 550, 676, 904, 5, 11, 71, 113, 115,
355, 473, 563, 883, 1, 2, 8, 9, 10, 12,
22, 29, 32, 50, 57, 69, 81, 122, 138, 200,
296, 514, 656, 682, 778, 881, 4, 8, 12, 24,
48, 52, 64, 84, 96, 1, 3, 9, 13, 15,
17, 19, 23, 47, 57, 67, 73, 77, 81, 83,
191, 301, 321, 435, 867, 869, 917, 3, 4, 7,
10, 15, 18, 19, 24, 27, 39, 60, 84, 111, /* 700 */
171, 192, 222, 639, 954, 2, 6, 26, 32, 66,
128, 170, 288, 320, 470, 1, 9, 45, 177, 585,
1, 4, 5, 7, 8, 11, 19, 25, 28, 35,
65, 79, 212, 271, 361, 461, 10, 18, 54, 70,
3, 7, 11, 19, 63, 75, 95, 127, 155, 163,
171, 283, 563, 2, 3, 5, 6, 8, 9, 25,
32, 65, 113, 119, 155, 177, 299, 335, 426, 462,
617, 896, 10, 12, 18, 24, 28, 40, 90, 132,
214, 238, 322, 532, 858, 940, 9, 149, 177, 419,
617, 8, 14, 74, 80, 274, 334, 590, 608, 614, /* 800 */
650, 1, 3, 11, 13, 19, 21, 31, 49, 59,
69, 73, 115, 129, 397, 623, 769, 12, 16, 52,
160, 192, 216, 376, 436, 1, 3, 21, 27, 37,
43, 91, 117, 141, 163, 373, 421, 2, 4, 44,
182, 496, 904, 25, 113, 2, 14, 34, 38, 42,
78, 90, 178, 778, 974, 3, 11, 15, 19, 31,
59, 75, 103, 163, 235, 375, 615, 767, 2, 18,
38, 62, 1, 5, 7, 9, 15, 19, 21, 35,
37, 39, 41, 49, 69, 111, 115, 141, 159, 181,
201, 217, 487, 567, 677, 765, 811, 841, 917, 2, /* 900 */
4, 6, 8, 12, 18, 26, 32, 34, 36, 42,
60, 78, 82, 84, 88, 154, 174, 208, 256, 366,
448, 478, 746, 5, 13, 15, 31, 77, 151, 181,
245, 445, 447, 883, 4, 16, 48, 60, 240, 256,
304, 5, 221, 641, 2, 8, 14, 16, 44, 46,
82, 172, 196, 254, 556, 806, 1, 5, 33, 121,
125, 305, 445, 473, 513, 2, 6, 18, 22, 34,
54, 98, 122, 146, 222, 306, 422, 654, 682, 862,
3, 31, 63, 303, 4, 6, 8, 10, 16, 32,
38, 42, 52, 456, 576, 668, 1, 5, 11, 17, /* 1000 */
67, 137, 157, 203, 209, 227, 263, 917, 2, 4,
6, 16, 32, 50, 76, 80, 96, 104, 162, 212,
230, 260, 480, 612, 1, 3, 9, 21, 23, 41,
47, 57, 69, 83, 193, 249, 291, 421, 433, 997,
8, 68, 108, 3, 5, 7, 9, 11, 17, 23,
31, 35, 43, 47, 83, 85, 99, 101, 195, 267,
281, 363, 391, 519, 623, 653, 673, 701, 2, 6,
10, 18, 26, 40, 46, 78, 230, 542, 1, 17,
21, 53, 253, 226, 3, 15, 27, 63, 87, 135,
543, 2, 16, 20, 22, 40, 82, 112, 178, 230, /* 1100 */
302, 304, 328, 374, 442, 472, 500, 580, 694, 1,
5, 7, 15, 19, 23, 25, 27, 43, 65, 99,
125, 141, 165, 201, 211, 331, 369, 389, 445, 461,
/* obtain our required libs */
* lucas_chk - check the lucas function on known primes
* This function tests entries in the above h_p, n_p table
* when n_p is below a given limit.
* high_n skip tests on n_p[i] > high_n
* [quiet] if given and != 0, then do not print individual test results
local result; /* 0 => non-prime, 1 => prime, -1 => bad test */
local error; /* number of errors and bad tests found */
ldebug("test_lucas", "high_n is non-int");
quit "FATAL: bad args: high_n must be an integer";
* scan thru the above prime table
for (i=0; i < prime_cnt; ++i) {
/* skip primes where h>=2^n */
if (highbit(h_p[i]) >= n_p[i]) {
print "h>=2^n skip:", h_p[i]:"*2^":n_p[i]:"-1";
/* test the prime if it is small enough */
/* test the table value */
result = lucas(h_p[i], n_p[i]);
print "ERROR, bad primality test of",\
h_p[i]:"*2^":n_p[i]:"-1";
} else if (result == 1) {
print h_p[i]:"*2^":n_p[i]:"-1 is prime";
} else if (result == -1) {
print "ERROR, failed to compute v(1) for",\
h_p[i]:"*2^":n_p[i]:"-1";
print "ERROR, bogus return value:", result;
/* return the full status */
print "lucas_chk(":high_n:") passed";
print "lucas_chk(":high_n:") failed", error, "test";
print "lucas_chk(":high_n:") failed", error, "tests";
print "lucas_chk(high_n) defined";
projects / unix-history / blob