This repository has been archived by the owner on Aug 5, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 4
/
mp_prime.pas
2220 lines (1981 loc) · 68.3 KB
/
mp_prime.pas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
unit mp_prime;
{Basic 16/32 bit prime numbers and functions}
interface
{$i STD.INC}
{$ifdef VER70}
{$A+,N+}
{$endif}
uses
mp_types;
{$i mp_conf.inc}
(*************************************************************************
DESCRIPTION : Basic 16/32 bit prime numbers and functions
REQUIREMENTS : BP7, D1-D7/D9-D10/D12/D17-D18, FPC, VP
EXTERNAL DATA : mp_prm16.inc, mp_pbits.inc
MEMORY USAGE : 4 KB pbits, 13KB Primes16
heap for prime sieve routines
DISPLAY MODE : ---
REFERENCES : [3] Knuth, D.E.: The Art of computer programming. Vol 2
Seminumerical Algorithms, 3rd ed., 1998
[4] Forster, O.: Algorithmische Zahlentheorie, 1996
[5] (HAC) Menezes,A., von Oorschot,P., Vanstone, S: Handbook of
Applied Cryptography, 1996, www.cacr.math.uwaterloo.ca/hac
[7] P. Ribenboim: The New Book of Prime Number Records, 3rd ed., 1995.
[8] Marcel Martin: NX - Numerics library of multiprecision
numbers for Delphi and Free Pascal, 2006-2009
www.ellipsa.eu/public/nx/index.html
[10] R.Crandall, C.Pomerance: Prime Numbers, A Computational
Perspective, 2nd ed., 2005
[12] PARI/GP at http://pari.math.u-bordeaux.fr/
[24] H. Cohen, A Course in Computational Algebraic Number Theory
4th printing, 2000
[39] D.H. Lehmer: On the exact number of primes less than a given limit,
Illinois Journal of Mathematics, vol. 3, pp. 381-388, 1959;
available from http://projecteuclid.org/euclid.ijm/1255455259
[40] H. Riesel, Prime Numbers and Computer Methods for Factorization,
Vol. 126 of Progress in Mathematics, Boston, 2nd ed. 1994.
Paperback reprint 2012 in Modern Birkh„user Classics Series.
Version Date Author Modification
------- -------- ------- ------------------------------------------
1.9.00 26.12.08 W.Ehrhardt Initial version: IsPrime16/32,is_spsp32,is_spsp32A from mp_base
1.9.01 29.12.08 we prime_sieve routines
1.9.02 30.12.08 we prime_sieve with mp_getmem, prime_sieve_reset
1.9.03 04.01.09 we Primes16 array as global data, Primes16Index
1.9.04 04.01.09 we 32 bit first/next prime routines from mp_numth
1.12.00 05.07.09 we Removed prime residue classes mod 30
1.20.00 14.01.12 we BIT64: _spsp32, is_spsp32A, IsPrime32, modnpd2
1.20.01 21.01.12 we IsPrime16 inline if supported
1.21.00 08.07.12 we Assert psieve<>nil and paux<>nil in prime_sieve_next
1.21.01 14.07.12 we lsumphi32 and primepi32
1.21.02 15.07.12 we improved Primes16Index
1.21.03 23.07.12 we TFactors32 and PrimeFactor32
1.21.04 23.07.12 we EulerPhi32 and Carmichael32
1.21.05 24.07.12 we Moebius32
1.21.06 25.07.12 we primroot32
1.21.07 27.07.12 we order32
1.22.00 06.08.12 we SIEVE_MAXPRIME = MaxLongint if not MPC_SmallSieve
1.22.01 07.08.12 we prime32
1.22.02 03.09.12 we is_fundamental32
1.22.03 04.09.12 we is_squarefree32, improved Moebius32
1.22.04 05.09.12 we primepi32 uses icbrt32
1.22.05 05.09.12 we core32, quaddisc32, rad32
1.24.00 04.01.13 we BIT64 changed to PurePascal
1.29.00 14.07.14 we tau32
1.29.01 15.07.14 we safeprime32
1.30.00 28.09.14 we is_primepower32
1.30.01 30.09.14 we dlog32
1.30.02 01.10.14 we dlog32_ex, expanded and tuned version
1.30.03 03.10.14 we is_primroot32
1.30.04 04.10.14 we dlog32/ex without safe prime requirement
1.33.00 08.09.16 we is_Carmichael32
1.38.00 25.06.18 we prime32 uses doubles
**************************************************************************)
(*-------------------------------------------------------------------------
(C) Copyright 2004-2018 Wolfgang Ehrhardt
This software is provided 'as-is', without any express or implied warranty.
In no event will the authors be held liable for any damages arising from
the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software in
a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
----------------------------------------------------------------------------*)
{#Z+}
{---------------------------------------------------------------------------}
{------------------------ Basic prime functions ----------------------------}
{---------------------------------------------------------------------------}
{#Z-}
const
NumPrimes16 = 6542; {Number of 16 bit primes}
const
Primes16: array[1..NumPrimes16] of word {array of all 16 bit primes}
= {#Z+}({$i mp_prm16.inc}){#Z-};
const
pbits16: array[0..4095] of byte = {#Z+}({$i mp_pbits.inc}){#Z-};
{Bit array for 16 bit primes. Only odd integers are represented}
pmask16: array[0..15] of byte = (0,1,0,2,0,4,0,8,0,16,0,32,0,64,0,128);
{Mask array: odd i is prime if pbits16[i shr 4] and pmask16[i and $f] <> 0)}
function IsPrime16(N: word): boolean; {$ifdef HAS_INLINE} inline;{$endif}
{-Test if N is prime}
function Primes16Index(n: word): integer;
{-Return index of largest prime <= n in Primes16 array; 1 if n<=2}
{ Since Primes16[1]=2, this is identical to primepi(n) for n>=2. }
function IsPrime32(N: longint): boolean;
{-Test if longint (DWORD) N is prime}
function is_primepower32(n: longint; var p: longint): integer;
{-Test if n is a prime power: return 0 if not, n = p^result otherwise.}
{ Note: contrary to mp_is_primepower a prime n is returned as n = n^1!}
function lsumphi32(x: longint; a: integer): longint;
{-Return the partial sieve function Phi(x,a), the number of integers in}
{ (0,x] which are co-prime to the first a primes, aka 'Legendre's sum'}
function primepi32(x: longint): longint;
{-Return the prime counting function pi(x) using Lehmer's formula}
function is_spsp32(N, A: longint): boolean;
{-Strong probable prime test for N with base A, calls is_spsp32A}
function is_spsp32A(N: longint; const bases: array of longint): boolean;
{-Strong probable prime test for N with a number of bases. }
{ Return true if N is a probable prime for all bases, or if}
{ N=2. Negative numbers are treated as unsigned (=cardinal)}
{#Z+}
{---------------------------------------------------------------------------}
{--------------------- 32-bit number-theoretic functions -------------------}
{---------------------------------------------------------------------------}
{#Z-}
type
{Record for prime factorization, 10 primes are enough for 32 bit numbers,}
{because numbers with 10 different prime factors are >= 6469693230 > 2^32}
{number = product(primes[i]^pexpo[i], i=1..pcount), if n > 1}
TFactors32 = record
number: longint; {original number}
primes: array[1..10] of longint; {prime factors in ascending order}
pexpo : array[1..10] of byte; {exponents of prime factors}
pcount: integer; {number of different primes}
end;
function Carmichael32(n: longint): longint;
{-Return the Carmichael function lambda(|n|), lambda(0)=0. For n > 0 this}
{ is the least k such that x^k = 1 (mod n) for all x with gcd(x,n) = 1. }
function core32(n: longint): longint;
{-Return the squarefree part of n, n<>0, i.e. the unique squarefree integer c with n=c*f^2}
function dlog32(a,b,p: longint): longint;
{-Compute the discrete log_a(b) mod p using Pollard's rho algorithm: i.e.}
{ solve a^x = b mod p, with p prime, a > 1, b > 0; return -1 for failure.}
{ If a is no primitive root mod p, solutions may not exist or be unique. }
function dlog32_ex(a,b,p,JMAX: longint): longint;
{-Compute the discrete log_a(b) mod p using Pollard's rho }
{ algorithm; p prime, a > 1, b > 0; return -1 for failure. }
{ Expanded version with variable trial parameter JMAX >= 0.}
function EulerPhi32(n: longint): longint;
{-Return Euler's totient function phi(|n|), phi(0)=0. For n > 0 }
{ this the number of positive integers k <= n with gcd(n,k) = 1.}
function is_Carmichael32(n: longint): boolean;
{-Test if |n| is a Carmichael number}
function is_fundamental32(d: longint): boolean;
{-Return true, if d is a fundamental discriminant (either d=1 mod 4 and }
{ squarefree, or d=0 mod 4, d/4 = 2,3 mod and squarefree), false if not.}
function is_primroot32(const g,n: longint): boolean;
{-Test if g is primitive root mod n}
function is_squarefree32(n: longint): boolean;
{-Return true if n is squarefree, i.e. not divisible by a square > 1}
function Moebius32(n: longint): integer;
{-Return the Moebius function mu(abs(n)), mu(0)=0, mu(1)=1. mu(n)=(-1)^k }
{ if n > 1 is the product of k different primes; otherwise mu(n)=0. }
function order32(n,m: longint): longint;
{-Return the order of n mod m, m > 1, i.e. the smallest integer}
{ e with n^e = 1 mod m; if gcd(n,m) <> 1 the result is 0.}
function prime32(k: longint): longint;
{-Return the kth prime if 1 <= k <= 105097565, 0 otherwise}
procedure PrimeFactor32(n: longint; var FR: TFactors32);
{-Return the prime factorization of n > 1, FR.pcount=0 if n < 2}
function primroot32(n: longint): longint;
{-Compute the smallest primitive root mod n, 0 if n does not have a prim.root}
function quaddisc32(n: longint): longint;
{-Return the discriminant of the quadratic field Q(sqrt(n))}
function rad32(n: longint): longint;
{-Return the radical rad(n) of n = product of the distinct prime factors of n.}
function tau32(n: longint): integer;
{-Return the number of positive divisors of n (including 1 and n)}
{#Z+}
{---------------------------------------------------------------------------}
{------------------------ Prime sieve functions ----------------------------}
{---------------------------------------------------------------------------}
{#Z-}
type
TPrimeContext = record {Context for FindFirst/NextPrime32}
prime: longint; {last found prime}
next : longint; {next value to test}
index: integer; {index in residue class table}
end;
procedure FindFirstPrime32(n: longint; var ctx: TPrimeContext);
{-Find first prime >= n and initialize ctx for FindNextPrime32}
procedure FindNextPrime32(var ctx: TPrimeContext);
{-Find next 32 bit prime (DWORD interpretation, see note), success if ctx.prime<>0}
function nextprime32(n: longint): longint;
{-Next 32 bit prime >= n (DWORD interpretation, see note)}
procedure nextprime32_array(n,cmax: longint; var a: array of longint);
{-Fill an array with the next 32 bit primes >= n (DWORD interpretation).}
{ If a[i0]=4294967291 then a[i]=0 for i>i0. If cmax<=0 fill complete array.}
function prevprime32(n: longint): longint;
{-Previous 32 bit prime <= n, prevprime32(0)=0, (DWORD interpretation)}
function safeprime32(n: longint): longint;
{-Return the next safe prime p >= n, i.e. p and (p-1)/2 prime; 0 if n > 2147483579}
{$ifdef MPC_SmallSieve}
const
SIEVE_MAXPRIME = 1908867043; {maximum prime that will be generated, see note}
SIEVE_PRIMES = 4550; {= primepi(sqrt(SIEVE_MAXPRIME))-1}
SIEVE_BLOCK = 32768; {sieve size, may be set to 16384 for BIT16 if memory is low}
type
TAux = array[0..SIEVE_PRIMES-1] of word; {prime offsets in sieve}
TFlags = array[0..SIEVE_BLOCK-1] of boolean; {prime flags in sieve}
TSieve = record {prime sieve context}
paux : ^TAux; {pointer to offsets }
psieve : ^TFlags; {pointer to flags }
curr_blk: longint; {current sieve block}
curr_off: longint; {offset in curr_blk }
end;
{$else}
const
SIEVE_MAXPRIME = MaxLongint; {maximum prime that will be generated, see note}
SIEVE_PRIMES = 4792; {= primepi(sqrt(SIEVE_MAXPRIME))-1}
SIEVE_BLOCK = 32768; {sieve size, may be set to 16384 for BIT16 if memory is low}
type
TAux = array[0..SIEVE_PRIMES-1] of longint;{prime offsets in sieve}
TFlags = array[0..SIEVE_BLOCK-1] of boolean; {prime flags in sieve}
TSieve = record {prime sieve context}
paux : ^TAux; {pointer to offsets }
psieve : ^TFlags; {pointer to flags }
curr_blk: longint; {current sieve block}
curr_off: longint; {offset in curr_blk }
end;
{$endif}
procedure prime_sieve_clear(var sieve: TSieve);
{-Release memory allocated by prime_sieve_init}
procedure prime_sieve_init(var sieve: TSieve; first_prime: longint);
{-Allocate/initialize sieve to return primes >= first_prime; first_prime <= SIEVE_MAXPRIME}
function prime_sieve_next(var sieve: TSieve): longint;
{-Return next prime from sieve, 1 if done}
procedure prime_sieve_reset(var sieve: TSieve; first_prime: longint);
{-Initialize already allocated sieve to return primes >= first_prime; first_prime <= SIEVE_MAXPRIME}
{#Z+}
{Tables for prime residue classes mod 210, calculated with t_rcnp.pas}
{Must be interfaced because they are used by mp_nextprime/mp_prevprime.}
{NPRC_OddIdx has the number of the prime residue classes of an odd value }
{n mod NPRC_MOD, or -1 if there is no corresponding prime residue class. }
{NPRC_Diff contains the differences: NPRC_Diff[i] = prctab[i+1]-prctab[i]}
{prime residue classes mod 210. Cnt=48=1*2*4*6, numbered from 0 to 47}
{ 1, 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, 121, 127, 131, 137, 139, 143, 149, 151, 157,
163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209}
const
NPRC_MOD = 210;
NPRC_NRC = 48;
const
NPRC_OddIdx: array[0..(NPRC_MOD div 2)-1] of shortint =
( 0, -1, -1, -1, -1, 1, 2, -1, 3, 4, -1, 5, -1, -1, 6,
7, -1, -1, 8, -1, 9, 10, -1, 11, -1, -1, 12, -1, -1, 13,
14, -1, -1, 15, -1, 16, 17, -1, -1, 18, -1, 19, -1, -1, 20,
-1, -1, -1, 21, -1, 22, 23, -1, 24, 25, -1, 26, -1, -1, -1,
27, -1, -1, 28, -1, 29, -1, -1, 30, 31, -1, 32, -1, -1, 33,
34, -1, -1, 35, -1, -1, 36, -1, 37, 38, -1, 39, -1, -1, 40,
41, -1, -1, 42, -1, 43, 44, -1, 45, 46, -1, -1, -1, -1, 47);
const
NPRC_Diff: array[0..NPRC_NRC-1] of mp_digit =
(10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4,
2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4,
6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2);
{#Z-}
implementation
uses
{$ifdef MPC_PurePascal} BTypes,{$endif}
mp_base, mp_prng;
{---------------------------------------------------------------------------}
{------------------------ Basic prime functions ----------------------------}
{---------------------------------------------------------------------------}
{---------------------------------------------------------------------------}
function IsPrime16(N: word): boolean; {$ifdef HAS_INLINE} inline;{$endif}
{-Test if N is prime}
begin
IsPrime16 := (N=2) or (pbits16[N shr 4] and pmask16[N and $0F] <> 0);
end;
{---------------------------------------------------------------------------}
function Primes16Index(n: word): integer;
{-Return index of largest prime <= n in Primes16 array; 1 if n<=2}
{ Since Primes16[1]=2, this is identical to primepi(n) for n>=2. }
var
i,j,k,m: integer;
p: word;
const
js: array[0..32] of integer = (1,309,564,801,1028,1254,1469,1681,1900,2110,
2312,2517,2725,2918,3124,3314,3512,3716,3908,
4098,4288,4495,4678,4858,5051,5239,5432,5616,
5814,6003,6179,6363,6541);
begin
{Binary search in Primes16 array}
Primes16Index := 1; {Avoid FPC warning}
if n<=2 then exit
else if n>=65521 then Primes16Index := NumPrimes16
else begin
k := n shr 11;
{Get initial lower and upper bounds better than 1 and NumPrimes16}
i := js[k];
j := js[k+1];
m := 0;
repeat
k := (i+j) shr 1;
if k=m then break;
p := Primes16[k];
if n < p then j := k-1
else if n > p then i := k+1
else break;
m := k;
until false;
Primes16Index := k;
end;
end;
{$ifdef BIT16}
{$undef USEA5} {smaller but slower lsumphi32}
{$else}
{$define USEA5} {faster but larger lsumphi32}
{$endif}
{---------------------------------------------------------------------------}
function lsumphi32(x: longint; a: integer): longint;
{-Return the partial sieve function Phi(x,a), the number of integers in}
{ (0,x] which are co-prime to the first a primes, aka 'Legendre's sum'}
var
sum: longint; {accumulates contributions from recursive calls}
{$ifdef USEA5}
const
BaseA = 5;
ProdA = 2310;
PrA_2 = 1155;
PrA_4 = 578;
PhiPA = 480;
const
{compressed table of Phi(x mod 2310, 5), TabPhiA[i] = lphi(2*i,5)}
TabPhiA: array[0..PrA_4] of byte = (
0, 1,1,1,1,1,1,2,2,3,4,4,5,5,5,6,7,7,7,8,8,9,10,10,11,11,11,
12,12,12,13,14,14,14,15,15,16,17,17,17,18,18,19,19,19,20,20,
20,20,21,21,22,23,23,24,25,25,26,26,26,26,26,26,26,27,27,28,
28,28,29,30,30,30,30,30,31,32,32,32,33,33,33,34,34,35,36,36,
37,37,37,38,39,39,39,39,39,40,41,41,42,43,43,43,43,43,43,44,
44,44,44,44,45,46,46,47,48,48,49,49,49,50,51,51,51,52,52,53,
53,53,54,54,54,55,55,55,56,57,57,57,58,58,59,60,60,60,61,61,
62,62,62,63,63,63,63,64,64,65,66,66,67,67,67,68,68,68,68,69,
69,69,70,70,70,70,70,71,72,72,73,73,73,74,75,75,75,76,76,76,
77,77,78,79,79,80,80,80,81,82,82,82,83,83,84,85,85,85,86,86,
86,86,86,87,88,88,88,88,88,89,90,90,91,92,92,93,93,93,94,94,
94,94,95,95,96,97,97,98,98,98,98,98,98, 99, 100,100,100,101,
101,102,103,103,103,104,104,105,105,105,106,106,106,106,106,
106,107,108,108,109,110,110,111,111,111,111,112,112,112,113,
113,114,114,114,115,116,116,117,117,117,118,119,119,119,120,
120,120,120,120,121,122,122,123,123,123,124,125,125,125,126,
126,127,128,128,129,130,130,130,130,130,131,132,132,132,132,
132,133,134,134,135,135,135,136,136,136,137,138,138,138,139,
139,139,140,140,141,141,141,142,142,142,143,144,144,144,145,
145,146,147,147,147,148,148,149,149,149,150,150,150,150,151,
151,152,153,153,153,154,154,155,155,155,155,156,156,156,157,
157,158,158,158,159,160,160,161,161,161,162,162,162,162,163,
163,163,164,164,165,166,166,166,166,166,167,168,168,168,169,
169,170,171,171,172,173,173,173,173,173,174,175,175,175,175,
175,176,177,177,178,179,179,180,180,180,180,181,181,181,182,
182,183,184,184,185,185,185,186,186,186,187,188,188,188,189,
189,190,190,190,190,191,191,192,192,192,193,193,193,193,194,
194,195,196,196,197,198,198,199,199,199,199,200,200,200,201,
201,202,202,202,203,203,203,204,204,204,205,206,206,206,207,
207,207,208,208,209,210,210,211,211,211,212,213,213,213,214,
214,215,216,216,217,218,218,218,218,218,219,220,220,220,220,
220,221,222,222,222,223,223,224,224,224,225,226,226,226,227,
227,228,229,229,230,230,230,231,231,231,232,232,232,232,233,
233,234,235,235,235,236,236,236,236,236,237,237,237,237,238,
238,239,240,240);
{$else}
const
BaseA = 4;
ProdA = 210;
PrA_2 = 105;
PrA_4 = 53;
PhiPA = 48;
const
{compressed table of Phi(x mod 210, 4), TabPhiA[i] = lphi(2*i,4)}
TabPhiA: array[0..PrA_4] of shortint = (
0,1,1,1,1,1,2,3,3,4,5,5,6,6,6,7,8,8,8,9,9,10,11,11,
12,12,12,13,13,13,14,15,15,15,16,16,17,18,18,18,19,
19,20,20,20,21,21,21,21,22,22,23,24,24);
{$endif}
procedure phirec(x: longint; a, sign: integer);
{-Internal function for Phi(x,a), uses the accumulator variable sum}
begin
{This is an optimized Pascal implementation of the C snippet given in}
{Table II of T. Oliveira e Silva, Computing pi(x): the combinatorial }
{method, Revista do DETUA, vol. 4, no. 6, pp. 759-768, March 2006. }
{Available online from http://www.ieeta.pt/~tos/bib/5.4.pdf }
repeat
{T. Oliveira e Silva uses recursion downto a=0 in the demo, I stop }
{at a=5 or a=4, and use a second function for smaller values of a. }
if a=BaseA then begin
{Compute Phi(x,BaseA) using the division and symmetry properties }
{of Phi(x,a), see Lehmer[39], p.384 or Riesel[40], p.14-17 }
a := x mod ProdA;
x := (x div ProdA)*PhiPA;
if a<PrA_2 then inc(x, integer(TabPhiA[succ(a) shr 1]))
else inc(x, PhiPA - integer(TabPhiA[(ProdA-a) shr 1]));
if sign<0 then dec(sum,x) else inc(sum,x);
exit;
end
else if x < primes16[a+1] then begin
inc(sum,sign);
exit;
end
else begin
{Use recursive property of Phi}
dec(a);
phirec(x div primes16[a+1], a, -sign);
end;
until false;
end;
function slowphi(x: longint; a: integer): longint;
{-Internal function for weird cases when Phi is called with a<BaseA}
begin
{This is just the recursive property of Phi}
if a<=0 then slowphi := x
else slowphi := slowphi(x, a-1) - slowphi(x div Primes16[a], a-1);
end;
begin
if x<0 then lsumphi32 := -lsumphi32(-x,a)
else if a<BaseA then lsumphi32 := slowphi(x,a)
else begin
sum := 0;
phirec(x,a,1);
lsumphi32 := sum;
end;
end;
{---------------------------------------------------------------------------}
function primepi32(x: longint): longint;
{-Return the prime counting function pi(x) using Lehmer's formula}
var
sum,z: longint;
a,b,c,i,j,k: integer;
begin
if x <= $FFFF then begin
if x < 2 then primepi32 := 0
else primepi32 := Primes16Index(word(x));
exit;
end;
z := isqrt32(x);
b := Primes16Index(z); {b = pi(x^(1/2) <= 4792}
c := Primes16Index(icbrt32(x)); {c = pi(x^(1/3) <= 209}
a := Primes16Index(isqrt32(z)); {a = pi(x^(1/4) <= 47}
sum := longint(b+a-2)*(b-a+1) shr 1;
for i:=a+1 to c do begin
z := x div Primes16[i];
sum := sum - primepi32(z);
k := Primes16Index(isqrt32(z));
for j:=i to k do begin
{z div Primes16[j] <= x / Primes16[a+1]^2 <= x^(1/2) < 2^16}
sum := sum - (Primes16Index(z div Primes16[j]) - j + 1);
end;
end;
for i:=c+1 to b do sum := sum - primepi32(x div Primes16[i]);
primepi32 := lsumphi32(x,a) + sum;
end;
{$ifdef MPC_PurePascal}
{---------------------------------------------------------------------------}
function _spsp32(a,N,d,s: uint32): boolean;
{-Strong probable prime test for N with base a, N-1=2^s*d, s>0}
var
dk,p,N1: uint32;
i: word;
begin
_spsp32 := true;
N1 := N-1;
{d>0 is odd, so first iteration can be done manually without MulMod}
p := a;
dk := d shr 1;
while dk<>0 do begin
a := uint64(a)*a mod N;
if odd(dk) then p := uint64(p)*a mod N;
dk := dk shr 1;
end;
{here p=bases[k]^d mod N}
if (p=1) or (p=N1) then exit;
{calculate p^i for i<s, note shifted for loop variable}
for i:=2 to s do begin
p := uint64(p)*p mod N;
if p=N1 then exit;
end;
_spsp32 := false;
end;
{-----------------------------------------------------------------------------}
function is_spsp32A(N: longint; const bases: array of longint): boolean;
{-Strong probable prime test for N with a number of bases. }
{ Return true if N is a probable prime for all bases, or if}
{ N=2. Negative numbers are treated as unsigned (=cardinal)}
var
s,d,N1: uint32;
k: integer;
begin
N1 := uint32(N)-1;
if odd(N1) or (N1=0) then begin
is_spsp32A := N=2;
exit;
end;
{calculate d*2^s = N-1, this is done once for all bases}
d := N1;
s := 0;
while not odd(d) do begin
d := d shr 1;
inc(s);
end;
{default no spsp}
is_spsp32A := false;
{now loop through all the bases}
for k:=low(bases) to high(bases) do begin
{spsp test for bases[k]}
if not _spsp32(uint32(bases[k]),uint32(N),d,s) then exit;
end;
{All tests passed, this is a spsp}
is_spsp32A := true;
end;
{$else}
{$ifdef BIT32}
{---------------------------------------------------------------------------}
function _spsp32(a,N,d,s: longint): boolean;
{-Strong probable prime test for N with base a, N-1=2^s*d, s>0}
var
res: boolean;
begin
res := true;
asm
pushad
mov esi,[N] {esi=N}
mov ecx,[a] {ecx= a mod N}
cmp ecx,esi
jb @@0
mov eax,ecx {reduce a mod N}
sub edx,edx {to avoid overflow}
div esi
mov ecx,edx {ecx=a mod N}
@@0: mov ebx,ecx {ebx=p=a}
mov edi,[d] {edi=dk}
shr edi,1 {dk=dk shr 1}
{while dk<>0 do }
@@1: or edi,edi
jz @@3
mov eax,ecx {a=a^2 mod N}
mul eax
div esi
mov ecx,edx
shr edi,1 {dk=dk shr 1}
jnc @@1
mov eax,ebx {p=p*a mod N}
mul ecx
div esi
mov ebx,edx
jmp @@1
{end while}
@@3: dec ebx {if p=1 then goto done}
jz @@5
inc ebx
mov edi,[N]
dec edi {edi=N1}
cmp ebx,edi {if p=N1 then goto done}
jz @@5
mov ecx,s {remember: s > 0}
dec ecx {if s=1 then N not prime}
jz @@np
{for i=2 to s do}
@@4: mov eax,ebx {p=p^2 mod N}
mul ebx
div esi
mov ebx,edx
cmp ebx,edi {if p=N1 then goto done}
jz @@5
dec ecx
jnz @@4
{end for}
@@np: mov [res],cl {not prime, here cl=0!}
@@5:
popad
end;
_spsp32 := res;
end;
{$else}
{---------------------------------------------------------------------------}
function _spsp32(a,N,d,s: longint): boolean;
{-Strong probable prime test for N with base a, N-1=2^s*d, s>0}
var
res: boolean;
begin
res := true;
asm
db $66; mov si,word ptr [N] {si=N}
db $66; mov cx,word ptr [a] {cx=a mod N}
db $66; cmp cx,si
jb @@0
db $66; mov ax,cx
db $66; sub dx,dx {to avoid overflow}
db $66; div si
db $66; mov cx,dx {cx=a mod N fixed in V1.2.12}
@@0: db $66; mov bx,cx {bx=p=a}
db $66; mov di,word ptr [d] {di=dk}
db $66; shr di,1 {dk=dk shr 1}
{while dk<>0 do }
@@1: db $66; or di,di
jz @@3
db $66; mov ax,cx {a=a^2 mod N}
db $66; mul ax
db $66; div si
db $66; mov cx,dx
db $66; shr di,1 {dk=dk shr 1}
jnc @@1
db $66; mov ax,bx {p=p*a mod N}
db $66; mul cx
db $66; div si
db $66; mov bx,dx
jmp @@1
{end while}
@@3: db $66; dec bx {if p=1 then goto done}
jz @@5
db $66; inc bx
db $66; mov di,word ptr [N] {di=N1}
db $66; dec di
db $66; cmp bx,di {if p=N1 then goto done}
jz @@5
mov cx,word ptr [s] {remember: s > 0}
dec cx {if s=1 then N not prime}
jz @@np
{for i=2 to s do}
@@4: db $66; mov ax,bx {p=p^2 mod N}
db $66; mul bx
db $66; div si
db $66; mov bx,dx
db $66; cmp bx,di {if p=N1 then goto done}
jz @@5
dec cx
jnz @@4
{end for}
@@np: mov [res],cl {not prime, here cl=0!}
@@5:
end;
_spsp32 := res;
end;
{$endif}
{-----------------------------------------------------------------------------}
function is_spsp32A(N: longint; const bases: array of longint): boolean;
{-Strong probable prime test for N with a number of bases. }
{ Return true if N is a probable prime for all bases, or if}
{ N=2. Negative numbers are treated as unsigned (=cardinal)}
var
s,d,N1: longint;
k: integer;
begin
N1 := N-1;
if odd(N1) or (N1=0) then begin
is_spsp32A := N=2;
exit;
end;
{calculate d*2^s = N-1, this is done once for all bases}
{$ifdef BIT16}
asm
db $66; mov ax,word ptr [N1]
db $66, $0F,$BC,$C8 {bsf ecx,eax}
db $66; mov word ptr [s],cx
db $66; shr ax,cl
db $66; mov word ptr [d],ax
end;
{$else}
asm
mov eax,[N1]
bsf ecx,eax
mov [s],ecx
shr eax,cl
mov [d],eax
end;
{$endif}
{default no spsp}
is_spsp32A := false;
{now loop through all the bases}
for k:=low(bases) to high(bases) do begin
{spsp test for bases[k]}
if not _spsp32(bases[k],N,d,s) then exit;
end;
{All tests passed, this is a spsp}
is_spsp32A := true;
end;
{$endif}
{---------------------------------------------------------------------------}
function IsPrime32(N: longint): boolean;
{-Test if longint (DWORD) N is prime}
const
a1: array[0..1] of longint = (31,73);
a2: array[0..2] of longint = (2,7,61);
type
LH = packed record L,H: word; end;
begin
if LH(N).H=0 then begin
{$ifndef BIT16}
IsPrime32 := (N=2) or (pbits16[N shr 4] and pmask16[N and $0F] <> 0);
{$else}
IsPrime32 := (word(N)=2) or (pbits16[word(N) shr 4] and pmask16[word(N) and $0F] <> 0);
{$endif}
end
else begin
{Normally here N>=2^16, IN ANY CASE: N MUST NOT BE IN [2,7,31,61,73]!}
if odd(N) then begin
{First test N and $80000000 <> 0 selects N>MaxLongint,}
{second test is the upper bound for spsp set (31,73)] }
if (N and $80000000 <> 0) or (N>=9080191) then IsPrime32 := is_spsp32A(N, a2)
else IsPrime32 := is_spsp32A(N, a1);
end
else begin
{N is even, N<>2 therefore N not prime}
IsPrime32 := false;
end;
end;
end;
{---------------------------------------------------------------------------}
function is_primepower32(n: longint; var p: longint): integer;
{-Test if n is a prime power: return 0 if not, n = p^result otherwise.}
{ Note: contrary to mp_is_primepower a prime n is returned as n = n^1!}
var
i,r: integer;
begin
is_primepower32 := 0;
if n>1 then begin
if n and 1 = 0 then begin
{even n can only be powers of 2}
if n and (n-1) = 0 then begin
p := 2;
is_primepower32 := bitsize32(n)-1;
end;
end
else if isprime32(n) then begin
p := n;
is_primepower32 := 1;
end
else begin
for i:=2 to NumPrimes16 do begin
p := Primes16[i];
if n mod p = 0 then begin
r := 1;
n := n div p;
while (n >= p) and (n mod p = 0) do begin
n := n div p;
inc(r);
end;
if n=1 then is_primepower32 := r;
exit;
end;
end;
end;
end;
end;
{---------------------------------------------------------------------------}
function is_spsp32(N, A: longint): boolean;
{-Strong probable prime test for N with base A, calls is_spsp32A}
begin
is_spsp32 := is_spsp32A(N, A);
end;
{$ifdef MPC_PurePascal}
{------------------------------------------------------------}
function modnpd2(n: longint): integer; {$ifdef HAS_INLINE} inline;{$endif}
{-Calculate (n mod NPRC_MOD) div 2 (DWORD interpretation)}
begin
modnpd2 := (uint32(n) mod NPRC_MOD) shr 1;
end;
{$else}
{$ifdef BIT32}
{------------------------------------------------------------}
function modnpd2(n: longint): integer; assembler; {&frame-}
{-Calculate (n mod NPRC_MOD) div 2 (DWORD interpretation)}
asm
{$ifdef LoadArgs}
mov eax,[n]
{$endif}
mov ecx,NPRC_MOD
sub edx,edx
div ecx
mov eax,edx
shr eax,1
end;
{$else}
{------------------------------------------------------------}
function modnpd2(n: longint): integer; assembler;
{-Calculate (n mod NPRC_MOD) div 2 (DWORD interpretation)}
asm
db $66; mov ax, word ptr [n]
db $66; sub dx,dx
db $66; mov cx,NPRC_MOD; dw 0;
db $66; div cx
mov ax,dx
shr ax,1
end;
{$endif}
{$endif}
{---------------------------------------------------------------------------}
function nextprime32(n: longint): longint;
{-Next 32 bit prime >= n (DWORD interpretation, see note)}
var
id,k: integer;
begin
{easy outs and assure valid range for mod MP_MOD calculation}
{note: (n>=-4) CANNOT be omitted because of DWORD interpretation}
{for m=-4=$FFFFFFFC=4294967292, nextprime(m) is greater 2^32 and}
{will be set to 0 (nextprime(m) mod 2^32 = 15 would be strange!)}
if (n>=-4) and (n<=7) then begin
if n<0 then nextprime32 := 0
else if n<=2 then nextprime32 := 2
else if n<=3 then nextprime32 := 3
else if n<=5 then nextprime32 := 5
else nextprime32 := 7;
exit;
end;
{make n odd}
if n and 1 = 0 then inc(n);
{$ifndef BIT16}
{avoid warning, id WILL always be initialized (for bug-free modnpd2!)}
id := 0;
{$endif}
{move n to next prime residue class mod MP_MOD and index id into diff array}
for k:=modnpd2(n) to (NPRC_MOD div 2)-1 do begin
id := NPRC_OddIdx[k];
{note: loop terminates via break because NPRC_OddIdx[(NPRC_MOD div 2)-1]<>-1}
if id<>-1 then break;
inc(n,2);
end;
repeat
{loop through possible primes}
if IsPrime32(n) then begin
nextprime32 := n;
exit;
end;
{move to next candidate}
inc(n,longint(NPRC_Diff[id]));
{get next increment index}
inc(id); if id>=NPRC_NRC then id:=0;
until false;
end;
{---------------------------------------------------------------------------}
procedure FindFirstPrime32(n: longint; var ctx: TPrimeContext);
{-Find first prime >= n and initialize ctx for FindNextPrime32}
begin
with ctx do begin
next := n;
index := -1;
end;
FindNextPrime32(ctx);
end;
{---------------------------------------------------------------------------}
procedure FindNextPrime32(var ctx: TPrimeContext);
{-Find next 32 bit prime (DWORD interpretation, see note), success if ctx.prime<>0}
var
k: integer;
found: boolean;
const
MP32 = longint($FFFFFFFB); {4294967291 = prevprime(2^32)}