This repository was archived by the owner on Aug 5, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathmp_base.pas
11198 lines (9914 loc) · 306 KB
/
mp_base.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_base;
{Multi precision integer arithmetic basic routines}
interface
{$ifdef VirtualPascal}
{$X+} {needed for pchars/RESULT}
{$endif}
{$i STD.INC}
{$ifdef BIT16}
{$N+}
{$X+} {needed for pchars}
{$endif}
uses
BTypes, mp_types;
{$i mp_conf.inc}
(*************************************************************************
DESCRIPTION : Multi precision integer arithmetic basic routines
REQUIREMENTS : BP7, D1-D7/D9-D10/D12/D17-D18, FPC, VP
EXTERNAL DATA : (mp_types)
MEMORY USAGE : lots of heap
DISPLAY MODE : ---
REFERENCES : [1] LibTomMath V0.30+ by Tom St Denis
[2] MPI by M.J. Fromberger
[3] D.E. Knuth, The Art of computer programming:
Volume 1, Fundamental Algorithms, 3rd ed., 1997;
Volume 2, Seminumerical Algorithms, 3rd ed., 1998;
http://www-cs-faculty.stanford.edu/~knuth/taocp.html
[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
[8] Marcel Martin: NX - Numerics library of multiprecision
numbers for Delphi and Free Pascal, 2006-2009
www.ellipsa.eu/public/nx/index.html
[10] Crandall,R., C.Pomerance: Prime Numbers, A Computational
Perspective, 2nd ed., 2005
[15] The GNU Multiple Precision Arithmetic Library, http://gmplib.org/
[24] H. Cohen, A Course in Computational Algebraic Number Theory
4th printing, 2000
[29] V. Shoup, A Computational Introduction to Number Theory and
Algebra, Version 2, 2008, from http://shoup.net/ntb/
[30] J. v. zur Gathen, J. Gerhard, Modern computer algebra, 2nd ed., 2003
http://math-www.uni-paderborn.de/mca/
[33] C. Burnikel, J. Ziegler: Fast Recursive Division. MPI fr Informatik,
Forschungsbericht MPI-I-98-1-022 (1998); available via
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.47.565
[34] P. Zimmermann, Karatsuba Square Root, INRIA Research Report RR-3805;
available from http://hal.inria.fr/inria-00072854/en/
[35] R.P. Brent, P. Zimmermann: Modern Computer Arithmetic, Cambridge University Press, 2010.
A preliminary version (V0.5.9, Oct. 2010) of the book is available from
http://maths-people.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf
or http://arxiv.org/abs/1004.4710 (V0.5.1)
[36] M. Bodrato, A.Zanoni, What About Toom-Cook Matrices Optimality?
available from http://bodrato.it/papers/#CIVV2006
See also M. Bodrato's page http://bodrato.it/software/
Version Date Author Modification
------- -------- ------- ------------------------------------------
0.0.01 09.05.04 W.Ehrhardt Initial version: BP7
0.0.02 10.05.04 we mp_copyprim
0.0.03 11.05.04 we mp_expand, mp_copy
0.0.04 12.05.04 we mp_newXY, mp_delete, LSign
0.0.05 22.07.04 we new design a la MPI/LibTomMath
mp_init, mp_clear, mp_exch, mp_zero
0.0.06 23.07.04 we mp_shrink, mp_grow, mp_init_size, mp_set,
mp_init_set, mp_clamp, mp_lshd
0.0.07 23.07.04 we cleanup/reordering, mp_iseven/odd/zero
0.0.08 23.07.04 we mp_copy
0.0.09 23.07.04 we debug/magic, mp_reverse, mp_toradix_n
0.0.10 24.07.04 we mp_div_d (only a=0,b=1 quick outs)
0.0.11 24.07.04 we mp_div_2d
0.0.12 24.07.04 we mp_rshd
0.0.13 24.07.04 we mp_mod_2d
0.0.14 24.07.04 we mp_mul_d
0.0.15 24.07.04 we mp_set_int
0.0.16 24.07.04 we mp_mul_2d
0.0.17 24.07.04 we fixed severe bug in mp_lshd
0.0.18 25.07.04 we debug code checking initialized mp_ints
0.0.19 25.07.04 we mp_read_radix
0.0.20 25.07.04 we mp_add_d
0.0.21 26.07.04 we mp_sub_d
0.0.22 26.07.04 we mp_div_2, mp_mul_2
0.1.00 26.07.04 we const parameters, single digit functions
complete, code clean up
0.1.01 27.07.04 we mp_cmp_d, mp_cmp_mag, mp_cmp
0.1.02 27.07.04 we InitCheck / RunError(210)
0.1.03 27.07.04 we s_mp_add/sub, mp_add/sub
0.1.04 28.07.04 we mp_rand, bugfix: mp_mul_d, s_mp_add
0.1.05 28.07.04 we $Q- in s_mp_sub/mp_sub_d
0.1.06 29.07.04 we s_mp_mul_digs {+tbd}, mp_mul {+tbd}
0.1.07 29.07.04 we mp_count_bits, mp_div
0.1.08 30.07.04 we mp_abs, mp_chs, mp_2expt
0.1.09 30.07.04 we mp_sqr {*tbd}, mp_expt, mp_expt_d
0.1.10 01.08.04 we mp_show_plus, mp_radix_size, mp_toradix
0.1.11 01.08.04 we mp_uppercase, MAXRadix, InitCheck->mp_argchk
0.1.12 01.08.04 we RunError with MP_RTE_xx constants
0.1.13 01.08.04 we mp_cmp_z, mp_cmp_int
0.1.14 01.08.04 we fix compiler warnings, formatting cleanup
0.1.15 01.08.04 we mp_sqrt, mp_mod
0.1.16 02.08.04 we mp_karatsuba_mul, mp_mul completed
0.1.17 02.08.04 we Bug? in Delphi mem mgt, fixed with mp_realloc in mp_grow
0.1.18 02.08.04 we Bug fixed 16 bit mp_realloc, mp_abs interfaced
0.1.19 03.08.04 we Fix Delphi32 with signs in Karatsuba
0.1.20 03.08.04 we BugFix mp_div (negative zero)
0.1.21 03.08.04 we small change in s_mp_mul_digs
0.1.22 04.08.04 we conditional define: use move in Karatsuba
0.1.23 04.08.04 we s_mp_sqr, mp_sqr_cutoff
0.1.24 04.08.04 we mp_expt_d with b: word
0.1.25 04.08.04 we mp_clear with fillchar(a,sizeof(a),0)
0.1.26 04.08.04 we {$i mp_conf.inc), MP_VERSION from mp_types
0.1.27 04.08.04 we mp_karatsuba_sqr, mp_conf.inc removed
0.1.28 05.08.04 we 32Bit fix s_mp_sqr
0.1.29 06.08.04 we Bug search cleanups marked {*0.1.29}
0.1.30 06.08.04 we "Karatsuba" bug fixed in s_mp_add
0.1.31 07.08.04 we last digit optimization in mp_expt, 0 is even
0.2.00 08.08.04 we "release" version 0.2
0.2.01 08.08.04 we Debug: check mp_digits <= MP_DIGIT_MAX
0.2.02 08.08.04 we mp_unsigned_bin_size, mp_radix_size with size: longint
0.2.03 09.08.04 we mp_to_(un)signed_bin_n, mp_read_(un=signed_bin
0.2.04 09.08.04 we mp_mod_d, changed mp_div_d parameter list
0.2.05 09.08.04 we mp_rand renamed to mp_rand_radix, new mp_rand
0.2.04 09.08.04 we bugfix mp_div_2d
0.2.05 09.08.04 we optimize mp_init_copy
0.2.06 09.08.04 we LTM031 change applied to mp_2expt
0.2.07 10.08.04 we mp_init/clear_multi(_p)
0.2.08 11.08.04 we easy outs for mp_expt(_d)
0.2.09 14.08.04 we mp_init_set_int, mp_get_int
0.2.10 18.08.04 we optimized mp_mod_d
0.2.11 24.08.04 we mp_mul_int
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_radix_size: longint, remove D6+ warnings
0.3.03 27.08.04 we mp_result
0.3.04 30.08.04 we mp_radix_size with ln
0.3.05 31.08.04 we mp_is_<xy> functions, mp_memused, lograd table
0.3.06 27.02.05 we mp_freemem in mp_clear, mp_memused with mp_memstat
0.3.07 16.05.05 we FPC: $mode objfpc/$goto on
0.3.08 09.08.05 we Bugfix mp_toradix_n if a has only one digit
0.3.09 15.08.05 we mp_isbit
0.3.10 17.08.05 we severe bug fix mp_expt(a,b,c) for @b=@c
0.4.00 20.08.05 we use mp_set_error
0.4.01 20.08.05 we some mp_digit typecasts, radix/digits word in mp_rand_radix
0.4.02 21.08.05 we mp_conf.inc, mp_memused to mp_base
0.4.03 21.08.05 we MPC_ArgCheck, MPC_Assert
0.4.04 22.08.05 we removed mp_argcheck, mp_errchk
0.4.05 22.08.05 we mp_init<i>, mp_clear<i>, i=2..5
0.4.06 22.08.05 we mp_2expt with longint b
0.4.07 22.08.05 we bugfix mp_init5
0.4.08 23.08.05 we MPC_HaltOnArgCheck
0.4.09 23.08.05 we IsPrime32, IsSPP32, IsSPP32A
0.4.10 24.08.05 we use MaxDigits related codes MP_MAXDIGITS, MP_RTE_OTHER
0.4.11 24.08.05 we fix IsPrime32 for BP7/{$R+,Q+}
0.4.12 24.08.05 we fix mp_init_set/_int
0.4.13 26.08.05 we uses mp_prng (mp_rand, mp_rand_radix)
0.4.14 27.08.05 we argcheck for digit in some mp_xxx_d functions
0.4.15 27.08.05 we functions mp_cmp_mag_d, mp_read_decimal
0.4.16 28.08.05 we usage of exceptions implemented
0.4.17 29.08.05 we mp_set_short, mp_isone, mp_is0, mp_is1
0.4.18 07.09.05 we mp_n_root, bugfix mp_mul_d
0.4.19 08.09.05 we some functions moved from mp_supp
0.4.20 10.09.05 we mp_reduce_2k functions
0.4.21 10.09.05 we mp_core routines integrated
0.4.22 10.09.05 we optimized mp_reduce_is_2k
0.4.23 16.09.05 we optimized mp_cmp_mag and s_mp_sub_d
0.4.24 18.09.05 we mp_let (alternative for mp_read_decimal)
0.4.25 21.09.05 we use mp_clearzero
0.4.26 21.09.05 we $argcheck pc<>pd in mp_div
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 02.10.05 we changed 'SPP' to more conventional 'spsp'
0.5.03 07.10.05 we is_spsp32A: optimized BASM16
0.5.04 08.10.05 we more is_spsp32A optimization
0.5.05 13.10.05 we mp_isbit with BASM16
0.5.06 14.10.05 we BugFix is_spsp32A: reduce a := bases[k] mod N in
0.5.07 15.10.05 we BugFix is_spsp32A: no mod N, use MulMod32 or div
0.5.08 15.10.05 we BugFix N>=9080191 in IsPrime32
0.5.09 20.11.05 we ArgCheck in mp_isbit, new name: mp_montgomery_calcnorm
0.5.10 21.11.05 we mp_rand_bits, mp_rand: check digits<MAXDigits
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 31.12.05 we mp_set_int via array[0..3] of byte
0.6.03 31.12.05 we mp_mul_w, mp_set_w
0.6.04 31.12.05 we popcount16/32, mp_popcount
0.6.05 10.01.06 we Halt on BAD_ARG in mp_montgomery_setup
0.6.06 28.01.06 we s_mp_add_d; changes in mp_sub_d, s_mp_sqr, s_mp_chs
0.7.00 19.03.06 we mp_div_w
0.7.01 04.08.06 we Bugfix in mp_mod_2d, improve mp_reduce_2k
0.7.02 08.08.06 we mp_makeodd
0.7.03 09.08.06 we mp_clear6, mp_init6; removed mp_let, mp_isone
0.7.04 10.08.06 we mp_inc, mp_dec
0.7.05 11.08.06 we bugfix mp_set_short
0.7.06 11.08.06 we avoid FPC warnings: mp_set_int/mp_set_w
0.7.07 11.08.06 we fixed and improved mp_reduce
0.7.08 12.08.06 we s_mp_mul_high_digs: error if id<0
0.7.09 13.08.06 we mp_set_pow, mp_expt_int, mp_expt_d uses mp_expt_int
0.7.10 15.08.06 we mp_[x]_int with [x]: add,dec,div,inc,mod,sub
0.7.11 26.08.06 we rewrite is_spsp32A, new internal function _spsp32
0.7.12 27.08.06 we BIT32: bigalloc; FPC: use ReturnNilIfGrowHeapFails
0.7.13 28.08.06 we mp_n_root with longint parameter
0.7.14 28.08.06 we mp_clrbit, mp_setbit, mp_2expt uses mp_setbit
0.7.15 30.08.06 we mp_gr_mod, mp_gr_setup
0.7.16 30.08.06 we mp_shr, mp_shl, mp_div_2d(..,nil) replaced by mp_shr
0.7.17 06.09.06 we better initial approximation for mp_sqrt
0.7.18 07.09.06 we mp_mod_w
0.8.00 18.09.06 we mp_alloc interfaced
0.9.00 26.12.06 we some minor changes related to mp_clamp/assert
0.9.01 27.12.06 we mp_lshd2
0.9.02 27.12.06 we s_mp_div: Knuth's q calculation from Alg. D
0.9.03 28.12.06 we s_mp_div: separate treatment of single digit b
0.9.04 29.12.06 we $ifdef MPC_USE_Assert
0.9.05 01.01.07 we mp_prod_int, bigalloc renamed to IAlloc
0.9.06 01.01.07 we mp_mul_int optimized for small longints
0.9.07 01.01.07 we mp_mul optimized for single digit factors
0.9.08 02.01.07 we mp_toradix10_n, speed up 2.5 .. 4 for mp_toradix(..,10,..)
0.9.09 03.01.07 we changed mp_read_radix to accept uppercase and lowercase
0.9.10 03.01.07 we new mp_toradix_n (generalization of mp_toradix10_n)
0.9.11 03.01.07 we mp_toradix_n: local TRadixCMap
0.9.12 04.01.07 we bugfixed/improved mp_is_power_of_two, renamed to mp_is_pow2_d
0.9.13 04.01.07 we mp_is_pow2
0.9.14 07.01.07 we improved mp_read_radix
1.0.00 11.04.07 we mp_init_prim: use mp_precision if size=0
1.0.01 11.04.07 we mp_mul_d/w with factor w/d=0/1
1.0.02 12.04.07 we s_mp_read_radix
1.0.03 12.04.07 we Bugfix EstimateQDigit, improved normalization in s_mp_div
1.0.04 13.04.07 we Easy outs in mp_expt_int
1.0.05 14.04.07 we s_mp_toradix_n, off by 1 bugfix mp_radix_astr
1.0.06 15.04.07 we s_mp_write_radix
1.0.07 16.04.07 we mp_is1 without mp_cmp_d
1.0.08 01.05.07 we mp_todouble, ldexpd
1.0.09 01.05.07 we DblPosInf,DblNegInf,DblNaN; improved mp_todouble
1.0.10 05.05.07 we s_mp_read_radix with SignAllowed, bugfix mp_set_w/int
1.0.11 05.05.07 we rewrite mp_set_w; frexpd, bugfix mp_todouble
1.0.12 07.05.07 we 0^0=1 in mp_expt
1.0.13 11.05.07 we Removed MulMod32, mp_shrink
1.0.14 11.05.07 we Bugfix s_mp_write_radix
1.0.15 13.05.07 we Corrected some exception strings
1.0.16 13.05.07 we MPAF prefix in assert strings
1.1.00 27.06.07 we mp_abs: Arg check done in mp_copy
1.1.01 01.07.07 we EstimateQDigit: removed = from q>=MP_MASK test
1.1.02 10.07.07 we mp_writeln, mp_clear[x]/mp_init[x] (x=7..9)
1.1.03 15.07.07 we mp_reduce: easy out if x<m, allow x<0
1.1.04 21.07.07 we improved mp_shl, easy out in mp_shr
1.1.05 21.07.07 we mp_sqrt: Initial double approximation based on highest mp_digit(s)
1.1.06 26.07.07 we isqrt32
1.1.07 26.07.07 we isqrt32 with FPU
1.1.08 26.07.07 we mp_sqrt: new recursive integer square root algorithm
1.1.09 29.07.07 we improved: mp_reduce, s_mp_mul_high_digs, s_mp_mul_digs
1.1.10 04.08.07 we renamed/new: s_mp_mod_2d, mp_mod_2d
1.2.00 17.08.07 we gcd32/gcd32u
1.2.01 19.08.07 we changed mp_mod_int to use mp_mod
1.2.02 19.08.07 we Bugfix mp_mod: don't adjust sign if result=0
1.2.03 19.08.07 we mp_mod_int without temporary mp_int, local BASM in loop
1.2.04 26.08.07 we gcd32/gcd32u call internal __gcd32
1.2.05 02.09.07 we new s_mp_mul_int used in mp_prod_int
1.2.06 04.09.07 we mp_rand_bits_ex
1.2.07 05.09.07 we s_mp_expt_dl, s_mp_expt_wl, changed mp_set_pow
1.2.08 05.09.07 we MP_32BIT/MP_16BIT versions of mp_set_int
1.2.09 06.09.07 we popcount16 and popcount32 return integer
1.2.10 07.09.07 we mp_checksum, s_mp_checksum
1.2.11 10.09.07 we use MP_INV_MASK to avoid warnings if DIGIT_BIT=31
1.2.12 11.09.07 we Bugfix in BIT16 version of _spsp32
1.2.13 17.09.07 we IAlloc, mp_alloc, mp_freemem inline $ifdef HAS_INLINE
1.2.14 17.09.07 we mp_init_multi/_p use mp_init_prim, mp_abs inline
1.2.15 17.09.07 we Bugfix s_mp_read_radix for DIGIT_BIT<10
1.2.16 20.09.07 we Removed inline (D9 bug with mp_freemem(pointer(pchar)...)
1.2.17 21.09.07 we mp_sqrt uses at most one local mp_int
1.3.00 03.11.07 we mp_toextended, frexpx, ldexpx
1.3.01 05.11.07 we mp_todouble_ex, mp_toextended_ex
1.3.02 11.11.07 we s_mp_read_radix: sep now a string
1.3.03 14.11.07 we mp_radix_astr: prefill result with #0
1.3.04 19.11.07 we merged routines from mp_supp
1.3.05 22.11.07 we complete rewrite of mp_to_unsigned_bin_n,
mp_read_unsigned_bin, and mp_rand
1.3.06 26.11.07 we add32_ovr
1.3.07 26.11.07 we fix mask/bit logic in mp_rand_bits_ex
1.3.08 29.11.07 we Removed some arg checks (done in mp_copy)
1.3.09 09.12.07 we s_mp_add_ovr -> add32_ovr, undo mp_supp merging
1.3.10 17.12.07 we moved lograd table to interface
1.5.00 24.01.08 we mp_and, mp_or, mp_xor
1.5.01 31.01.08 we {$x+} for VP and D1
1.6.00 23.05.08 we Optimized Argcheck in: mp_mod, mp_cmp_int, mp_todouble_ex, mp_toextended_ex
1.6.01 24.05.08 we mp_not_init_multi
1.6.02 25.05.08 we mp_hex/mp_ahex
1.6.03 03.06.08 we mp_is1a, bugfix mp_is1
1.6.04 06.06.08 we mp_is? routines return false if mp_error<>MP_OKAY
1.6.05 08.06.08 we improved carry propagation in s_mp_add_d,s_mp_sub_d
1.6.06 10.06.08 we fix mp_mod_d for a<0
1.6.07 11.06.08 we mp_init_prim: use size=4 if size=mp_allocprec=0
1.7.00 23.08.08 we Avoid FPC222 warning in isqrt32
1.7.01 24.08.08 we s_mp_toradix_n improved if radix is power of 2
1.7.02 14.09.08 we new function IsPow2_w
MP_16BIT: improved mp_div_d, mp_div_w, mp_mod, s_mp_div
1.7.03 14.09.08 we Fix FPC RTE 201 if R+ in mp_read_unsigned_bin
1.7.04 15.09.08 we mp_get/set_allocprec, mp_allocprec local,
mp_init_prim rounds up to multiple of mp_allocprec}
1.7.05 17.09.08 we mp_rand_ex, mp_is_longint; BIT16: (f)LeftShiftAdd
1.7.06 17.09.08 we mp_mod with mp_is_longint, BASM16 for 'shr DIGIT_BIT'
1.7.07 21.09.08 we mp_shr1
1.7.08 24.09.08 we mp_sign, removed mp_cmp_z, renamed ??_2d to ??_2k
1.7.09 24.09.08 we string replaced by mp_string
1.7.10 26.09.08 we mp_read_radix_str, mp_read_decimal_str
1.7.11 27.09.08 we invmod32
1.7.12 02.10.08 we improved mp_reduce
1.8.00 04.10.08 we BASM16 in s_mp_mul_(high)_digs
1.8.01 04.10.08 we BASM16 in s_mp_sqr
1.8.02 04.10.08 we BASM16 in mp_montgomery_reduce
1.8.03 05.10.08 we BASM16 in mp_shl, mp_shr
1.8.04 05.10.08 we BASM16 in mp_div_w and new s_mp_div_d, fixed IsPow2_w
1.8.05 06.10.08 we Simplified MP_32BIT power of two code in mp_div_w
1.8.06 07.10.08 we Check for power of two in mp_mod_int
1.8.07 09.10.08 we mp_set1
1.8.08 11.10.08 we Improved mp_reduce_2k_setup
1.8.09 18.10.08 we Improved mp_expt_int, mp_popcount, mp_gr_setup
1.8.10 19.10.08 we s_mp_mod_w
1.8.11 23.10.08 we mp_n_root with Halley's iteration or bisection method
1.8.12 24.10.08 we mp_n_root: check startup convergence of Halley steps
1.8.13 25.10.08 we s_mp_n_root2, mp_n_root2
1.8.14 26.10.08 we Fix check for exact root in Halley
1.8.15 27.10.08 we Halley: optimized to keep remainder
1.8.16 31.10.08 we Halley: bisection count = bitsize32(d)
1.8.17 01.11.08 we mp_mod_int: check if a is longint
1.8.18 04.11.08 we s_mp_n_root2: b=1 if n>=mp_bitsize(a)
1.9.00 08.11.08 we mp_rshd2, avoid some warnings
1.9.01 08.11.08 we Halley renamed to iroot, improved initial approximation
1.9.02 30.11.08 we mp_read_radix_arr
1.9.03 02.12.08 we Uses BTypes: char8, pchar8
1.9.04 06.12.08 we s_mp_is_le0
1.9.05 26.12.08 we IsPrime16/32,is_spsp32,is_spsp32A move to mp_prime
1.9.06 02.01.09 we s_mp_sqrtrem, mp_sqrtrem, renamed mp_sqrt to s_mp_sqrt
1.9.07 03.01.09 we mp_n_root2 with mp_sqrtrem
1.9.08 06.01.09 we improved mp_reduce_2k
1.9.09 06.01.09 we Replaced CHAR_BIT
1.10.00 06.01.09 we skip final sqr in mp_expt
1.10.01 21.01.09 we mp_shl1, renamed (s)mp_div to (s)mp_divrem, new mp_div
1.10.02 25.01.09 we improved mp_add/dec/inc/sub_int
1.10.03 01.02.09 we s_mp_ln, s_mp_log2, s_mp_set_ext
1.10.04 04.02.09 we s_mp_mod_is0
1.10.05 16.02.09 we mp_shlx, mp_shrx
1.10.06 18.02.09 we mp_divrem_newton
1.11.00 07.03.09 we s_mp_divrem renamed to s_mp_divrem_basecase, new s_mp_divrem with Burnikel/Ziegler
1.11.01 08.03.09 we improved bz_d3n2n, changed bz_divrem_pos argument list
1.11.02 14.03.09 we s_mp_toom3_mul, renamed s_mp_karatsuba_mul/sqr
1.11.03 15.03.09 we s_mp_toom3_mul: use s_mp_mod_2k, mp_shl
1.11.04 15.03.09 we s_mp_toom3_sqr, changed mp_sqr
1.11.05 16.03.09 we s_mp_toom3_mul: split with B ~ 2^(bitsize(max(a,b)/3)
1.11.06 16.03.09 we s_mp_fakeinit (used in s_mp_toom3_mul/sqr)
1.11.07 17.03.09 we complete rewrite of s_mp_karatsuba_mul/sqr
1.11.08 20.03.09 we mp_mul: separate handling of unbalanced factors
1.11.09 21.03.09 we mp_mul: fix sign for unbalanced part
1.11.10 25.03.09 we fix sign in mp_lshd2
1.11.11 25.03.09 we s_mp_toom3_mul uses Bodrato algorithm
1.11.12 26.03.09 we s_mp_fakeinit with mp_clamp
1.11.13 26.03.09 we s_mp_toom3_sqr uses Bodrato algorithm
1.11.14 27.03.09 we mp_clamp without init check, improved fake init in Toom-3
1.11.15 27.03.09 we removed MPC_UseToom3 (and mem check for BP7)
1.11.16 29.03.09 we improved s_mp_mul_high_digs
1.11.17 30.03.09 we mp_init_size2, removed mp_divrem_newton
1.11.18 30.03.09 we reduce temporary memory in s_mp_karatsuba_mul/sqr
1.12.00 20.06.09 we Fix mp_rand_bits_ex: add mp_clamp(a)
1.14.00 13.02.10 we MPC_MAXRadix64 adjustments
1.16.00 05.06.10 we mp_read_decimal_astr, mp_read_radix_astr
1.16.01 25.07.10 we const pa in mp_read_radix_arr
1.19.00 03.11.11 we Memory routines for MAXDigits=32000 with MP_32BIT
1.20.00 12.01.12 we Include files mp_bas64/32/16.inc
1.20.01 13.01.12 we rewrite (ld/fr)exp(x/d)
1.20.02 13.01.12 we adjust mp_toextended_ex for EXT64
1.20.03 15.01.12 we Full length for mp_radix_astr
1.20.04 16.01.12 we New add32_ovr in mp_bas64.inc
1.20.05 17.01.12 we Fix border cases in mp_lshd/2
1.20.06 19.01.12 we Subquadratic mp_radix_astr
1.20.07 20.01.12 we s_mp_write_radix uses mp_radix_astr for BIT32/64
1.20.08 21.01.12 we s_mp_radix_astr
1.21.00 19.07.12 we Improved mp_div_int for 0 < b < MP_DIGIT_MAX
1.21.01 19.07.12 we s_mp_mod_int, mp_mod_d calls s_mp_mod_int if MP_32BIT
1.21.02 21.07.12 we mp_product
1.21.03 24.07.12 we isqrt32 without FPU, new is_square32/ex
1.21.04 25.07.12 we exptmod32/jacobi32 from mp_numth
1.21.05 30.07.12 we Improved mp_set_pow, mp_expt_int: 0^0 = 1
1.22.00 26.08.12 we mp_mul: call mp_sqr if @a=@b
1.22.01 29.08.12 we kronjac32: used by kronecker32 and jacobi32
1.22.02 02.09.12 we mp_radix_astr etc with truncation for BIT16
1.22.03 04.09.12 we icbrt32
1.23.00 23.09.12 we improved BASM16 s_mp_div_d, toom3 with s_mp_div_d
1.23.01 25.09.12 we Fix: uses production code invmod32 in mp_bas16.inc
1.23.02 01.10.12 we xgcd32
1.23.03 02.10.12 we Avoid -0 in s_mp_divrem_basecase after mp_div_w
1.23.04 03.10.12 we mp_set_pow: special code if b is a power of two
1.23.05 04.10.12 we BASM16: improved mp_mul_d
1.23.06 04.10.12 we mp_div_int with s_mp_divrem_basecase
1.23.07 05.10.12 we Removed second s_mp_mod_2k in mp_mod_2k for a<0
1.23.08 05.10.12 we Reintroduce s_mp_shrink (mainly for BP7 real mode)
1.23.09 12.10.12 we Improved s_mp_n_root2: Use Newton if overflow probable, allow deeper recursion
1.24.00 12.12.12 we MPC_PurePascal
1.24.01 15.12.12 we Some word/integer types changed to longint or TNInt
1.24.02 03.01.13 we Fix some more words to TNInt
1.24.03 04.01.13 we Fix mp_set_int for MP_16bit and b=$80000000
1.24.04 04.01.13 we Renamed mp_expt_d to mp_expt_w
1.24.05 05.01.13 we Fix debug mode icbrt32 for VPC, D2/D3
1.24.06 05.01.13 we Fix FPC242/244 bug in debug mode s_mp_div_d
1.30.00 27.09.14 we mp_random
1.30.01 30.09.14 we Reintroduced mulmod32 function
1.30.02 03.10.14 we mlinsolve32
1.36.00 12.09.17 we Fix obscure 16-Bit (MSDOS) FP error in s_mp_n_root2
1.37.00 10.05.18 we Renamed mp_reverse to s_mp_reverse
1.38.00 25.06.18 we frexpx/ldxepx only if not EXT64
1.38.01 26.06.18 we s_mp_set_ext -> s_mp_set_dbl, s_mp_ln/log2 return double
*************************************************************************)
(*-------------------------------------------------------------------------
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
See the file '3rdparty.mpa' for the licenses.
----------------------------------------------------------------------------*)
(*-------------------------------------------------------------------------
(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+}
{----------------------------------}
{32 bit and floating point routines}
{----------------------------------}
{#Z-}
function add32_ovr(x,y: longint; var z: longint): boolean;
{-Add z=x+y with overflow detection}
function bitsize32(a: longint): integer;
{-Return the number of bits in a (index of highest bit), 0 if no bit is set}
function gcd32(A, B: longint): longint;
{-Calculate GCD of two longints}
function gcd32u(A, B: longint): longint;
{-Calculate GCD of two longints (DWORD interpretation)}
function exptmod32(a,b,c: longint): longint;
{-Calculate a^b mod c if a>=0, b>=0, c>0; result=0 otherwise}
function icbrt32(a: longint): longint;
{-Return the integer cube root sign(a)*floor(|a|^(1/3))}
function invmod32(a,b: longint): longint;
{-Return a^-1 mod b, b>1. Result is 0 if gcd(a,b)<>1 or b<2}
function isqrt32(a: longint): longint;
{-Return floor(sqrt(abs(a))}
function is_square32(a: longint): boolean;
{-Test if a is square, i.e. test if a=sqr(isqrt32(a)), false if a<0}
function is_square32ex(a: longint; var b: longint): boolean;
{-Test if a is square, false if a<0. If yes, b = sqrt(a) else b is undefined}
function jacobi32(a,b: longint): integer;
{-Compute the Jacobi/Legendre symbol (a|b), b: odd and > 2}
function kronecker32(a,b: longint): integer;
{-Compute the Kronecker symbol (a|b)}
function mlinsolve32(a,b,m: longint; var x0,d: longint): boolean;
{-Solve ax=b mod m; m>1. If result=true then x0 >= 0, d=gcd(a,m) > 0}
{ and there are d solutions: x_k = x0 + k*(m/d) mod m, for k=0..d-1.}
function mulmod32(a,b,n: longint): longint;
{-Return a*b mod n, assumes n>0, a,b>=0}
function popcount16(w: word): integer;
{-Get population count = number of 1-bits in a word}
function popcount32(l: longint): integer;
{-Get population count = number of 1-bits in a longint}
procedure xgcd32(a,b: longint; var u1,u2,u3: longint);
{-Extended gcd algorithm, calculate a*u1 + b*u2 = u3 = gcd(a,b); a,b <> -2^31}
function DblPosInf: double;
{-Return positive double infinity}
function DblNegInf: double;
{-Return negative double infinity}
function DblNaN: double;
{-Return double NaN (Not a Number)}
procedure frexpd(d: double; var m: double; var e: longint);
{-Return m,e with d=m*2^e and 0.5 <= abs(m) < 1}
function ldexpd(d: double; e: longint): double;
{-Return d*2^e}
{$ifndef EXT64}
procedure frexpx(x: extended; var m: extended; var e: longint);
{-Return m,e with x=m*2^e and 0.5 <= abs(m) < 1}
function ldexpx(x: extended; e: longint): extended;
{-Return x*2^e}
{$endif}
{#Z+}
{---------------------------------------}
{mp functions are sorted alphabetically }
{---------------------------------------}
{#Z-}
procedure mp_2expt(var a: mp_int; b: longint);
{-Compute a = 2^b, b>=0, error if b<0 or b>MP_MAXBIT}
procedure mp_abs(const a: mp_int; var b: mp_int);
{-Absolute value, b = |a|}
procedure mp_add(const a,b: mp_int; var c: mp_int);
{-High level addition (handles signs)}
procedure mp_add_d(const a: mp_int; b: mp_digit; var c: mp_int);
{-Single digit addition}
procedure mp_add_int(const a: mp_int; b: longint; var c: mp_int);
{-Calculate c = a + b}
function mp_adecimal(const a: mp_int): ansistring;
{-Convert to decimal ansistring}
function mp_ahex(const a: mp_int): ansistring;
{-Convert to hex ansistring}
function mp_radix_astr(const a: mp_int; radix: word): ansistring;
{-Convert to radix representation ansistring (set mp_show_plus to show '+')}
procedure mp_read_decimal_astr(var a: mp_int; const s: ansistring);
{-Read an mp_int from a decimal ansistring}
procedure mp_read_radix_astr(var a: mp_int; const s: ansistring; radix: word);
{-Read an mp_int from an ansistring in given radix}
{$ifdef BIT32or64}
function s_mp_radix_astr(const a: mp_int; radix: word; plus: boolean): ansistring;
{-Convert to radix representation ansistring, plus to show '+'}
{$endif}
procedure mp_and(const a,b: mp_int; var c: mp_int);
{-Calculate c = a and b}
function mp_bitsize(const a: mp_int): longint;
{-Return the number of bits in a (index of highest bit), 0 if no bit is set}
function mp_checksum(const a: mp_int): longint;
{-Return a checksum for a, -1 if mp_error<>MP_OKAY, -2 if not initialized}
procedure mp_chs(const a: mp_int; var b: mp_int);
{-Change sign, b = -a}
procedure mp_clear(var a: mp_int);
{-Free an mp_int}
procedure mp_clear2(var a,b: mp_int);
{-Clear 2 mp_ints}
procedure mp_clear3(var a,b,c: mp_int);
{-Clear 3 mp_ints}
procedure mp_clear4(var a,b,c,d: mp_int);
{-Clear 4 mp_ints}
procedure mp_clear5(var a,b,c,d,e: mp_int);
{-Clear 5 mp_ints}
procedure mp_clear6(var a,b,c,d,e,f: mp_int);
{-Clear 6 mp_ints}
procedure mp_clear7(var a,b,c,d,e,f,g: mp_int);
{-Clear 7 mp_ints}
procedure mp_clear8(var a,b,c,d,e,f,g,h: mp_int);
{-Clear 8 mp_ints}
procedure mp_clear9(var a,b,c,d,e,f,g,h,i: mp_int);
{-Clear 9 mp_ints}
procedure mp_clear_multi(var vi: array of mp_int);
{-Clear a vector of mp_ints}
procedure mp_clear_multi_p(const pv: array of pmp_int);
{-Clear a list of mp_ints given as a ptr vector}
procedure mp_clrbit(var a: mp_int; n: longint);
{-Clear bit n of a, no action if out of range, (1 = bit 0)}
function mp_cmp(const a,b: mp_int): integer;
{-Compare two mp_ints (signed), return sign(a-b)}
function mp_cmp_d(const a: mp_int; b: mp_digit): integer;
{-Compare a with an mp_digit, return sign(a-b)}
function mp_cmp_int(const a: mp_int; b: longint): integer;
{-Compare a with a longint, return sign(a-b)}
function mp_cmp_mag(const a,b: mp_int): integer;
{-Compare magnitude of two mp_ints (unsigned), return sign(|a|-|b|)}
function mp_cmp_mag_d(const a: mp_int; b: mp_digit): integer;
{-Compare |a| with a digit, return sign(|a|-b)}
function mp_cnt_lsb(const a: mp_int): longint;
{-Count the number of least significant bits which are zero}
procedure mp_copy(const a: mp_int; var b: mp_int);
{-Copy an mp_int, b = a}
procedure mp_dec(var a: mp_int);
{-Decrement an mp_int by 1}
procedure mp_dec_int(var a: mp_int; b: longint);
{-Calculate a = a - b}
function mp_decimal(const a: mp_int): mp_string;
{-Convert to decimal, max 255 digits}
procedure mp_div(const a,b: mp_int; var c: mp_int);
{-Integer signed division, c = a div b}
procedure mp_divrem(const a,b: mp_int; pc,pd: pmp_int);
{-Integer signed division, pc^ = a div b, pd^ = a rem b; sign(pd^)=sign(a)}
procedure mp_div_2(const a: mp_int; var b: mp_int);
{-Divide by 2, b = a/2}
procedure mp_div_2k(const a: mp_int; b: longint; var c: mp_int; pd: pmp_int);
{-Divide by 2^b; quotient in c, optional remainder in pd^, sign(pd^)=sign(a)}
procedure mp_div_d(const a: mp_int; b: mp_digit; pc: pmp_int; var d: mp_digit);
{-Single digit division, pc^ = a div b, d = a mod b}
procedure mp_div_int(const a: mp_int; b: longint; pc: pmp_int; var d: longint);
{-Integer signed division, pc^ = a div b, d = a rem b; sign(d)=sign(a)}
procedure mp_div_w(const a: mp_int; b: word; pc: pmp_int; var r: word);
{-Divide a by a single word b, pc^=sign(a)(|a| div b), r = |a| mod b}
procedure mp_exch(var a,b: mp_int);
{-Exchange two mp_ints}
procedure mp_expt(const a,b: mp_int; var c: mp_int);
{-Calculate c = a^b, b>=0}
procedure mp_expt_w(const a: mp_int; b: word; var c: mp_int);
{-Calculate c = a^b}
procedure mp_expt_int(const a: mp_int; b: longint; var c: mp_int);
{-Calculate c = a^b, b>=0}
function mp_gcd_int(const a: mp_int; b: longint): longint;
{-Return gcd(a,b), b<>0}
function mp_get_int(const a: mp_int): longint;
{-Get the lower signed 31 bits of an mp_int}
procedure mp_gr_mod(var x: mp_int; const N, R: mp_int);
{-Reduce x to x mod N, N > 1, using generalized reciprocal iteration.}
{ Result is >= 0. R is from mp_gr_setup. Compared to the similar}
{ Barrett reduction the restricted range 0<x<N^2 is not required.}
procedure mp_gr_setup(var RN: mp_int; const N: mp_int);
{-Calculate the generalized reciprocal for N>0}
function mp_hex(const a: mp_int): mp_string;
{-Convert to hex string, max 255 digits}
procedure mp_inc(var a: mp_int);
{-Increment an mp_int by 1}
procedure mp_inc_int(var a: mp_int; b: longint);
{-Calculate a = a + b}
procedure mp_init(var a: mp_int);
{-Initialize an mp_int}
procedure mp_init2(var a,b: mp_int);
{-Initialize 2 mp_ints}
procedure mp_init3(var a,b,c: mp_int);
{-Initialize 3 mp_ints}
procedure mp_init4(var a,b,c,d: mp_int);
{-Initialize 4 mp_ints}
procedure mp_init5(var a,b,c,d,e: mp_int);
{-Initialize 5 mp_ints}
procedure mp_init6(var a,b,c,d,e,f: mp_int);
{-Initialize 6 mp_ints}
procedure mp_init7(var a,b,c,d,e,f,g: mp_int);
{-Initialize 7 mp_ints}
procedure mp_init8(var a,b,c,d,e,f,g,h: mp_int);
{-Initialize 8 mp_ints}
procedure mp_init9(var a,b,c,d,e,f,g,h,i: mp_int);
{-Initialize 9 mp_ints}
procedure mp_init_copy(var a: mp_int; const b: mp_int);
{-Create a, then copy b into it}
procedure mp_init_multi(var vi: array of mp_int);
{-Initialize a vector of mp_ints}
procedure mp_init_multi_p(var pv: array of pmp_int);
{-Initialize a list of mp_ints given as a ptr vector}
procedure mp_init_set(var a: mp_int; b: mp_digit);
{-Initialize and set a digit}
procedure mp_init_set_int(var a: mp_int; b: longint);
{-Initialize and set a to a longint}
procedure mp_init_size(var a: mp_int; size: longint);
{-Initialize a to size digits, rounded up to multiple of mp_allocprec}
procedure mp_init_size2(var a,b: mp_int; size: longint);
{-Initialize a and b to size digits, rounded up to multiple of mp_allocprec}
function mp_isbit(const a: mp_int; n: longint): boolean;
{-Test if bit n of a is set, (1 = bit 0)}
function mp_iseven(const a: mp_int): boolean;
{-Initialized and even}
function mp_isodd(const a: mp_int): boolean;
{-Initialized and odd}
function mp_is_eq(const a,b: mp_int): boolean;
{-Return a = b}
function mp_is_ge(const a,b: mp_int): boolean;
{-Return a >= b}
function mp_is_gt(const a,b: mp_int): boolean;
{-Return a > b}
function mp_is_le(const a,b: mp_int): boolean;
{-Return a <= b}
function mp_is_lt(const a,b: mp_int): boolean;
{-Return a < b}
function mp_is_ne(const a,b: mp_int): boolean;
{-Return a <> b}
function mp_is_longint(const a: mp_int; var b: longint): boolean;
{-Test if a fits into longint, if true set b := a}
function mp_is_pow2(const a: mp_int; var n: longint): boolean;
{-Check if |a| is a power of 2, if true, return n with |a|=2^n}
function mp_is_pow2_d(d: mp_digit; var n: integer): boolean;
{-Check if d is power of 2, if true, return n with d=2^n}
function IsPow2_w(w: word; var n: integer): boolean;
{-Check if w is power of 2, if true, return n with w=2^n}
function mp_iszero(const a: mp_int): boolean;
{-Initialized and zero}
function mp_is0(const a: mp_int): boolean;
{-Initialized and = 0}
function mp_is1(const a: mp_int): boolean;
{-Initialized and a = 1}
function mp_is1a(const a: mp_int): boolean;
{-Initialized and abs(a) = 1}
procedure mp_lshd(var a: mp_int; b: longint);
{-Shift left a certain amount of digits}
procedure mp_lshd2(const a: mp_int; var b: mp_int; cnt: longint);
{-Set b to a shifted left by cnt digits}
procedure mp_makeodd(const a: mp_int; var b: mp_int; var s: longint);
{-Return b,s with a = 2^s*b if a<>0, b=0,s=-1 otherwise}
procedure mp_mod(const a,b: mp_int; var c: mp_int);
{-Calculate c = a mod b, 0 <= c < b}
procedure mp_mod_2k(const a: mp_int; b: longint; var c: mp_int);
{-Calculate c = a mod 2^b, 0 <= c < 2^b}
procedure mp_mod_d(const a: mp_int; b: mp_digit; var c: mp_digit);
{-Calculate c = a mod b, 0 <= c < b (digit version)}
procedure mp_mod_int(const a: mp_int; b: longint; var c: longint);
{-Calculate c = a mod b}
procedure mp_mod_w(const a: mp_int; b: word; var r: word);
{-Calculate r = a mod b for a single word b}
procedure mp_montgomery_calcnorm(var R: mp_int; const m: mp_int);
{-Calculate R = B^n mod m, n=number of digits in m, B=2^DIGIT_BIT}
procedure mp_montgomery_reduce(var x: mp_int; const n: mp_int; rho: mp_digit);
{-Calculate x = xR^-1 (mod n) via Montgomery reduction}
procedure mp_montgomery_setup(const n: mp_int; var rho: mp_digit);
{-Calculate rho = -1/n mod B for Montgomery reduction, B=2^DIGIT_BIT}
procedure mp_mul(const a,b: mp_int; var c: mp_int);
{-High level multiplication, c = a*b}
procedure mp_mul_2(const a: mp_int; var b: mp_int);
{-Multiply by 2, b = 2*a}
procedure mp_mul_2k(const a: mp_int; b: longint; var c: mp_int);
{-Shift left a, c = a*2^b; c=a if b<=0 [synonym for mp_shl]}
procedure mp_mul_d(const a: mp_int; b: mp_digit; var c: mp_int);
{-Multiply by a digit}
procedure mp_mul_int(const a: mp_int; b: longint; var c: mp_int);
{-Multiply by a 32 bit integer}
procedure mp_mul_w(const a: mp_int; b: word; var c: mp_int);
{-Multiply by a word}
function mp_not_init(const a: mp_int): boolean;
{-Sanity check if a is initialized, does not catch all cases!}
function mp_not_init_multi(const a: array of mp_int): boolean;
{-Sanity check if all elements of a are initialized, does not catch all cases!}
procedure mp_n_root(const a: mp_int; n: longint; var b: mp_int);
{-Calculate the n'th root of a, a must >=0 if n is even; b=0 if n<1}
procedure mp_n_root2(const a: mp_int; n: longint; var b: mp_int; pr: pmp_int);
{-Calculate the n'th root of a, pr^=a-b^n, a must be >=0 if n is even; b=0,pr^=0 if n<1}
procedure mp_output_decimal(const a: mp_int);
{-Write decimal representation to output}
procedure mp_output_radix(const a: mp_int; radix: word);
{-Write radix representation to output}
procedure mp_or(const a,b: mp_int; var c: mp_int);
{-Calculate c = a or b}
function mp_popcount(const a: mp_int): longint;
{-Get population count = number of 1-bits in a}
procedure mp_product(var a: mp_int; n,m: longint);
{-Return the product a = n*(n+1)*...(m-1)*m, a=1 if m<n}
procedure mp_prod_int(var a: mp_int; const b: array of longint; n: longint);
{-Calculate a = product of first n elements of longint array b}
function mp_radix_size(const a: mp_int; radix: word): longint;
{-Return size of ASCII representation (incl. sign and #0)}
function mp_radix_str(const a: mp_int; radix: word): mp_string;
{-Convert to radix representation, max 255 digits}
procedure mp_rand(var a: mp_int; digits: longint);
{-Make a pseudo-random mp_int of a given digit size}
procedure mp_random(const a: mp_int; var b: mp_int);
{-Make a pseudo-random mp_int b with 0 <= b < |a|}
procedure mp_rand_ex(var a: mp_int; digits: longint; sethi: boolean);
{-Make a pseudo-random mp_int of a given digit size, if not sethi}
{ then a[digits-1] may be zero (and a.used will be decremented)}
procedure mp_rand_bits(var a: mp_int; bits: longint);
{-Make a pseudo-random mp_int of a given bit size}
procedure mp_rand_bits_ex(var a: mp_int; bits: longint; sethi: boolean);
{-Make pseudo-random a with bitsize <= bits, if sethi highest bit is set}
procedure mp_rand_radix(var a: mp_int; radix: word; digits: longint);
{-Make a pseudo-random mp_int of order radix^digits}
procedure mp_read_decimal(var a: mp_int; str: pchar8);
{-Read an mp_int from a decimal ASCII pchar}
procedure mp_read_decimal_str(var a: mp_int; const s: mp_string);
{-Read an mp_int from a decimal ASCII string[255]}
procedure mp_read_radix(var a: mp_int; str: pchar8; radix: word);
{-Read an mp_int from a ASCII pchar in given radix}
procedure mp_read_radix_arr(var a: mp_int; const pa: array of pchar8; radix: word);
{-Read an mp_int from concatenated pchar array pa in given radix,}
{ max 65000 chars. Mainly used for 16 bit compatibility with max.}
{ length of string literals = 255 and line length = 127 chars }
procedure mp_read_radix_str(var a: mp_int; const s: mp_string; radix: word);
{-Read an mp_int from a ASCII string[255] in given radix}
procedure mp_read_signed_bin(var a: mp_int; const b; numbytes: longint);
{-Read signed bin, big endian, first byte is 0=positive or 1=negative}
procedure mp_read_unsigned_bin(var a: mp_int; const b; numbytes: longint);
{-Reads a unsigned mp_int, assumes the msb is stored first [big endian]}
procedure mp_reduce(var x: mp_int; const m, mu: mp_int);
{-Reduce x mod m via Barrett, assumes x<m^2, mu is from mp_reduce_setup}
procedure mp_reduce_setup(var mu: mp_int; const a: mp_int);
{-Pre-calculate the value required for Barrett reduction}
procedure mp_reduce_2k(var a: mp_int; const n: mp_int; d: mp_digit);
{-Reduce a mod n where n is of the form 2^p-d, @a<>@n}
procedure mp_reduce_2k_setup(const a: mp_int; var d: mp_digit);
{-Determine setup value d for unrestricted diminished radix reduction, a>=0}
function mp_reduce_is_2k(const a: mp_int): boolean;
{-Determine if mp_reduce_2k can be used}
function mp_result: integer;
{-Return and reset mp_error}
procedure mp_rshd(var a: mp_int; b: longint);
{-Shift right a certain amount of digits}
procedure mp_rshd2(const a: mp_int; var b: mp_int; cnt: longint);
{-Set b to a shifted right by cnt digits}
procedure mp_set(var a: mp_int; b: mp_digit);
{-Set a to digit b}