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_numth.pas
8414 lines (7518 loc) · 247 KB
/
mp_numth.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_numth;
{MPArith number theoretic functions}
interface
{$i STD.INC}
{$ifndef FPC}
{$N+}
{$endif}
uses
mp_types;
{$i mp_conf.inc}
(*************************************************************************
DESCRIPTION : MPArith number theoretic functions
REQUIREMENTS : BP7, D2-D7/D9-D10/D12/D17-D18/D25S, FPC, VP, WDOSX
EXTERNAL DATA : (mp_types)
MEMORY USAGE : heap
DISPLAY MODE : ---
REFERENCES : [1] LibTomMath 0.30+ by Tom St Denis
[2] MPI by M.J. Fromberger
[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
[6] R. P. Brent, Factor: an integer factorization program for
the IBM PC, Report TR-CS-89-23, October 1989, 7 pp.
http://maths-people.anu.edu.au/~brent/pub/pub117.html
http://maths-people.anu.edu.au/~brent/ftp/rpb117/rpb117.exe
[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
[9] Wei Dai: Lucas Sequences in Cryptography,
http://www.weidai.com/lucas.html
[10] Crandall,R., C.Pomerance: Prime Numbers, A Computational
Perspective, 2nd ed., 2005
[11] Peter Luschny's Homepage of Factorial Algorithms,
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
[12] PARI/GP at http://pari.math.u-bordeaux.fr/
[15] The GNU Multiple Precision Arithmetic Library, http://gmplib.org/
[24] H. Cohen, A Course in Computational Algebraic Number Theory
4th printing, 2000
[25] IEEE P1363/Draft. Standard Specifications for Public Key Cryptography.
Annex A (informative). Number-Theoretic Background.
[26] A. Adler, J.E. Coury: The Theory of Numbers, 1995
[32] S.C. Lindhurst, Computing Roots in Finite Fields and Groups, with a
Jaunt through Sums of Digits, University of Wisconsin, Madison 1997
http://scott.lindhurst.com/papers/thesis.ps.gz
[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
------- -------- ------- ------------------------------------------
0.0.01 09.08.04 W.Ehrhardt Initial version: mp_addmod, mp_submod
0.0.02 09.08.04 we mp_mulmod, mp_sqrmod
0.0.03 09.08.04 we raw versions of mp_pollard_rho, mp_gcd
0.0.04 09.08.04 we mp_pollard_rho with Barrett reduction
0.0.05 09.08.04 we mp_gcd (binary version), mp_exteuclid (xgcd)
0.0.06 10.08.04 we mp_pollard_rho: Barrett after x*x+2
0.0.07 11.08.04 we mp_lcm
0.0.08 11.08.04 we mp_invmod using xgcd
0.0.08 12.08.04 we mp_invmodb using binary gcd
0.0.09 12.08.04 we fast_mp_invmod, inv_mod -> inv_mod_euclid, inv_modb -> inv_mod
0.0.10 12.08.04 we code cleanup
0.0.11 13.08.04 we s_mp_exptmod
0.0.12 13.08.04 we basic mp_exptmod
0.0.13 14.08.04 we WEXP_ constants
0.0.14 14.08.04 we mp_pollard_rho with cnt parameter
0.0.15 17.08.04 we mp_is_square
0.0.16 17.08.04 we mp_small_factor
0.0.17 17.08.04 we mp_small_factor with mod 30 increments
0.0.18 18.08.04 we mp_is_SPP
0.0.19 22.08.04 we small prime tables
0.0.20 22.08.04 we mp_miller_rabin
0.0.21 23.08.04 we mp_show_progress variable, mp_small_factor with f0
0.0.22 23.08.04 we Prime127 removed, reordering of functions
0.0.23 24.08.04 we s_mp_ppexpo, mp_pollard_pm1
0.0.24 24.08.04 we code cleanup, updated references
0.3.00 26.08.04 we "release" version 0.3
0.3.01 26.08.04 we mp_error, most functions now procedures
0.3.02 26.08.04 we mp_progress, mp_set_progress
0.3.03 16.02.05 we Montgomery Reduction
0.3.04 16.02.05 we temp. $R- for mp_montgomery_setup (for 8-bit)
0.3.05 17.02.05 we general mp_exptmod_win, removed s_mp_exptmod
0.3.06 20.02.05 we mp_jacobi
0.3.07 20.02.05 we mp_jacobi: fix up D6/FPC warnings
0.3.08 27.03.05 we Montgomery routines to mp_core
0.3.09 29.03.05 we mp_jacobi: non-recursive version
0.3.10 08.05.05 we mp_jacobi: clean up, remove old version
0.3.11 16.05.05 we FPC: $mode objfpc/$goto on
0.3.12 10.08.05 we mp_lucasv.., debug check in mp_xgcd
0.3.13 15.08.05 we bugfix mp_lucasv2p if v same var as p or q
0.3.13 15.08.05 we mp_williams_pp1: William's p+1 method
0.3.15 16.08.05 we mp_coshmult compatible with Forster
0.3.16 17.08.05 we s_mp_coshmult easy out for b=1, loop error check
0.3.17 18.08.05 we mp_williams_pp1: numtest parameter
0.4.00 20.08.05 we use mp_set_error
0.4.01 20.08.05 we longint bit count in mp_exptmod_win
0.4.02 21.08.05 we MPC_ArgCheck
0.4.03 22.08.05 we removed mp_argcheck, mp_errchk
0.4.04 22.08.05 we now functions: mp_is_SPP, mp_is_square
0.4.05 22.08.05 we mp_fermat
0.4.06 22.08.05 we pollard, williams, coshmult with init<i>
0.4.07 23.08.05 we mp_mersenne, MPC_HaltOnArgCheck
0.4.08 23.08.05 we new values for MaxFermat, MaxMersenne
0.4.09 23.08.05 we IsPrime15 (quick and dirty prime test)
0.4.10 23.08.05 we mp_exptmod_d
0.4.11 24.08.05 we Fibonacci numbers, mp_fib, mp_fib2
0.4.12 25.08.05 we Lucas numbers, mp_lucas, mp_lucas2
0.4.13 25.08.05 we mp_fib2 for n=0
0.4.14 26.08.05 we miller_rabin: generate a in the range 2<=a<n-1
0.4.15 27.08.05 we easy outs in mp_small_factor
0.4.16 28.08.05 we usage of exceptions implemented
0.4.17 10.09.05 we function IsMersennePrime
0.4.18 12.09.05 we mp_fact, MaxFact
0.4.19 12.09.05 we optimized 3-level mp_fact, FSIVec
0.4.20 13.09.05 we mp_fact optimization (balanced factors)
0.4.22 16.09.05 we changed mp_gcd easy out, comment in mp_xgcd
0.4.23 17.09.05 we mp_ecm_brent (Brent's ECM factoring method)
0.4.24 18.09.05 we mp_ecpwr (faster version of mp_ecpwr2)
0.4.25 20.09.05 we mp_lucasv2p without goto, changed TProgressProc
0.4.25 22.09.05 we $argcheck for pointers (mp_lucasv2p, mp_xgcd)
0.4.26 25.09.05 we Bugfix: k now longint in IsMersennePrime
0.4.27 26.09.05 we $argcheck for identical args in mp_fib2,mp_lucas2
0.5.00 29.09.05 we 'internal' functions moved to end of interface
0.5.01 30.09.05 we IsPrime16, pbits16, pmask16
0.5.02 01.10.05 we easy out for IsMersennePrime
0.5.03 02.10.05 we changed 'SPP' to more conventional 'spsp'
0.5.04 03.10.05 we function nextprime32
0.5.05 03.10.05 we function prevprime32
0.5.06 09.10.05 we next/prev prime residues classes module via $ifdef
0.5.07 09.10.05 we mp_is_spsp_d, mp_is_pprime, mp_small_factor with fmax
0.5.08 09.10.05 we small changes in mp_miller_rabin
0.5.08 13.10.05 we mp_is_slpsp
0.5.09 16.10.05 we mp_is_pprime with BPSW probable prime
0.5.10 14.11.05 we mp_nextprime/mp_prevprime
0.5.11 15.11.05 we mp_rand_prime
0.5.12 18.11.05 we mp_rand_prime improved for safe primes
0.5.13 20.11.05 we TPrimeType = (pt_normal, pt_3mod4, pt_safe)
0.5.14 20.11.05 we mp_is_pprime_ex, mp_rand_prime with mp_is_pprime_ex
0.5.15 21.11.05 we mp_rand_prime with mp_rand_bits
0.5.16 22.11.05 we New name: mp_isMersennePrime, mp_prng removed from uses
0.6.00 30.12.05 we mp_count_bits/CountBits32 renamed to mp_bitsize/bitsize32
0.6.01 30.12.05 we MP_8BIT removed
0.6.02 05.01.06 we MaxMersenne = DIGIT_BIT*MAXDigits-1;
0.6.03 10.01.06 we mp_fact using Recursive Split
0.7.00 28.04.06 we mp_xgcd: u2,v2,t2 suppressed
0.7.01 04.08.06 we TRedType MR_Reduce2k, mp_reduce_2k in mp_exptmod_win
0.7.03 05.08.06 we mp_show_progress initialized with false
0.7.04 08.08.06 we mp_makeodd in mp_miller_rabin, mp_is_slpsp, mp_jacobi
0.7.05 08.08.06 we mp_sqrtmod initial version
0.7.06 09.08.06 we mp_sqrtmod arg/init checks
0.7.07 09.08.06 we mp_lucasvmod
0.7.08 09.08.06 we mp_sqrtmod18_lucas, TPrimeType with pt_1mod8
0.7.09 10.08.06 we mp_sqrtmodp2, mp_sqrtmethod
0.7.10 11.08.06 we mp_miller_rabin and mp_is_spsp modified
0.7.11 11.08.06 we Avoid FPC warnings: nextprime32/prevprime32
0.7.12 11.08.06 we experimental mp_sqrtmod18_fp
0.7.13 11.08.06 we check prime/q-residue after 100 trials in mp_sqrtmod18_shanks/fp
0.7.14 15.08.06 we small tweak in mp_isMersennePrime
0.7.15 15.08.06 we $ifdef MPC_sqrtmod_fp
0.7.16 16.08.06 we bigprime_pm1: big prime stage for Pollard p-1
0.7.17 17.08.06 we bigprime_pm1: Barret and bitsize dependent constants
0.7.17 20.08.06 we const mp_max_small / mp_small_factor
0.7.18 21.08.06 we mp_ecm_factor (first complete version using 7KB stack)
0.7.19 22.08.06 we mp_ecm_factor (dynamic vectors, progress, intermeditate GCDs)
0.7.20 23.08.06 we mp_ecm_factor,mp_ecm_brent: check n>1, @n<>@f
0.7.21 23.08.06 we s_mp_ppexpo with mp_mul_w
0.7.22 23.08.06 we mp_ecm_factor with Barrett: speed more than doubled
0.7.23 24.08.06 we mp_ecm_factor: check seed>5 in ecm_setup, more comments
0.7.24 26.08.06 we nextprime32_array
0.7.25 27.08.06 we Prime table in mp_ecm_factor $ifdef MPC_ECM_Primetable
0.7.26 28.08.06 we basic mp_is_power
0.7.27 07.09.06 we mp_is_square2
0.7.28 07.09.06 we nextprime32_array with parameter cmax
0.7.29 07.09.06 we mp_ecm_factor with parameters C1,phase; mp_ecm_simple
0.7.30 07.09.06 we function ProgressAssigned
0.8.00 19.09.06 we mp_primorial, MaxPrimorial
0.8.01 19.09.06 we mp_isMersennePrime: search small factor before Lucas/Lehmer
0.8.02 23.09.06 we MaxFibonacci, MaxLucas
0.8.03 25.09.06 we mp_kronecker, mp_jacobi uses kron_intern
0.9.00 25.12.06 we minor changes in mp_isMersennePrime
0.9.01 28.12.06 we Bugfix mp_is_power for a=-4
0.9.02 30.12.06 we FindFirst/NextPrime32, rewrite nextprime32_array
0.9.03 01.01.07 we bigalloc renamed to IAlloc
0.9.04 01.01.07 we mp_primorial with mp_prod_int
0.9.05 01.03.07 we mp_primorial with mp_primor_cutoff
1.0.00 07.05.07 we mp_exptmod returns 0 for c=1, error for c<=0
1.0.01 08.05.07 we Removed Fp[sqrt(D)] arithmetic for MPC_sqrtmod_fp
1.0.02 09.05.07 we Removed mp_ecpwr2
1.0.03 13.05.07 we Corrected some exception strings
1.1.00 19.05.07 we break mp_gcd loop if u=v
1.1.01 22.05.07 we Binary Kronecker symbol algorithm
1.1.02 22.05.07 we Binary Kronecker: skip operations after mp_sub
1.1.03 23.05.07 we Binary Kronecker: s_mp_sub, break if x=y
1.1.04 31.05.07 we mp_is_power: check mod 4 before mp_expt_int
1.1.05 27.06.07 we mp_gcd_initial_mod in mp_gcd
1.1.06 01.07.07 we removed useless first iteration in mp_fact
(thanks to Marcel Martin)
1.1.07 02.07.07 we mp_miller_rabin uses IsPrime16 if possible
1.1.08 09.07.07 we mp_gcd_initial_mod renamed to mp_initial_mod
and also used in kron_intern
1.1.09 12.07.07 we mp_pell and mp_pell4
1.1.10 14.07.07 we s_mp_lucasvmod1 (used by mp_is_slpsp)
1.1.11 14.07.07 we mp_sqrtmod14_mueller
1.1.12 14.07.07 we Mueller sqrtmod: Bugfix use 1/t mod p
1.1.13 15.07.07 we Mueller sqrtmod: case t=1 handled in Jacobi loop
1.1.14 15.07.07 we Mueller sqrtmod: don't compute PP in Jacobi loop
1.1.15 21.07.07 we mp_xgcd: easy out for trivial cases
1.1.17 30.07.07 we mp_sqrtmodpk
1.1.17 31.07.07 we mp_sqrtmod14_mueller for mp_sqrtmethod=3 and
in auto mode instead of Lucas
1.1.18 01.08.07 we fast_mp_invmod: no init checks, y removed
1.1.19 04.08.07 we mp_sqrtmod2k
1.1.20 05.08.07 we special cases k=1,2 of mp_sqrtmod2k
easy out in mp_exptmod if b=1
1.1.21 05.08.07 we mp_sqrtmod_ex with optional Jacobi check
1.1.22 05.08.07 we mp_sqrtmodpq, removed mp_invmod_euclid
1.2.00 16.08.07 we mp_lcm: special cases, divide larger arg by gcd
1.2.01 17.08.07 we mp_is_primepower
1.2.02 21.08.07 we jacobi32
1.2.03 23.08.07 we mp_jacobi_lm
1.2.04 24.08.07 we mp_jacobi_lm used in: mp_is_slpsp, mp_sqrtmod18_shanks_intern
1.2.05 26.08.07 we improved error handling in mp_sqrtmod18_* routines
1.2.06 30.08.07 we t=1 in mp_sqrtmod14_mueller, improved jacobi32
1.2.07 02.09.07 we updated mp_primor_cutoff values
1.2.08 03.09.07 we mp_OddProd, mp_dfact, mp_fact with local mp stack
1.2.09 04.09.07 we small tables for mp_lucas and mp_fib
1.2.10 04.09.07 we Bugfix mp_fib2 for n=k*$10000000, k=1..7
1.2.11 06.09.07 we mp_binomial, s_mp_binom_w
1.2.12 07.09.07 we mp_binomial for negative arguments
1.2.13 08.09.07 we bugfix mp_binomial for small k and n>65535
1.2.14 08.09.07 we s_mp_binom_l, case k=1 in mp_binomial
1.2.15 09.09.07 we mp_binomial: init X[s] when needed
1.2.16 10.09.07 we mp_sqrtmod916_kong
1.2.17 10.09.07 we avoid warning in mp_miller_rabin if MP_16BIT
1.2.18 11.09.07 we mp_jacobi_ml
1.2.19 12.09.07 we Bugfix mp_is_primepower for a=2^1
1.2.20 17.09.07 we mp_is_power uses residues mod 8
1.3.00 26.10.07 we s_mp_sqrtmod2k for odd a, mp_sqrtmod2k for a>=0
1.3.01 12.11.07 we Fix memory leak(s) if MPC_HaltOnError is not defined
1.3.02 10.12.07 we improved mp_is_square2
1.5.00 13.01.08 we mp_sqrtmod2k: bugfix due to error in Adler/Coury [26]
1.6.00 24.05.08 we mp_crt_solve, mp_crt_setup, mp_crt_single
1.6.01 24.05.08 we functions mp_invmodf, mp_crt_setupf
1.6.02 04.06.08 we mp_is_square2: psqrt^ := a if a=0 or 1
1.6.03 04.06.08 we mp_cornacchia
1.6.04 05.06.08 we mp_cornacchia4
1.6.05 06.06.08 we function mp_gcd1
1.6.06 08.06.08 we Cosmetic changes / corrected exception strings
1.6.07 11.06.08 we Arg check in mp_binomial
1.7.00 27.07.08 we Improved trial factors in mp_isMersennePrime
1.7.01 17.09.08 we use mp_is_longint after mp_initial_mod in mp_gcd
1.7.02 17.09.08 we mp_gcd_euclid with mp_mod, more uses of mp_is_longint
1.7.03 20.09.08 we mp_xgcd with Lehmer, mp_gcd_ml (modified Lehmer)
1.7.04 21.09.08 we $ifdef MPC_UseGCD32 in mp_gcd
1.7.05 22.09.08 we mp_xgcd_bin
1.7.06 22.09.08 we mp_invmodf via xgcd, removed fast_mp_invmod
1.7.07 22.09.08 we extended range for Lehmer routine
1.7.08 23.09.08 we improved inner loop of (binary) mp_gcd
1.7.09 24.09.08 we use mp_sign for trivial cases in mp_xgcd_bin
1.7.10 24.09.08 we mp_is_pprime_ex skips 1-digit test, fixed mp_is_square2
1.7.11 24.09.08 we mp_invmodf with invmod32
1.7.12 02.10.08 we bugfix in mp_gcd $ifdef MPC_UseGCD32
1.8.00 09.10.08 we use mp_set1; mp_is_longint in mp_miller_rabin
1.8.01 21.10.08 we exptmod32, s_mp_is_pth_power, improved mp_is_power
1.8.02 22.10.08 we changed mp_is_primepower
1.8.03 25.10.08 we mp_is_power with s_mp_n_root2
1.8.04 26.10.08 we Check a mod k^3 and a mod q^2 in s_mp_is_pth_power
1.8.05 27.10.08 we improved mp_is_power
1.8.06 01.11.08 we mp_xlcm, mp_rnr
1.8.07 02.11.08 we sieve test loop in s_mp_is_pth_power, new: mp_is_pth_power
1.9.00 06.12.08 we mp_rnr2
1.9.01 23.12.08 we s_mp_mca_alg1816
1.9.02 25.12.08 we improved mp_small_factor
1.9.03 29.12.08 we uses mp_prime
1.9.04 03.01.09 we mp_is_square2, s_mp_pell4 with s_mp_sqrtrem
1.9.05 04.01.09 we mp_is_square2: s_mp_sqrtrem only (removed trick branch)
removed mp_ecm_brent, mp_ecpwr
1.9.06 04.01.09 we 32 bit first/next prime routines to mp_prime
1.9.07 06.01.09 we minor changes in mp_isMersennePrime
1.9.08 06.01.09 we explicit range check in mp_fermat
1.10.00 21.01.09 we changes related to (s)mp_divrem
1.10.01 22.01.09 we mp_cbrtmod
1.10.02 25.01.09 we s_mp_is_pprime_ex
1.10.03 28.01.09 we improved mp_is_square2
1.10.04 02.02.09 we mp_sqrtmodpk: check invmod, store result only if Err=0
1.10.05 03.02.09 we s_mp_sqrtmodpk with red parameter, mp_sqrtmodpk with red=true
1.10.06 03.02.09 we mp_cbrtmodpk, s_mp_cbrtmodpk, mp_cbrtmod3k
1.10.07 04.02.09 we sqrt/cbrtmodpk functions: handle a mod p = 0
1.10.08 04.02.09 we mp_sqrtmodpk handles p=2
1.10.09 07.02.09 we mp_cbrtmodpq
1.10.10 07.02.09 we s_mp_is_cubres
1.10.11 08.02.09 we mp_cbrtmod_ex
1.10.12 12.02.09 we improved mp_cbrtmod_ex
1.10.13 19.02.09 we improved s_mp_is_cubres
1.10.14 27.02.09 we mp_isMersennePrime: check for spurious factor
1.12.00 20.06.09 we mp_is_power_max
1.12.01 21.06.09 we mp_provable_prime
1.12.02 05.07.09 we new: s_mp_npfirst, s_mp_npnext, updated: mp_nextprime, s_mp_is_pprime_ex
1.12.04 29.07.09 we Increased trace level in mp_provable_prime
1.13.00 15.08.09 we mp_crt_single/solve: improved and check n>0
1.13.01 23.08.09 we allow p=q in mp_sqrtmodpq/mp_cbrtmodpq
1.16.00 04.06.10 we mp_4sq, mp_4sq_sa, mp_4sq_sd
1.17.00 30.12.10 we mp_binomial for k < 0
1.19.00 03.11.11 we Adjust MaxFermat, MaxFact for MAXDigits=32000 with MP_32BIT
1.20.00 14.01.12 we BIT64: exptmod32
1.20.01 21.01.12 we Adjusted MaxPrimorial, MaxPrimorial, NMAX in mp_primorial
1.21.00 15.07.12 we mp_isconsole: boolean [used for write('.')]
1.21.01 18.07.12 we mp_catalan, MaxCatalan
1.21.02 19.07.12 we improved mp_small_factor
1.21.03 20.07.12 we mp_poch
1.21.04 20.07.12 we mp_perm
1.21.05 21.07.12 we mp_poch and mp_perm use mp_product
1.21.06 25.07.12 we exptmod32/jacobi32 moved to mp_base
1.21.07 27.07.12 we mp_is_spsp returns false for n<=1
1.21.08 29.07.12 we mp_is_square2 uses is_square32ex
1.21.09 29.07.12 we mp_squfof
1.21.10 30.07.12 we mp_sigmak
1.21.11 31.07.12 we mp_sigmak with geometric sum to avoid overflow
1.22.00 10.08.12 we mp_nextprime_ex
1.22.01 11.08.12 we mp_small_factor with inline IsPrime16
1.22.02 12.08.12 we mp_pollard_pm1 : max. bound = 65000, fix while condition bug
mp_williams_pp1: max. bound = 65000
1.22.03 14.08.12 we Borwein/Schoenhage s_mp_borsch_fact
1.22.04 15.08.12 we mp_pollard_brent/ex
1.22.05 16.08.12 we Changed mp_lucasv2p; new: mp_lucasuv, mp_lucasu
1.22.06 18.08.12 we Removed mp_coshmult
1.22.07 19.08.12 we mp_squad_mod
1.22.08 21.08.12 we mp_cornacchia_ex
1.22.09 27.08.12 we mp_qnr
1.22.10 01.09.12 we special treatment for d=12 in s_mp_pell4
1.22.11 06.09.12 we mp_rqffu and s_mp_rqffu
1.22.12 08.09.12 we mp_qnr for all n > 0
1.22.13 11.09.12 we mp_powerd
1.22.14 12.09.12 we mp_ecm_factor uses Primes16
1.23.00 16.09.12 we fix s1 in mp_small_factor
1.23.01 24.09.12 we move s_mp_mca_alg1816 to mp_rsa
1.23.02 24.09.12 we s_mp_nextprime_sieve
1.23.03 24.09.12 we basic modular arithmetic, GCD/LCM moved to mp_modul
1.23.04 24.09.12 we factorization routines moved to mp_pfu
1.23.05 28.09.12 we mp_cornacchia_ex renamed to mp_cornacchia_pk
1.23.06 29.09.12 we s_mp_cornacchia_ex, mp_cornacchia_pq
1.23.07 30.09.12 we improved mp_is_power
1.23.08 02.10.12 we BIT: 16fix case in s_mp_is_pth_power with p > 2^16
1.23.09 03.10.12 we mp_val, mp_valrem, s_mp_valrem_d
1.23.10 03.10.12 we mp_is_power with mp_valrem
1.23.11 04.10.12 we Fix bug for test d=5,12 in s_mp_pell4
1.23.12 14.10.12 we mp_is_power with flexible trial division
1.24.00 17.12.12 we Some word types changed to longint
1.24.01 19.12.12 we mp_fact with longint arguments, s_mp_borsch_fact uses sieve
1.24.02 20.12.12 we MaxFact, MaxFermat computed in initialization
1.24.03 20.12.12 we mp_binomial with prime_sieve for 32+ bits
1.24.04 20.12.12 we mp_dfact and mp_OddProd with longint arguments
1.24.05 21.12.12 we s_mp_gbinom (used in mp_binomial and mp_catalan)
1.24.06 03.01.13 we Improved mp_primorial
1.24.07 04.01.13 we Improved s_mp_sqrtmod2k; mp_sqrtmod2k: No error if a is a multiple of 2^k
1.24.08 04.01.13 we Handle k=0 in s_mp_cbrtmodpk, error if k<0
1.27.00 04.09.13 we Allow a < 0 in mp_sqrtmod2k
1.27.01 08.02.14 we mp_val/rem: fix for a=0, changed iterative logic
1.29.00 15.07.14 we s_mp_nextprime_sieve can return safe primes
1.29.01 16.07.14 we mp_safeprime
1.29.02 16.07.14 we MPC_SieveNextPrime: use sieve for mp_next_prime
1.29.03 17.07.14 we s_mp_is_psp2
1.29.04 17.07.14 we s_mp_nextprime_sieve with s_mp_is_psp2
1.29.05 18.07.14 we mp_is_spsp with Barret reduction if s >= SBMin
1.29.06 18.07.14 we mp_rand_prime with mp_safeprime for pt_safe type
1.30.00 29.07.14 we Improved s_mp_sqrtmodpk for k a power of 2
1.30.01 11.09.14 we Improved s_mp_cbrtmodpk for k a power of 2
1.30.02 20.09.14 we mp_is_power_modpk
1.30.03 21.09.14 we s_mp_nroot_modp
1.30.04 21.09.14 we s_mp_nroot_modpk
1.30.05 22.09.14 we s_mp_nroot_modpq
1.30.06 25.09.14 we Fixed mp_is_power_modpk with n: mp_int
1.30.07 28.09.14 we mp_safeprime_ex
1.30.08 29.09.14 we Improved mp_perm
1.30.09 29.09.14 we Check/handle n<0 in mp_fact
1.31.00 17.10.14 we small change in s_mp_is_pth_power
1.31.01 05.11.14 we Fix/extend mp_dfact for n<0
1.31.02 22.11.14 we Extended s_mp_nroot_modp
1.32.00 03.03.15 we Exception strings corrected in mp_cbrtmod3k, mp_cbrtmod_ex
1.33.00 14.09.16 we Allow n > p-1 in s_mp_nroot_modp
1.34.00 09.12.16 we mp_is_pcarmichael/ex
1.35.00 25.07.17 we mp_fact: check n > MaxFact
1.35.01 26.07.17 we mp_fact_swing
1.35.02 26.07.17 we mp_swing
1.36.00 10.09.17 we mp_is_psp
1.36.01 11.09.17 we mp_is_psp_d
1.36.02 12.09.17 we mp_nextpower/_ex
1.36.03 18.09.17 we mp_digitsum
1.36.04 21.09.17 we s_mp_is_lpsp with Selfridge A parameters
1.36.05 22.09.17 we Rename existing strong Lucas test to mp_is_slpsp_alt
1.36.06 22.09.17 we mp_is_lpsp, mp_is_slpsp with s_mp_is_lpsp
1.37.00 15.10.17 we mp_reverse
**************************************************************************)
(*-------------------------------------------------------------------------
This code uses material/ideas from the following 3rd party libraries:
- LibTomMath 0.30+ by Tom St Denis
- MPI 1.8.6 by Michael J. Fromberger
- NX V0.18 and V0.9+ by Marcel Martin
See the file '3rdparty.mpa' for the licenses.
----------------------------------------------------------------------------*)
(*-------------------------------------------------------------------------
(C) Copyright 2004-2017 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.
----------------------------------------------------------------------------*)
const
mp_max_small = mp_digit(MP_DIGIT_MAX and $7FFF); {max. tested small factor}
var
MaxFact: longint; {MaxFact depends on DIGIT_BIT and MAXDigits. It }
{is an approx. value for max. valid mp_fact args.}
{Since no checks are performed some larger values}
{might work. With 16 bit real mode even smaller }
{values may result in "Out of memory" errors. }
var
MaxFermat: integer; {Maximum Fermat number with this implementation }
var
mp_sqrtmethod : integer; {Selects algorithm for sqrt mod p, p=1 mod 8.}
{0:Auto, 1:Shanks, 2:Lucas, 3:Mueller.}
{Default=0: auto select between Shanks/Mueller.}
var
mp_primor_cutoff: longint; {Cutoff point for recursive product in mp_primorial,}
{initialization will set rough value. For better}
{results use t_tunepr}
type
TPrimeType = (pt_normal, pt_3mod4, pt_safe, pt_1mod8);
{pt_normal: prime without special features.}
{ pt_3mod4: prime p with p modulo 4 = 3. }
{ pt_1mod8: prime p with p modulo 8 = 1. }
{ pt_safe: safe prime, ie p=2q+1, q prime.}
const
MaxMersenne = DIGIT_BIT*MAXDigits-1; {Approx. maximum Mersenne number}
MaxPrimorial = trunc(MaxMersenne*0.6932); {Approx. maximum arg for mp_primorial, Note: for real}
{mode "out of memory" will occur for smaller values!}
MaxFibonacci = trunc(MaxMersenne*1.4404); {Approx. maximum arg for mp_fib}
MaxLucas = MaxFibonacci-2; {Approx. maximum arg for mp_lucas}
MaxCatalan = trunc(MaxMersenne*0.4999); {Approx. maximum arg for mp_catalan}
procedure mp_4sq(const m: mp_int; var a,b,c,d: mp_int);
{-Decompose m>=0 into 4 squares: m = a^2 + b^2 + c^2 + d^2.}
{ If no solutions for m>0 was found, a=b=c=d=0 is returned.}
procedure mp_4sq_sa(const m: mp_int; var a,b,c,d: mp_int);
{-Decompose m>=0 into 4 squares: m = a^2 + b^2 + c^2 + d^2, a<=b<=c<=d.}
{ If no solutions for m>0 was found, a=b=c=d=0 is returned.}
procedure mp_4sq_sd(const m: mp_int; var a,b,c,d: mp_int);
{-Decompose m>=0 into 4 squares: m = a^2 + b^2 + c^2 + d^2, a>=b>=c>=d.}
{ If no solutions for m>0 was found, a=b=c=d=0 is returned.}
procedure mp_binomial(n,k: longint; var a: mp_int);
{-Calculate the binomial coefficient a = (n choose k)}
procedure mp_catalan(n: longint; var c: mp_int);
{-Return the n'th Catalan number C_n = binomial(2n,n)/(n+1), 0 if n < 0}
procedure mp_cbrtmod(const a, p: mp_int; var x: mp_int; var Err: integer);
{-Compute a cube root x of a with x^3 = a mod p, p prime. }
{ Err in [1..3], if a is not a 3rd power or p is not prime.}
procedure mp_cbrtmod3k(const a: mp_int; k: longint; var b: mp_int; var Err: integer);
{-Compute a cube root b of a with b^3 = a mod 3^k. Error codes: 1<Err<4 code from}
{ mp_cbrtmod, Err=4 if p=3, Err=5: no inverse mod 3^(2^i), Err=6: (internal d<>0)}
procedure mp_cbrtmodpk(const a,p: mp_int; k: longint; var b: mp_int; var Err: integer);
{-Compute a cube root b of a with b^3 = a mod p^k, p prime. b will be }
{ reduced mod p^k. Error codes: 1<Err<4: from mp_cbrtmod, Err=4 if p=3,}
{ Err=5: no inverse mod p^(2^i), Err=6: internal check for p=3 failed. }
procedure mp_cbrtmodpq(const a,p,q: mp_int; var x: mp_int; var Err: integer);
{-Compute a cube root x of a mod (pq); p,q primes. If p=q, mp_cbrtmodpk is}
{ used. For p<>q: Err=8 if gcd(p,q)<>1, otherwise error code from mp_cbrtmod.}
procedure mp_cbrtmod_ex(const a, p: mp_int; var x: mp_int; pru: pmp_int; var Err: integer);
{-Compute a cube root x of a with x^3 = a mod p, p prime. If pru<>nil, pru^ }
{ is set to a 3rd root of unity (or 1 if x is unique), the other roots are }
{ x*pru^, x*pru^*pru^. Err <> 0, if a is not a 3rd power or p is not prime. }
procedure mp_cornacchia(const d,p: mp_int; var x,y: mp_int; var Err: integer);
{-Solve x^2 + d*y^2 = p, p prime, 0<d<p with Cornacchia's algorithm. Check}
{ @x<>@y, but not p prime. Err<>0 if no solution exist. x >= y if d=1.}
procedure mp_cornacchia_pk(const a,b,p: mp_int; k: longint; var x,y: mp_int; var Err: integer);
{-Solve a*x^2 + b*y^2 = p^k, p prime; k,a,b > 0, a + b < p^k, (a,p)=1 with Cornacchia's }
{ algorithm. Check @x<>@y, but not p prime. Err <> 0 if no solution exist. x >= y if a=b.}
procedure mp_cornacchia_pq(const a,b,p,q: mp_int; var x,y: mp_int; var Err: integer);
{-Solve a*x^2 + b*y^2 = p*q, p,q prime; a,b > 0, a + b < p*q, (a,p*q)=1 with Cornacchia's }
{ algorithm. Check @x<>@y, but not p,q prime. Err <> 0 if no solution exist. x >= y if a=b.}
procedure mp_cornacchia4(const d,p: mp_int; var x,y: mp_int; var Err: integer);
{-Solve x^2 + |d|*y^2 = 4p, p prime, -4p<d<0 with the modified Cornacchia}
{ algorithm. Check @x<>@y, but not p prime. Err<>0 if no solution exist.}
procedure mp_crt_setup(n: integer; const m: array of mp_int; var c: array of mp_int);
{-Calculate CRT coefficients c[i] for pairwise co-prime moduli m[i], i=0..n-1.}
function mp_crt_setupf(n: integer; const m: array of mp_int; var c: array of mp_int): boolean;
{-Calculate CRT coefficients c[i] for pairwise co-prime moduli m[i], i=0..n-1.}
{ Return true if c[i] are successfully calculated.}
procedure mp_crt_single(n: integer; const m,v: array of mp_int; var x: mp_int);
{-Calculate x with x mod m[i] = v[i], pairwise co-prime moduli m[i], i=0..n-1.}
{ Use mp_crt_setup/calc if more than one system shall be solved}
procedure mp_crt_solve(n: integer; const m,c,v: array of mp_int; var x: mp_int);
{-Calculate x with x mod m[i] = v[i], pairwise co-prime moduli m[i], i=0..n-1.}
{ Coefficients c[i] must be precalculated with mp_crt_setup}
function mp_digitsum(const a: mp_int; radix: word): longint;
{-Return the digitsum of a in base radix (2..MAXRadix)}
procedure mp_dfact(n: longint; var a: mp_int);
{-Calculate double factorial a = n!!, a=n*(n-2)..., lowest term 1 or 2; a=0 for n < -3}
procedure mp_fact(n: longint; var a: mp_int);
{-Calculate a = factorial(n) using Recursive Split or Borwein/Schoenhage}
{ prime factorization method, error if n > MaxFact}
procedure mp_fact_swing(n: longint; var a: mp_int);
{-Calculate a = factorial(n) using Luschny's prime swing method}
procedure mp_fermat(n: word; var fn: mp_int);
{-Return nth Fermat number, fn = 2^(2^n)+1 (MP_RANGE error for n>MaxFermat)}
procedure mp_fib(n: longint; var fn: mp_int);
{-Calculate Fibonacci number fn=fib(n), fib(-n)=(-1)^(n-1)*fib(n)}
procedure mp_fib2(n: longint; var fn,f1: mp_int);
{-Calculate two Fibonacci numbers fn=fib(n), f1=fib(n-1), n>=0}
function mp_isMersennePrime(p: longint): boolean;
{-Lucas-Lehmer test for Mersenne number m=2^p-1, HAC 4.37}
function mp_is_lpsp(const a: mp_int): boolean;
{-Lucas probable prime test for n using Selfridge A method}
function mp_is_pcarmichael(const a: mp_int): boolean;
{-Test if a is is a Carmichael number with high probability, 100 tests, prime check}
function mp_is_pcarmichael_ex(const a: mp_int; numtests: integer; chkprime: boolean): boolean;
{-Test if a is a Carmichael number with high probability, perform numtests}
{ tests (100 if numtests < 1) and optionally checks if a is prime.}
procedure mp_is_power(const a: mp_int; var b: mp_int; var p: longint);
{-Calculate smallest prime p with a=b^p; p=1,b=a if a is no power}
procedure mp_is_power_max(const a: mp_int; var b: mp_int; var k: longint);
{-Calculate largest k with a=b^k; k=1,b=a if a is no power}
function mp_is_power_modpk(a,n,p: mp_int; k: longint): boolean;
{-Return true if a is a nth power residue mod p^k;}
{ p prime (not checked); gcd(a,p)=1; n, k > 0. }
function mp_is_pprime(const a: mp_int): boolean;
{-Test if a is prime (BPSW probable prime if a>2^32)}
function mp_is_pprime_ex(const a: mp_int; smax: mp_digit): boolean;
{-Test if a is prime (BPSW probable prime if a>2^32); trial division up to smax}
function mp_is_primepower(const a: mp_int; var b: mp_int; var k: longint): boolean;
{-Return true if a=b^k, b prime, k>1, otherwise false and a=b^k, k=1 if no power}
function mp_is_pth_power(const a: mp_int; p: longint; var r: mp_int): boolean;
{-Return true if a is pth power, a>0, p prime. If true, calculate r with a=r^p}
function mp_is_psp(const n,a: mp_int): boolean;
{-Probable prime test of n to base a > 1, true if a^(n-1) mod n = 1}
function mp_is_psp_d(const n: mp_int; a: mp_digit): boolean;
{-Probable prime test of n to base a > 1, true if a^(n-1) mod n = 1}
function mp_is_slpsp_alt(const a: mp_int): boolean;
{-Strong Lucas probable prime test for a. Lucas test is }
{ done for the first p=2k+1 with mp_jacobi(p^2-4,a) = -1}
function mp_is_slpsp(const a: mp_int): boolean;
{-Lucas (strong) probable prime test for n using Selfridge A method}
function mp_is_spsp(const n,a: mp_int): boolean;
{-Strong probable prime test of n to base a > 1 from HAC p. 139 Alg.4.24}
function mp_is_spsp_d(const n: mp_int; a: mp_digit): boolean;
{-Strong probable prime test of n to mp_digit base a > 1 from HAC p. 139 Alg.4.24}
function mp_is_square(const a: mp_int): boolean;
{-Test if a is square}
function mp_is_square2(const a: mp_int; psqrt: pmp_int): boolean;
{-Test if a is square, return sqrt(a) if a is a square and psqrt<>nil}
procedure mp_lucas(k: longint; var lk: mp_int);
{-Calculate Lucas number lk=luc(k), luc(-k)=(-1)^k*luc(k)}
procedure mp_lucas2(k: longint; var lk,l1: mp_int);
{-Calculate two Lucas numbers lk=luc(k), l1=luc(k-1), k>=0}
procedure mp_lucasv(const p,q: mp_int; k: longint; var v: mp_int);
{-Calculate v[k] of Lucas V sequence for p,q, p^2-4q <>0, k>=0}
procedure mp_lucasuv(const p,q: mp_int; k: longint; var uk,vk: mp_int);
{-Calculate u[k], v[k] of Lucas sequence for p,q; p^2-4q<>0, k>=0}
procedure mp_lucasu(const p,q: mp_int; k: longint; var u: mp_int);
{-Calculate u[k] of Lucas sequence for p,q, p^2-4q <>0, k>=0}
procedure mp_lucasvmod(const p,q,n,k: mp_int; var vk: mp_int);
{-Calculate v[k] mod n of Lucas V sequence for p,q.}
{ Ranges n>1, 0 <= p,q,k < n (checked if MPC_ArgCheck).}
{ Note: p^2-4q<>0 is not tested, no proper Lucas sequence!!}
procedure mp_lucasv2(const p,q: mp_int; k: longint; var v,w: mp_int);
{-Calculate v=v[k],w=v[k+1] of Lucas V sequence for p,q, p^2-4q<>0, k>=0}
procedure mp_lucasv2p(const p,q: mp_int; k: longint; var vk: mp_int; p2: pmp_int; uk: boolean);
{-Calculate v[k] of Lucas V sequence for p,q; p^2-4q<>0, k>=0;}
{ if p2<>nil, p2^ will be set to v[k+1] or u[k] if uk is true.}
procedure mp_mersenne(n: longint; var mn: mp_int);
{-Return nth Mersenne number, mn = 2^n-1, MP_RANGE err for n>MaxMersenne}
procedure mp_miller_rabin(const n: mp_int; t: word; var prime: boolean);
{-Miller-Rabin test of n, security parameter t, from HAC p. 139 Alg.4.24}
{ if t<=0, calculate t from bit_count with prob of failure < 2^-96}
procedure mp_nextpower(var a: mp_int);
{-Return smallest perfect power >= a}
procedure mp_nextpower_ex(const a: mp_int; var b,c: mp_int; var p: longint);
{-Return smallest perfect power b := c^p >= a}
procedure mp_nextprime(var n: mp_int);
{-Next prime >= n, 2 if n<=2}
procedure mp_nextprime_ex(var n: mp_int; smax: mp_digit);
{-Next prime >= n, 2 if n<=2; trial division from 3 to smax}
procedure mp_OddProd(a,b: longint; var p: mp_int);
{-Calculate p=prod(2*i+1),i=a+1...b; p=1 if a>=b}
procedure mp_pell(const d: mp_int; var x,y: mp_int);
{-Calculate the smallest non-trivial solution of the Pell equation}
{ x^2 - d*y^2 = 1; error if d<2, if d is a square, or if @x = @y.}
procedure mp_pell4(const d: mp_int; var x,y: mp_int; var r: integer);
{-Calculate the smallest non-trivial solution of x^2 - d*y^2 = r,}
{ r in [-4,+4,-1,+1]. Uses continued fraction expansion of sqrt(d)}
{ from Forster [4]. Error if d<2, if d is a square, or if @x = @y.}
procedure mp_perm(n,r: longint; var a: mp_int);
{-Compute a=n!/(n-r)!, the number of permutations of n distinct objects}
{ taken r at a time, n,r >= 0.}
procedure mp_poch(n,k: longint; var a: mp_int);
{-Return the Pochhammer symbol a = (n)_k = n*(n+1)*...(n+k-1), a=1 if k<1}
{ (n)_k is sometimes called "rising factorial" or "Pochhammer function". }
procedure mp_powerd(const p,q,d: mp_int; n: longint; var x,y: mp_int);
{-Calculate x + y*sqrt(d) = (p + q*sqrt(d))^n, n >= 0}
procedure mp_prevprime(var n: mp_int);
{-Previous prime <= n, 0 if n<2}
procedure mp_primorial(n: longint; var a: mp_int);
{-Primorial of n; a = n# = product of primes <= n.}
procedure mp_provable_prime(bits: longint; var p: mp_int);
{-Generate a random provable prime p with bitsize bits using Maurer's algorithm}
function mp_qnr(const n: mp_int): longint;
{-Return a small quadratic nonresidue for n > 0, -1 if error or no QNR is}
{ found: the smallest prime with kronecker(p|n) = -1 is returned. If n is}
{ a prime, the result is the least positive QNR, for composite n this may}
{ differ from the least QNR: e.g. mp_qnr(15) returns 7 and not 2.}
procedure mp_rand_prime(bitsize: longint; pt: TPrimeType; var p: mp_int);
{-Generate random (probable BPSW) prime of bitsize > 3, pt: prime type of p}
procedure mp_reverse(const a: mp_int; radix: mp_digit; var b: mp_int);
{-Calculate b with the reversed digits of a in base radix (2..MAXRadix)}
procedure mp_rnr(const a,m: mp_int; var x,y: mp_int);
{-Rational number reconstruction: for m>0 calculate x,y with a*x=y mod m,}
{ gcd(x,y)=1, 0<=x,|y|<sqrt(m/2), x<>0. x=y=0 if no (unique) solution exists.}
procedure mp_rnr2(const a,m,NN,DD: mp_int; var n,d: mp_int);
{-Rational number reconstruction: for m,NN,DD > 0 calculate co-prime d,n}
{ with a*d=n mod m, |n|<=N, 0<d<=DD, i.e. a=n/d mod m. Return d=n=0 if }
{ no solution exists. The reconstruction is unique if 2*NN*DD < m.}
function mp_rqffu(d: longint; var u,v: mp_int): integer;
{-Return the fundamental unit e = (u + v*sqrt(d))/2 of the real quadratic field}
{ with discriminant d>4, d mod 4 = 0,1. Result is the norm +/- 1 or 0 if error.}
procedure mp_safeprime(var n: mp_int);
{-Compute the next safe prime >= n: new n and (n-1)/2 are prime, 5 if n<=5}
procedure mp_safeprime_ex(var p,g: mp_int; smallest: boolean);
{-Compute the next safe prime >= p and a generator g of Z_p^*, smallest g or random}
procedure mp_sigmak(k,n: longint; var a: mp_int);
{-Compute the sum of the kth powers of the divisors of n; zero if k<0 or n=0.}
{ Special cases k=0: number of divisors of n; k=1: sum of the divisors of n.}
procedure mp_small_factor(const a: mp_int; f0,fmax: mp_digit; var f: mp_digit);
{-Compute small digit prime factor or 0, f0..fmax, f will be <= min(fmax,$7FFF)}
procedure mp_sqrtmod(const a,p: mp_int; var b: mp_int; var Err: integer);
{-Calculate square root b of a with b*b = a mod p, p prime, with Jacobi check.}
{ Err=-1 if (a|p)<>1, Err=1 if failure for p=1 mod 8, Err=2 if p even and <>2}
procedure mp_sqrtmod_ex(const a,p: mp_int; ChkJ: boolean; var b: mp_int; var Err: integer);
{-Calculate square root b of a with b*b = a mod p, p prime.}
{ If ChkJ=true then check Jacobi symbol (a|p)=1, Err=-1 if check fails.}
{ Err=1 if failure for p=1 mod 8, Err=2 if p even and not 2.}
procedure mp_sqrtmod2k(const a: mp_int; k: longint; var b: mp_int; var Err: integer);
{-Calculate a square root b of an integer a with b*b = a mod 2^k.}
{ Return Err=-1 if there is no solution. For odd a the solutions }
{ are normalized to 0 < b < 2^k/4 for k>3.}
procedure mp_sqrtmodp2(const a,p: mp_int; var b: mp_int; var Err: integer);
{-Calculate square root b of a with b*b = a mod p^2, p prime.}
{ Alias for mp_sqrtmodpk(,,2,,). Err: error code from mp_sqrtmodpk.}
procedure mp_sqrtmodpk(const a,p: mp_int; k: longint; var b: mp_int; var Err: integer);
{-Calculate square root b < p^k of a with b*b = a mod p^k, p prime. }
{ Error codes: if p=2: Err=1 if a<0; Err=-1 if there is no solution }
{ if p<>2: Err=-1 if (a|p)<>1, Err=1: failure for p=1 mod 8, Err=2 if}
{ p is even, Err=4: no inverse mod p^(2^i)}
procedure mp_sqrtmodpq(const a,p,q: mp_int; var x,y: mp_int; var Err: integer);
{-Calculate square roots +x,-x,+y,-y of a mod (pq); p,q primes}
{ If p=q, x is the root from mp_sqrtmodp2 and y is not computed.}
{ For p<>q: Err=4 if gcd(p,q)<>1, otherwise error code from mp_sqrtmod}
function mp_squad_mod(const a,b,c,p: mp_int; var x1,x2: mp_int): integer;
{-Solve a^x + b*x + c = 0 mod p, p odd prime; returns number of roots: 0,1,2.}
{ Primality of p is not checked, if 2a has no inverse, b*x + c = 0 is solved.}
procedure mp_swing(n: longint; var a: mp_int);
{-Calculate Luschny swing number a = swing(n) = n!/((n div 2)!)^2}
function mp_val(const a: mp_int; r: longint): longint;
{-Return the valuation of a with respect to r, ie. the largest v with}
{ r^v divides a, v=0 if r=0,1,or -1, v=MaxLongint if a=0}
procedure mp_valrem(const a: mp_int; r: longint; var b: mp_int; var v: longint);
{-Return valuation v of a with respect to r and remainder b: a = b*r^v}
{#Z+}
{---------------------------------------------------------------------------}
{- 'Internal' functions, don't use them unless you know what you are doing -}
{---------------------------------------------------------------------------}
{#Z-}
function mp_sqrtmod14_mueller(const a,p: mp_int; var b: mp_int): boolean;
{-Calculate b with b*b = a mod p, p prime = 1 mod 4, return success status.}
{ *internal* assumes a>1,p>1, no argument checks}
function mp_sqrtmod18_shanks_intern(const a,p: mp_int; var b,q: mp_int; k: longint): boolean;
{-Calculate b with b*b = a mod p, p prime = 1 mod 8, return success status.}
{ *internal* q and k must be setup with p-1 = 2^k*q}
function mp_sqrtmod18_shanks(const a,p: mp_int; var b: mp_int): boolean;
{-Calculate b with b*b = a mod p, p prime = 1 mod 8, return success status; *internal*}
function mp_sqrtmod18_lucas(const a,p: mp_int; var b: mp_int): boolean;
{-Calculate b with b*b = a mod p, p prime = 1 mod 8, return success status}
{ *internal* assumes a,p>1, no arg check}
procedure mp_sqrtmod34(const g,p: mp_int; var z: mp_int);
{-Calculate z with z*z = g mod p, p prime > 0, p mod 4 = 3; *internal*}
procedure mp_sqrtmod58(const g,p: mp_int; var z: mp_int);
{-Calculate z with z*z = g mod p, p prime = 5 mod 8; *internal*}
function mp_sqrtmod916_kong(const a,p: mp_int; var b: mp_int): boolean;
{-Calculate b with b*b = a mod p, p prime = 9 mod 16, return success status}
{ *internal* assumes a,p>1, no arg check}
procedure s_mp_binom_l(n: longint; k: word; var a: mp_int);
{-Internal binomial coefficient for small k, no init check}
procedure s_mp_gbinom(n,m,k: longint; var a: mp_int);
{-Calculate the generalized binomial coefficient a = n!/m!/k!, using prime}
{ power decomposition of a. a must an integer and 0 <= m <= n, 0 <= k <= n}
procedure s_mp_recsplit_fact(n: word; var a: mp_int);
{-Calculate a = factorial(n) using Recursive Split, error if n > MaxFact}
procedure s_mp_borsch_fact(n: longint; var a: mp_int);
{-Calculate a = factorial(n) with Borwein/Schoenhage prime factorization method}
procedure s_mp_cbrtmodpk(const a,p: mp_int; k: longint; red: boolean; var b: mp_int; var Err: integer);
{-Compute a cube root b of a with b^3 = a mod p^k, p prime <> 3, k >= 0.}
{ If red=true, reduce b mod p^k, otherwise b may be >= p^k. Error codes:}
{ 1<Err<4: from mp_cbrtmod, Err=4 if p=3, Err=5: no inverse mod p^(2^i) }
procedure s_mp_cornacchia_ex(const a,b,p,q: mp_int; k: longint; var x,y: mp_int; var Err: integer);
{-Solve a*x^2 + b*y^2 = m, m=p^k if k>0 otherwise m=p*q, a,b > 0, a+b < m,}
{ (a,m)=1 with Cornacchia's algorithm. Check @x<>@y, but not p,q prime. }
{ Err <> 0 if error or no solution exist; normalisation x >= y if a=b. }
function s_mp_is_cubres(const a, p: mp_int): boolean;
{-Simple test if a is a cubic residue mod p, p prime. Primality of p}
{ is not checked, but some trivial non-prime cases are detected.}
function s_mp_is_lpsp(const n: mp_int; strong: boolean): boolean;
{-Lucas (strong) probable prime test for n using Selfridge A method}
function s_mp_is_pprime_ex(const a: mp_int; smin,smax: mp_digit): boolean;
{-Test if a is prime (BPSW probable prime if a>2^32); trial division}
{ from smin to smax; no init check, no check if a is less than 2<31 }
function s_mp_is_pth_power(const a: mp_int; p: longint; var r: mp_int): boolean;
{-Return true if a is pth power, then a=r^p. a>0, p>2 prime, no checks}
function s_mp_is_psp2(const n: mp_int): boolean;
{-Test if n>2 is a base-2 (Fermat) probable prime, no init check}
procedure s_mp_lucasvmod1(const p,n,k,mu: mp_int; var v: mp_int);
{-Calculate v=v[k] mod n of Lucas V sequence for p,q=1. mu: Barrett parameter}
{ for n, k>0. Internal use: no init check, assumes p<n, return v=2 for k<=0.}
procedure s_mp_npfirst(var a: mp_int; var idx: integer);
{-Setup idx and increment a to first prime candidate, no init check, a>7}
procedure s_mp_npnext(var a: mp_int; var idx: integer);
{-Update idx and increment a to next prime candidate, no init check, a>7}
procedure s_mp_nextprime_sieve(safe: boolean; var n: mp_int);
{-Compute the next prime >= n using a segmented sieve, safe prime if safe=true}
procedure s_mp_nroot_modp(const a,n,p: mp_int; var x: mp_int; var Err: integer);
{-Solve x^n = a mod p; p prime (not checked), gcd(n, p-1)=1 or gcd(n, (p-1)/n)=1}
procedure s_mp_nroot_modpq(const a,n,p,q: mp_int; var x: mp_int; var Err: integer);
{-Solve x^n = a mod pq; p,q primes (not checked), gcd(n, phi(pq))=1; gcd(a,p)=1 if p=q}
procedure s_mp_nroot_modpk(const a,n,p: mp_int; k: longint; var x: mp_int; var Err: integer);
{-Solve x^n = a mod p^k; p prime (not checked), k>0, gcd(a,p)=1 and 1/n mod p^k exists}
procedure s_mp_pell4(const d: mp_int; var x,y: mp_int; var r: integer);
{-Calculate the smallest non-trivial solution of x^2 - d*y^2 = r}
{ r in [-4,+4,-1,+1]; returns r=0 if d<2 or d is a square. Uses }
{ continued fraction expansion of sqrt(d) from Forster [4], ch.25}
{ 'Internal' procedure: Assumes d,x,y are initialized and @x<>@y.}
function s_mp_rqffu(const d: mp_int; var x,y: mp_int): integer;
{-Return the fundamental unit e = (u + v*sqrt(d))/2 of the real quadratic field}
{ with discriminant d>4, d mod 4 = 0,1. Result is the norm +/- 1 or 0 if error.}
{ Note: This function uses s_mp_pell4 and is much slower than mp_rqffu but the }
{ discriminant may be large. Beware of very long running times for large d!! }
procedure s_mp_sqrtmod2k(const a: mp_int; k: longint; var b: mp_int; var Err: integer);
{-Calculate unique square root b of an odd integer a with b*b = a mod 2^k and }
{ 0<b<2^k/4 for k<3. Err=1 if a<0; =2 for even a; =3 if a<>1 mod min(2^k,8)}
procedure s_mp_sqrtmodpk(const a,p: mp_int; k: longint; red: boolean; var b: mp_int; var Err: integer);
{-Calculate square root b of a with b*b = a mod p^k, p odd prime.}
{ Err=-1 if (a|p)<>1, Err=1 if failure for p=1 mod 8, Err=2 if p even.}
{ Err=4 if no inverse mod p^(2^i). If red=true, reduce b mod p^k, otherwise}
{ b may be >= p^k; ex: a=22,p=3,k=3 --> b=34 and 34^2 = 22 mod 27}
procedure s_mp_valrem_d(const a: mp_int; d: mp_digit; var b: mp_int; var v: longint);
{-Return valuation v of a with respect to d and remainder b: a = b*d^v}
{ internal iterative version, d no power of two; v=MaxLongint, b=0 if a=0}
implementation
uses
mp_base, mp_modul, mp_prime;
{---------------------------------------------------------------------------}
procedure s_mp_npfirst(var a: mp_int; var idx: integer);
{-Setup idx and increment a to first prime candidate, no init check, a>7}
var
k: integer;
m: mp_digit;
begin
idx := -1;
{make n odd}
if mp_iseven(a) then mp_inc(a);
mp_mod_d(a, NPRC_MOD, m);
if mp_error<>MP_OKAY then exit;
m := m shr 1;
{move a to next prime residue class mod MP_MOD and index idx into diff array}
for k:=m to (NPRC_MOD div 2)-1 do begin
idx := NPRC_OddIdx[k];
{note: loop terminates via break because NPRC_OddIdx[(NPRC_MOD div 2)-1]<>-1}
if idx<>-1 then break;
mp_add_d(a,2,a);
end;
end;
{---------------------------------------------------------------------------}
procedure s_mp_npnext(var a: mp_int; var idx: integer);
{-Update idx and increment a to next prime candidate, no init check, a>7}
begin
if mp_error<>MP_OKAY then exit;
{move to next candidate}
mp_add_d(a,NPRC_Diff[idx],a);
{get next increment index}
inc(idx);
if idx>=NPRC_NRC then idx:=0;
end;
{---------------------------------------------------------------------------}
procedure mp_nextprime(var n: mp_int);
{-Next prime >= n, 2 if n<=2}
begin
{$ifdef MPC_SieveNextPrime}
{use sieve for mp_next_prime}
s_mp_nextprime_sieve(false, n);