/
tst-prime.lisp
543 lines (496 loc) · 32.4 KB
/
tst-prime.lisp
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
;; -*- Mode:Lisp; Syntax:ANSI-Common-LISP; Coding:us-ascii-unix; fill-column:158 -*-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;; @file tst-prime.lisp
;; @author Mitch Richling <https://www.mitchr.me>
;; @brief Unit tests.@EOL
;; @std Common Lisp
;; @see use-prime.lisp
;; @copyright
;; @parblock
;; Copyright (c) 1992,1994,1997,1998,2004,2008,2013,2015, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
;;
;; Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
;;
;; 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
;;
;; 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
;; and/or other materials provided with the distribution.
;;
;; 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
;; without specific prior written permission.
;;
;; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
;; LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
;; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
;; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
;; DAMAGE.
;; @endparblock
;; @todo mjr_prime_pimep-miller-rabin-deterministic: Correct refs and fix name.@EOL@EOL
;; @todo mjr_prime_pimep-miller-rabin-deterministic: Correct refs.@EOL@EOL
;; @todo mjr_prime_probable-primep-lucas-selfridge: Finish documentation string.@EOL@EOL
;; @todo mjr_prime_probable-primep-lucas-selfridge: Add alternative selection for D.@EOL@EOL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defpackage :MJR_PRIME-TESTS (:USE :COMMON-LISP :LISP-UNIT :MJR_PRIME :MJR_PRNG :MJR_INTU))
(in-package :MJR_PRIME-TESTS)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Initialize the small prime list..
(if (null mjr_prime::*mjr_prime_small-list*)
(mjr_prime_init-small-prime-list))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prime_all-factors-naive (n)
"Return a list of all of the factors of the given number.
Prime numbers will have a return list that contains only two elements (the input number and the number 1).
This function is very slow, but provides a regression test for mjr_prime_all-factors -- use that function instead of this one."
(if (integerp n)
(let ((small-factors (loop for i from 1 upto (isqrt (abs n)) when (= 0 (mod n i)) collect i)))
(remove-duplicates (append small-factors (reverse (mapcar (lambda (x) (/ (abs n) x)) small-factors)))))
(error "mjr_prime_all-factors: input must be an integer")))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prime_jacobi-symbol-naive (a n)
(reduce #'* (mapcar (lambda (x) (expt (mjr_prime_legendre-symbol a (car x)) (cdr x))) (mjr_prime_prime-factorization n))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prime_lucas-sequence-naive (n p q modn retu)
""
(flet ((modit (k) (if modn (mod k modn) k)))
(let ((i 1) ;; Start with i=1
(ui (modit 1)) ;; $U_{1}=1$
(ui+1 (modit P)) ;; $U_{2}=P$
(vi (modit P)) ;; $V_{1}=P$
(vi+1 (modit (- (* P P) (* 2 Q))))) ;; $V_{2}=P^2-2Q$
(cond ((= n 0) (if retu 0 2)) ;; $U_{0}=0$, $V_{0}=2$
((= n 1) (if retu ui vi))
((= n 2) (if retu ui+1 vi+1))
((> n 1) (progn
(loop for qti = (if modn (mjr_intu_mod-expt q i modn) (expt q i))
while (< (* i 2) (1+ n))
do (psetq i (* i 2)
ui (modit (* ui vi)) ;; $$U_{2i} = U_i V_i$$
ui+1 (modit (- (* ui+1 vi) qti)) ;; $$U_{2i+1} = U_{i+1} V_i - Q^i$$
vi (modit (- (* vi vi) (* 2 qti))) ;; $$V_{2i} = V_i^2 - 2Q^i$$
vi+1 (modit (- (* vi+1 vi) (* p qti))))) ;; $$V_{2i+1} = V_{i+1} V_i - PQ^i$$
(loop while (< i n)
do (psetq i (+ i 1)
ui (if retu ui+1)
ui+1 (if retu (modit (- (* p ui+1) (* q ui)))) ;; $U_n = P U_{n-1} - Q U_{n-2}$
vi (if (not retu) vi+1)
vi+1 (if (not retu) (modit (- (* p vi+1) (* q vi)))))) ;; $V_n = P V_{n-1} - Q V_{n-2}$
(if retu
ui
vi)))))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defvar *mjr_prime_small-list-naive* nil
"An array of the first few prime numbers.")
(defvar *mjr_prime_small-bitmap-naive* nil
"An bitmap of the first few prime numbers.")
(defvar *mjr_prime_small-count-naive* 0
"Number of prime numbers in the *mjr_prime_small-list-naive* array. Normally 6 million, but may be smaller")
(defvar *mjr_prime_small-max-naive* 0
"Largest prime number in the *mjr_prime_small-list-naive* array")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prime_init-small-prime-list-naive (&optional (max-prime-can 104395301))
"This function uses a slow and RAM hungry implementation of the sieve of Atkin to do what MJR_PRIME_INIT-SMALL-PRIME-LIST can do.
This function provides a simplistic implementation of a sieve using an entirely different algorithm so that we can use this
function as a part of a regression test to help verify the correctness of MJR_PRIME_INIT-SMALL-PRIME-LIST.
References:
Atkin & Bernstein (2004); Prime sieves using binary quadratic forms; Mathematics of Computation 73."
(if (> max-prime-can array-dimension-limit) (error "mjr_prime_init-small-prime-list-naive: The array size limit is too small for such a large sieve!"))
(let ((abv (make-array (1+ max-prime-can) :element-type 'bit :initial-element 0)))
(mapc (lambda (n) (setf (aref abv n) 1)) (list 2 3))
(loop for x from 1 upto (isqrt max-prime-can) ; Solve the equation
do (loop for y from 1 upto (isqrt max-prime-can)
do (let* ((n (+ (* 4 x x) (* y y))))
(if (and (< n max-prime-can) (member (mod n 12) (list 1 5)))
(setf (aref abv n) (if (= 1 (aref abv n)) 0 1))))
do (let* ((n (+ (* 3 x x) (* y y))))
(if (and (< n max-prime-can) (= (mod n 12) 7))
(setf (aref abv n) (if (= 1 (aref abv n)) 0 1))))
do (let* ((n (- (* 3 x x) (* y y))))
(if (and (> x y) (< n max-prime-can) (= (mod n 12) 11))
(setf (aref abv n) (if (= 1 (aref abv n)) 0 1))))))
(loop for i from 5 upto max-prime-can ; Keep only square-free numbers
when (= 1 (aref abv i))
do (loop for j from (* i i) upto max-prime-can by (* i i)
do (setf (aref abv j) 0)))
(setq *mjr_prime_small-list-naive* (make-array (loop for b bit across abv when (= b 1) count 1)
:element-type 'fixnum
:initial-contents (loop for i from 2 upto max-prime-can when (= 1 (aref abv i)) collect i)))
(setq *mjr_prime_small-count-naive* (length *mjr_prime_small-list-naive*))
(setq *mjr_prime_small-max-naive* (aref *mjr_prime_small-list-naive* (1- *mjr_prime_small-count-naive*)))
(setq *mjr_prime_small-bitmap-naive* abv)
nil))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_init-small-prime-list
; Make sure small prime list gets created.
(assert-true mjr_prime::*mjr_prime_small-list*)
(assert-true mjr_prime::*mjr_prime_small-bitmap*)
(assert-true (< 0 mjr_prime::*mjr_prime_small-count*))
(assert-true (< 0 mjr_prime::*mjr_prime_small-max*))
;; atkin will not get a list of the same length given the same starting point, so we check the values that are computed by both
(mjr_prime_init-small-prime-list-naive)
(assert-equalp (length mjr_prime::*mjr_prime_small-list*) mjr_prime::*mjr_prime_small-count*)
(assert-equalp (length *mjr_prime_small-list-naive*) *mjr_prime_small-count-naive*)
(assert-equalp 5999999 (length *mjr_prime_small-list-naive*))
(assert-equalp 6000001 (length mjr_prime::*mjr_prime_small-list*))
(assert-equalp 27838752 (length mjr_prime::*mjr_prime_small-bitmap*))
(dotimes (i (min (length mjr_prime::*mjr_prime_small-list*) (length *mjr_prime_small-list-naive*)))
(assert-equalp (aref *mjr_prime_small-list-naive* i) (aref mjr_prime::*mjr_prime_small-list* i) i))
(setq *mjr_prime_small-list-naive* nil)
(setq *mjr_prime_small-bitmap-naive* nil)
(setq *mjr_prime_small-count-naive* 0)
(setq *mjr_prime_small-max-naive* 0)
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_nth-small-prime
; Too small
(assert-equal nil (mjr_prime_nth-small-prime -1))
(assert-equal 2 (mjr_prime_nth-small-prime 0))
; Good indexes
;; OEIS A000040
(loop with lst = '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149
151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271)
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_nth-small-prime i) i))
; Too big
(assert-equal nil (mjr_prime_nth-small-prime mjr_prime::*mjr_prime_small-count*))
(assert-equal nil (mjr_prime_nth-small-prime (1+ mjr_prime::*mjr_prime_small-count*)))
(assert-equal nil (mjr_prime_nth-small-prime (expt 2 64)))
;; Make sure mjr_prime_small-prime-index and mjr_prime_nth-small-prime are consistant
(dotimes (i 200)
(let ((r (random mjr_prime::*mjr_prime_small-count*)))
(assert-equal r (mjr_prime_small-prime-index (mjr_prime_nth-small-prime r)))))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_small-prime-index
;; XREF: Tested in mjr_prime_nth-small-prime
;; TODO: Code converage should be increased
1
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_all-factors
(dotimes (i 50)
(let ((j (random 1000000000)))
(assert-equalp (mjr_prime_all-factors j) (mjr_prime_all-factors-naive j))))
(assert-equal '(1) (mjr_prime_all-factors 1))
(assert-equal '(1 2) (mjr_prime_all-factors 2))
(assert-equal '(1 3) (mjr_prime_all-factors 3))
(assert-equal '(1 2 4) (mjr_prime_all-factors 4))
(assert-equal '(1 5) (mjr_prime_all-factors 5))
(assert-equal '(1 2 3 6) (mjr_prime_all-factors 6))
(assert-equal '(1 7) (mjr_prime_all-factors 7))
(assert-equal '(1 2 4 8) (mjr_prime_all-factors 8))
(assert-equal '(1 3 5 15) (mjr_prime_all-factors 15))
(assert-equal '(1 2 3 5 6 10 15 30) (mjr_prime_all-factors 30))
(assert-equal '(1 2 3 4 6 9 12 18 36) (mjr_prime_all-factors 36))
;; Errors
(assert-error 'error (mjr_prime_all-factors 't))
;; MJR TODO NOTE mjr_prime_all-factors: Fix next test
;(assert-error 'error (mjr_prime_all-factors 10.0))
(assert-error 'error (mjr_prime_all-factors #C(1 1)))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_deficient?
;; OEIS A005100
(loop with lst = '(1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 21 22 23 25 26 27 29 31 32 33 34 35 37 38 39 41 43 44 45 46 47 49 50 51
52 53 55 57 58 59 61 62 63 64 65 67 68 69 71 73 74 75 76 77 79 81 82 83 85 86)
with c = (car lst)
for n from 1
while c
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-true (mjr_prime_deficient? n))
(assert-true (mjr_prime_deficient? n)))))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_abundant?
;; OEIS A005101
(loop with lst = '(12 18 20 24 30 36 40 42 48 54 56 60 66 70 72 78 80 84 88 90 96 100 102 104 108 112 114 120 126 132 138 140 144
150 156 160 162 168 174 176 180 186 192 196 198 200 204 208 210 216 220 222 224 228 234 240 246 252 258 260
264 270)
with c = (car lst)
for n from 1
while c
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-true (mjr_prime_abundant? n))
(assert-true (mjr_prime_abundant? n)))))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_perfect?
;; OEIS A000396
(loop with lst = '(6 28 496 8128 33550336)
with c = (car lst)
for n from 1
while c
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-true (mjr_prime_perfect? n))
(assert-true (mjr_prime_perfect? n)))))
;; OEIS A000396
(loop for n in '(8589869056 137438691328 2305843008139952128)
do (assert-true (mjr_prime_perfect? n)))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_strong-probable-primep
(loop for (a . lst) in '((2 . (2047 3277 4033 4681 8321 15841 29341 42799 49141 52633 65281 74665 80581 85489 88357 90751
104653 130561 196093 220729 233017 252601 253241 256999 271951 280601 314821 357761 390937
458989 476971 486737))
(3 . (121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345 23521 31621 44287 47197 55969
63139 74593 79003 82513 87913 88573 97567 105163 111361 112141 148417 152551 182527 188191
211411 218791 221761 226801))
(5 . (781 1541 5461 5611 7813 13021 14981 15751 24211 25351 29539 38081 40501 44801 53971 79381
100651 102311 104721 112141 121463 133141 141361 146611 195313 211951 216457 222301 251521
289081 290629 298271 315121))
(7 . (25 325 703 2101 2353 4525 11041 14089 20197 29857 29891 39331 49241 58825 64681 76627
78937 79381 87673 88399 88831 102943 109061 137257 144901 149171 173951 178709 188191
197633 219781 227767 231793 245281))
(11 . (133 793 2047 4577 5041 12403 13333 14521 17711 23377 43213 43739 47611 48283 49601 50737
50997 56057 58969 68137 74089 85879 86347 87913 88831 102173 111055 114211 115231 137149
139231 171601 172369 193249 196555))
(13 . (85 1099 5149 7107 8911 9637 13019 14491 17803 19757 20881 22177 23521 26521 35371 44173
45629 54097 56033 57205 75241 83333 85285 86347 102719 110309 153401 184339 191959 222529
242845 253021 253927 269861))
(17 . (9 91 145 781 1111 2821 4033 4187 5365 5833 6697 7171 15805 19729 21781 22791 24211 26245
31621 33001 33227 34441 35371 38081 42127 49771 71071 74665 77293 78881 88831 96433 97921
98671 101101 102311 125563 129493))
(19 . (9 49 169 343 1849 2353 2701 4033 4681 6541 6697 7957 9997 12403 13213 13747 15251 16531
18769 19729 24761 30589 31621 31861 32477 41003 49771 63139 64681 65161 66421 68257 73555
96049 102831 118957 129961 137311))
;;((2 3 5 7) . (3215031751 118670087467 307768373641 315962312077 354864744877 457453568161
;; 528929554561 546348519181 602248359169))
;;((2 3 5) . (25326001 161304001 960946321 1157839381 3215031751 3697278427 5764643587 6770862367
;; 14386156093 15579919981 18459366157 19887974881 21276028621))
;;((2 3) . (1373653 1530787 1987021 2284453 3116107 5173601 6787327 11541307 13694761 15978007
;; 16070429 16879501 25326001 27509653 27664033 28527049 54029741 61832377 66096253))
((2 3) . (1373653 1530787 1987021))
)
do (loop with mx = (apply #'max lst)
with c = (car lst)
for n from 1 upto mx
for b = (and (not (mjr_prime_primep n)) (mjr_prime_strong-probable-primep n a))
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-true b))
(assert-false b))))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_pimep-miller-rabin-deterministic
(loop for i from 1 upto (min mjr_prime::*mjr_prime_small-max* 500000)
do (assert-equalp (mjr_prime_primep-small i) (mjr_prime_pimep-miller-rabin-deterministic i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_legendre-symbol
(assert-equalp 0 (mjr_prime_legendre-symbol 5 5))
(assert-equalp 0 (mjr_prime_legendre-symbol 0 5))
(assert-equalp -1 (mjr_prime_legendre-symbol 2 3))
(assert-equalp -1 (mjr_prime_legendre-symbol 2 5))
(assert-equalp 0 (mjr_prime_legendre-symbol 10 5))
(assert-equalp 1 (mjr_prime_legendre-symbol 10 3))
(assert-equalp -1 (mjr_prime_legendre-symbol -3 5))
(assert-equalp 1 (mjr_prime_legendre-symbol -3 2))
(assert-equalp -1 (mjr_prime_legendre-symbol 5 3))
(assert-equalp -1 (mjr_prime_legendre-symbol 5 7))
(assert-equalp -1 (mjr_prime_legendre-symbol 2 5))
(assert-equalp -1 (mjr_prime_legendre-symbol 3 5))
(assert-equalp 0 (mjr_prime_legendre-symbol 5 5))
(assert-equalp -1 (mjr_prime_legendre-symbol 7 5))
(assert-equalp 1 (mjr_prime_legendre-symbol 11 5))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_jacobi-symbol
(assert-equalp 1 (mjr_prime_jacobi-symbol 1236 20003))
(assert-equalp 1 (mjr_prime_jacobi-symbol 1236 1111))
(assert-equalp 0 (mjr_prime_jacobi-symbol 1236 21))
(assert-equalp -1 (mjr_prime_jacobi-symbol 1236 221))
(assert-equalp -1 (mjr_prime_jacobi-symbol 1001 9907))
;; mjr_prime_jacobi-symbol === mjr_prime_legendre-symbol when second argument is an odd prime...
(loop for i from 1 upto 1000
for p = (aref mjr_prime::*mjr_prime_small-list* i) ;; primes from 3 on
for a = (mjr_prng_int-cc 0 10000)
do (assert-equalp (mjr_prime_jacobi-symbol a p) (mjr_prime_legendre-symbol a p) a p))
;; Naive version should match production version
(loop for i from 1 upto 1000
for noe = (mjr_prng_int-cc 0 10000)
for n = (if (evenp noe) (1+ noe) noe)
for a = (mjr_prng_int-cc 0 10000)
do (assert-equalp (mjr_prime_jacobi-symbol-naive a n) (mjr_prime_jacobi-symbol a n) a n))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_probable-primep-lucas-selfridge
;; First 10 are from "Perfect, Amicable, and Sociable Numbers: A Computational Approach" By Song Y. Yan
(loop with lst = '(323 377 1159 1829 3827 5459 5777 9071 9179 10877 11419 11663 13919 14839 16109 16211 18407 18971 19043 22499
23407 24569 25199 25877 26069 27323 32759 34943 35207 39059 39203 39689 40309 44099 46979 47879)
with c = (car lst)
for n from 1
for bs = (mjr_prime_primep-small n)
for bl = (mjr_prime_probable-primep-lucas-selfridge n)
while c
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-false (equalp bs bl) n))
(assert-true (equalp bs bl) n)))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_probable-primep-bpsw
(assert-false (mjr_prime_probable-primep-bpsw 3825123056546413051))
(loop for i from 1 upto (min mjr_prime::*mjr_prime_small-max* 500000)
do (assert-equalp (mjr_prime_primep-small i) (mjr_prime_probable-primep-bpsw i) i))
;; A000668 Mersenne primes
(loop for i in '(3 7 31 127 8191 131071 524287 2147483647 2305843009213693951 618970019642690137449562111
162259276829213363391578010288127 170141183460469231731687303715884105727)
do (assert-true (mjr_prime_probable-primep-bpsw i)))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_lucas-sequence
;; Random tests
(dotimes (i 10000)
(let ((n (mjr_prng_int-cc 0 10000))
(p (mjr_prng_int-cc -10 10))
(q (mjr_prng_int-cc -10 10))
(m (mjr_prng_int-cc 100 1000)))
(assert-equalp (mjr_prime_lucas-sequence-naive n p q m 't) (mjr_prime_lucas-sequence n p q m 't) n p q m)))
;; A few lucas sequences with closed form solutions for the terms
(dotimes (k 1000)
(assert-equalp (mjr_prime_lucas-sequence k 3 2 nil 't) (1- (expt 2 k)) k) ;; p=3, q=2, U_k=2^k-1 (OEIS A000225)
(assert-equalp (mjr_prime_lucas-sequence k 3 2 nil nil) (1+ (expt 2 k)) k) ;; p=3, q=2, V_k=2^k+1 (OEIS A000051)
(assert-equalp (mjr_prime_lucas-sequence k 5 4 nil 't) (/ (1- (expt 4 k)) 3) k) ;; p=5, q=4, U_k=(4^k - 1)/3 (OEIS A002450)
(assert-equalp (mjr_prime_lucas-sequence k 5 4 nil nil) (1+ (expt 4 k)) k) ;; p=5, q=4, V_k=(4^k + 1) (OEIS A052539)
(assert-equalp (mjr_prime_lucas-sequence k 4 3 nil 't) (/ (1- (expt 3 k)) 2) k) ;; p=4, q=3, U_k=(3^n - 1)/2 (OEIS A003462)
(assert-equalp (mjr_prime_lucas-sequence k 4 3 nil nil) (1+ (expt 3 k)) k)) ;; p=4, q=3, V_k=(3^k + 1) (OEIS A034472)
;; NOTE: This function is tested heavily by other tests too: mjr_prime_fibonacci-number mjr_prime_lucas-number
;; mjr_prime_pell-lucas-number mjr_prime_pell-number mjr_prime_jacobsthal-number mjr_prime_jacobsthal-lucas-number
;; mjr_prime_probable-primep-lucas-selfridge mjr_prime_probable-primep-bpsw
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_fibonacci-number
;; Sloan A000045
(loop with lst = '(0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393
196418 317811 514229 832040 1346269 2178309 3524578 5702887 9227465 14930352 24157817 39088169)
for i from 0 upto (1- (length lst))
for f in lst
do (assert-equalp f (mjr_prime_fibonacci-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_lucas-number
;; Sloan A000204
(loop with lst = '(1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364 2207 3571 5778 9349 15127 24476 39603 64079 103682 167761
271443 439204 710647 1149851 1860498 3010349 4870847 7881196 12752043 20633239 33385282 54018521 87403803
141422324)
for i from 1 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_lucas-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_pell-lucas-number
;; Sloan A002203
(loop with lst = '(2 2 6 14 34 82 198 478 1154 2786 6726 16238 39202 94642 228486 551614 1331714 3215042 7761798 18738638 45239074
109216786 263672646 636562078 1536796802 3710155682 8957108166 21624372014 52205852194 126036076402
304278004998 )
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_pell-lucas-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_pell-number
;; Sloan A000129
(loop with lst = '(0 1 2 5 12 29 70 169 408 985 2378 5741 13860 33461 80782 195025 470832 1136689 2744210 6625109 15994428
38613965 93222358 225058681 543339720 1311738121 3166815962 7645370045 18457556052 44560482149 107578520350
259717522849)
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_pell-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_jacobsthal-number
;; Sloan A001045
(loop with lst = '(0 1 1 3 5 11 21 43 85 171 341 683 1365 2731 5461 10923 21845 43691 87381 174763 349525 699051 1398101 2796203
5592405 11184811 22369621 44739243 89478485 178956971 357913941 715827883 1431655765 2863311531)
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_jacobsthal-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_jacobsthal-lucas-number
;; Sloan A014551
(loop with lst = '(2 1 5 7 17 31 65 127 257 511 1025 2047 4097 8191 16385 32767 65537 131071 262145 524287 1048577 2097151 4194305
8388607 16777217 33554431 67108865 134217727 268435457 536870911 1073741825 2147483647 4294967297 8589934591)
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_jacobsthal-lucas-number i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_phi-func
;; OEIS A000010
(loop with lst = '(1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 12 18 12 28 8 30 16 20 16 24 12 36 18 24 16 40 12 42
20 24 22 46 16 42 20 32 24 52 18 40 24 36 28 58 16 60 30 36 32 48 20 66 32 44 )
for i from 1 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_phi-func i) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_pi-func
;; OEIS A006880 Number of primes < 10^n.
(loop with lst = '( 0 4 25 168 1229 9592 78498 664579 5761455 50847534 455052511 4118054813 37607912018 346065536839 3204941750802
29844570422669 279238341033925 2623557157654233 24739954287740860 234057667276344607 2220819602560918840)
for i from 0 upto (length lst)
for f in lst
do (assert-equalp f (mjr_prime_pi-func (expt 10 i)) i))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(define-test mjr_prime_square-free?
;; OEIS A005117
(loop with lst = '(1 2 3 5 6 7 10 11 13 14 15 17 19 21 22 23 26 29 30 31 33 34 35 37 38 39 41 42 43 46 47 51 53 55 57 58 59 61 62
65 66 67 69 70 71 73 74 77 78 79 82 83 85 86 87 89 91 93 94 95 97 101 102 103 105 106 107 109 110 111 113)
with c = (car lst)
for n from 1
while c
do (if (= n c)
(progn
(setq lst (cdr lst)
c (car lst))
(assert-true (mjr_prime_square-free? n))
(assert-true (mjr_prime_square-free? n)))))
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(run-tests
)
;; lucas primes. (sloan A005479)
;; Note: If Ln is prime, then n=0 || n is prime || n=2^m where m is 1,2,3,4
;;
;; (loop for n from 1
;; for l = (mjr_prime_lucas-number n)
;; until (> l mjr_prime::*mjr_prime_small-max*)
;; when (mjr_prime_primep-small l)
;; collect l)
;;
;; (3 7 11 29 47 199 521 2207 3571 9349 3010349 54018521)
;; lucas psuoprimes (composit n for which n|[L(n) - 1] (sloan A005845)
;;
;; (loop for n from 1 upto 7000
;; for ln = (mjr_prime_lucas-number n)
;; when (and (not (mjr_prime_primep-small n)) (zerop (rem (1- ln) n)))
;; collect n))
;;
;; (1 705 2465 2737 3745 4181 5777 6721)
;; ;; A000668 Mersenne primes
;; (loop for j from 0 upto 89
;; for p = (mjr_prime_nth-small-prime j)
;; for i = (1- (expt 2 p))
;; when (mjr_prime_probable-primep-bpsw i)
;; collect i)
;; ;; A000043 Mersenne exponents: primes p such that 2^p - 1 is prime.
;; '(2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9689 9941 11213 19937 21701 23209 44497 86243 110503
;; 132049 216091 756839 859433 1257787 1398269 2976221 3021377 6972593 13466917 20996011 24036583 25964951)