function NRootMod2k_2adic(
var A: IInteger;
const B,E: IInteger; K: Cardinal): Boolean;
// A = B^(1/E) mod 2^k in 2-adic metric such that (a^e * b) mod 2^k == 1 mod 2^k
var
R,S,T,U: IInteger;
procedure BRootMod2k_2adic(K: Cardinal);
// recursive newton method
begin
if K = 1
then NSet(R, 1)
else
begin
BRootMod2k_2adic((K +1)
shr 1);
NPowMod2k(T, R, S, K);
// T = R^(E +1) * B mod 2^k
NMulMod2k(T, B, K);
NMulMod2k(R, S, K);
// R = R * (E +1) - T
NSub(R, T);
NMulMod2k(R, U, K);
// R = R div E mod 2^k = R * C^-1 mod 2^k
end;
end;
resourcestring
sNRootMod2k = '
NRootMod2k(), Parameter must be > 0 and K > 0 and E > 0 and E = 2 or odd';
var
I: Integer;
begin
Result := False;
I := NSgn(E, True);
if (K = 0)
or (I <= 0)
or (NSgn(B) <= 0)
then NRaise(@sNRootMod2k);
if I = 1
then Result := NInvMod2k(A, B, K)
else
if I = 2
then Result := NSqrtMod2k_2adic(A, B, K)
else
if not NOdd(E)
then NRaise(@sNRootMod2k)
else
if K = 1
then
begin
Result := True;
NSet(A, 1);
end else
begin
Result := NOdd(B);
if Result
then
begin
NInvMod2k(U, E, K);
// U = E^-1 mod 2^k
NSet(S, E);
NInc(S);
// S = E +1
BRootMod2k_2adic(K);
NSwp(R, A);
end else NSet(A, 0);
end;
end;
function NRootMod2k(
var A: IInteger;
const E: IInteger; K: Cardinal): Boolean;
begin
Result := NRootMod2k(A, A, E, K);
end;
function NRootMod2k(
var A: IInteger;
const B,E: IInteger; K: Cardinal): Boolean;
// A = B^(1/E) mod (2^k)
begin
if NSgn(E, True) = 1
then
begin
NCut(A, B, K);
if PInt(A).FNeg
then NCpl(A, K, True);
Result := PInt(A).FCount <> 0;
end else
if NSgn(E, True) = 2
then Result := NSqrtMod2k(A, B, K)
else Result := NRootMod2k_2adic(A, B, E, K)
and NInvMod2k(A, K);
end;
function NIsPower(
const N: IInteger; B,E: Integer): Boolean;
overload;
// result := N = |B^E|
begin
if (Abs(B) = 2)
and (E >= 0)
then Result := NIsPowerOf2(N) = E
else Result := NIsPower(N, NInt(B), NInt(E));
end;
function NIsPower(
const N,B: IInteger; E: Integer): Boolean;
overload;
begin
if (NCmp(B, 2, True) = 0)
and (E >= 0)
then Result := NIsPowerOf2(N) = E
else Result := NIsPower(N, B, NInt(E));
end;
function NIsPower(
const N,B,E: IInteger): Boolean;
overload;
// result := N = |B^E|
var
K,S: Integer;
R,T: IInteger;
begin
Result := False;
K := 64;
S := NSize(N);
if NSgn(B) < 0
then NSet(T, B, True)
else T := B;
repeat
if K > S
then K := S;
if not NPowMod2k(R, T, E, K)
or (NCmp(R, N, K, True) <> 0)
then Exit;
if K >= S
then Break;
K := K * 8;
until False;
Result := True;
end;
function NIsPowerOf2(
const A: IInteger): Integer;
var
H: Cardinal;
begin
Result := NSize(A) -1;
if Result >= 0
then
begin
H := PInt(A).FNum[PInt(A).FCount -1];
if (H
and (H -1) <> 0)
or (NLow(A) <> Result)
then
Result := -1;
end;
end;
function NCheckSqrTable(R: Byte): Boolean;
const
T256:
array[0..7]
of Cardinal = ($02030213,$02020212,$02020213,$02020212,
$02030212,$02020212,$02020212,$02020212);
begin
Result := Odd(T256[R
div 32]
shr (R
mod 32));
end;
function NCheckSqr(
const A: IInteger): Boolean;
const
T193:
array[0.. 6]
of Cardinal = ($645AAC20,$3638B3EE,$4F85F6D4,$ADBE87C8,
$DF3471B0,$10D56899,$FFFFFFFE);
T97 :
array[0.. 3]
of Cardinal = ($74BAE4A0,$9F9867E4,$149D74B8,$FFFFFFFE);
T673:
array[0..21]
of Cardinal = ($C05A8C20,$7A08BB4E,$6545B0DE,$3DCE2A68,
$6722F3FE,$16BB3891,$266C16EA,$2DF1F0D4,
$DDADF754,$9E9DA604,$07484B83,$8196E5E7,
$ABBED6EC,$AC3E3ED0,$5DA0D990,$247375A1,
$FF3D139A,$5951CEF1,$EC368A98,$CB744179,
$10C5680D,$FFFFFFFE);
T257:
array[0.. 8]
of Cardinal = ($191854E8,$81E9ABE2,$6CED7CA6,$E0897EE3,
$1DFA441C,$94FADCDB,$1F565E05,$5CA86261,
$FFFFFFFE);
T241:
array[0.. 7]
of Cardinal = ($94EA6880,$C3985CEE,$B37056F6,$D02DE3E8,
$3B345F1E,$670DBDA8,$5CA5DCE8,$FFFE0459);
T17 = $FFFE5CE8;
T13 = $FFFFE9E4;
T9 = $FFFFFF6C;
T7 = $FFFFFFE8;
T5 = $FFFFFFEC;
P1 = $08A19417;
// 193 * 97 * 17 * 13 * 7 * 5
P2 = $165C5F19;
// 9 * 673 * 257 * 241
{$IFDEF ExpensiveSqrTest}
T41 :
array[0.. 1]
of Cardinal = ($7D4AF8C8,$FFFFFE4C);
T37 :
array[0.. 1]
of Cardinal = ($A1DEE164,$FFFFFFE9);
T61 :
array[0.. 1]
of Cardinal = ($F5A60DC4,$E8EC196B);
T59 :
array[0.. 1]
of Cardinal = ($C1846D44,$FDD49DE7);
T53 :
array[0.. 1]
of Cardinal = ($CCFC512C,$FFED228F);
T47 :
array[0.. 1]
of Cardinal = ($E4D8AC20,$FFFFFBCA);
T43 :
array[0.. 1]
of Cardinal = ($7C5C11AC,$FFFFFCA7);
T83 :
array[0.. 2]
of Cardinal = ($015CE164,$57F4EC8D,$FFFD978C);
T79 :
array[0.. 2]
of Cardinal = ($7902D0C8,$BF618AAE,$FFFFECF4);
T73 :
array[0.. 2]
of Cardinal = ($F472ECA0,$DD38BD86,$FFFFFE14);
T71 :
array[0.. 2]
of Cardinal = ($94E26880,$E9B8D68E,$FFFFFFFE);
T67 :
array[0.. 2]
of Cardinal = ($D81439AC,$A63D7E45,$FFFFFFFC);
T251:
array[0.. 7]
of Cardinal = ($650C4D44,$6BE4DD27,$848471C0,$C110A94F,
$E0D6AF77,$9FC71DED,$91B44D82,$FDD4DCF5);
T239:
array[0.. 7]
of Cardinal = ($14A86080,$8B30CEE8,$D254F662,$48ED8F82,
$D5B4BE0E,$F32EB990,$EAD7E88C,$FFFFFEF9);
T233:
array[0.. 7]
of Cardinal = ($09721C68,$2A61BF8C,$E5DDEA7A,$A4CC948B,
$5EEE9F44,$F6195179,$E13A40C7,$FFFFFE58);
T229:
array[0.. 7]
of cardinal = ($F1E425C4,$885483CD,$37D0A72E,$9DBF6E65,
$3942FB29,$F04A845D,$E909E3EC,$FFFFFFE8);
T227:
array[0.. 7]
of Cardinal = ($8156E164,$35DC66E9,$E949011C,$F8EC8E05,
$77F6D685,$899C453C,$978957E6,$FFFFFFFD);
T223:
array[0.. 6]
of Cardinal = ($0DF03C68,$2A597508,$BDB1A8C8,$2C6699CB,
$ECEA7242,$EF5165AB,$E9C3F04F);
T211:
array[0.. 6]
of Cardinal = ($BCC6958C,$B20507CB,$5F602D18,$9059D546,
$FB2E74BF,$CC22C1F5,$FFFCE569);
T199:
array[0.. 6]
of Cardinal = ($496A9848,$18C116E4,$A1BC7EB8,$81C27A2B,
$977CE7E2,$E6A96DD8,$FFFFFFED);
T197:
array[0.. 6]
of Cardinal = ($C836792C,$07157049,$CAD5EFBC,$7DEAD4CC,
$83AA380F,$279B04E4,$FFFFFFED);
T191:
array[0.. 5]
of Cardinal = ($B0684880,$E7A0966A,$EB9C16C4,$DC97C628,
$A996FA18,$FEEDE9F2);
T181:
array[0.. 5]
of Cardinal = ($D5EE05C4,$A66C8309,$BF7875B4,$994B6B87,
$EAE4304D,$FFE8E81D);
T179:
array[0.. 5]
of Cardinal = ($55A40DC4,$C4E4136F,$5C50C3A0,$8DCFA3CF,
$A550937D,$FFFDC4FD);
T173:
array[0.. 5]
of Cardinal = ($5C1E19AC,$EC257481,$68C59DF6,$A90DDBEE,
$1E0EA04B,$FFFFED66);
T167:
array[0.. 5]
of Cardinal = ($4492A420,$18B86BAC,$9C4DC6F8,$29E2E7E0,
$DAB6DDCA,$FFFFFFFB);
T163:
array[0.. 5]
of Cardinal = ($F89E39AC,$88153421,$5245DB5C,$BD357EEC,
$A6386E07,$FFFFFFFC);
T157:
array[0.. 4]
of Cardinal = ($35F481E4,$F8E42A45,$D9B9E766,$289509C7,$E9E04BEB);
T151:
array[0.. 4]
of Cardinal = ($5D80F0C8,$B379420A,$328CAACE,$45AFBD61,$FFECF0FE);
T149:
array[0.. 4]
of Cardinal = ($08A4FD0C,$5F9D1B45,$7E98EDC6,$4428B62E,$FFEC2FC9);
T139:
array[0.. 4]
of Cardinal = ($0CEED50C,$7D250983,$F5B41F50,$488CF3E6,$FFFFFCF5);
T137:
array[0.. 4]
of Cardinal = ($ADB03468,$46F9EF0A,$DE7D88CC,$B036D543,$FFFFFE58);
T131:
array[0.. 4]
of Cardinal = ($E5CE4544,$034C8521,$B5ECD3FC,$D5D8C587,$FFFFFFFD);
T127:
array[0.. 3]
of Cardinal = ($399054E8,$8FE96982,$BE69680E,$E8D5F663);
T113:
array[0.. 3]
of Cardinal = ($29BA1468,$0CC1EDEE,$7651DEDE,$FFFE58A1);
T109:
array[0.. 3]
of Cardinal = ($418E6D44,$4FFC97A3,$9C60B17A,$FFFFE8AD);
T107:
array[0.. 3]
of Cardinal = ($957681E4,$9CCC6841,$E91567DE,$FFFFFD87);
T103:
array[0.. 3]
of Cardinal = ($89701C68,$4269BDA8,$C7F16EEA,$FFFFFFE9);
T101:
array[0.. 3]
of Cardinal = ($3C049D8C,$FAAD57CD,$6E480F2C,$FFFFFFEC);
T89 :
array[0.. 2]
of Cardinal = ($FD88F0C8,$FD594A6A,$FE4C3C46);
T31 = $EDE2B848;
T29 = $EC2EDD0C;
T23 = $FFFACCA0;
T19 = $FFFCF50C;
T11 = $FFFFFDC4;
T3 = $FFFFFFFC;
I3 = $AAAAAAAB;
// modular inverse of 3 to 2^32, eg 3^-1 mod 2^32
I11 = $BA2E8BA3;
I19 = $286BCA1B;
I23 = $E9BD37A7;
I29 = $4F72C235;
I31 = $BDEF7BDF;
I37 = $914C1BAD;
I41 = $C18F9C19;
I43 = $2FA0BE83;
I47 = $677D46CF;
I53 = $8C13521D;
I59 = $A08AD8F3;
I61 = $C10C9715;
I67 = $07A44C6B;
I71 = $E327A977;
I73 = $C7E3F1F9;
I79 = $613716AF;
I83 = $2B2E43DB;
I89 = $FA3F47E9;
I101 = $7C32B16D;
I103 = $D3431B57;
I107 = $8D28AC43;
I109 = $DA6C0965;
I113 = $0FDBC091;
I127 = $EFDFBF7F;
I131 = $C9484E2B;
I137 = $077975B9;
I139 = $70586723;
I149 = $8CE2CABD;
I151 = $BF937F27;
I157 = $2C0685B5;
I163 = $451AB30B;
I167 = $DB35A717;
I173 = $0D516325;
I179 = $D962AE7B;
I181 = $10F8ED9D;
I191 = $EE936F3F;
I197 = $3D137E0D;
I199 = $EF46C0F7;
I211 = $6E68575B;
I223 = $DB43BB1F;
I227 = $9BA144CB;
I229 = $478BBCED;
I233 = $1FDCD759;
I239 = $437B2E0F;
I251 = $9A020A33;
P3 = $6A917C99;
// 41 * 37 * 31 * 29 * 23 * 19 * 3
P4 = $FCC0D7AD;
// 61 * 59 * 53 * 47 * 43 * 11
P5 = $87B81DA9;
// 83 * 79 * 73 * 71 * 67
P6 = $BEC8D631;
// 251 * 239 * 233 * 229
P7 = $7EB0F0B1;
// 227 * 223 * 211 * 199
P8 = $48A9A435;
// 197 * 191 * 181 * 179
P9 = $2C11944D;
// 173 * 167 * 163 * 157
P10 = $19899AC9;
// 151 * 149 * 139 * 137
P11 = $0C36CCA9;
// 131 * 127 * 113 * 109
P12 = $05E7A779;
// 89 * 101 * 103 * 107
IP3 = $10BB17A9;
// P3^-1 mod 2^32
IP4 = $FC8EA425;
// P4^-1 mod 2^32
IP5 = $75B3D699;
// P5^-1 mod 2^32
IP6 = $0138C2D1;
// P6^-1 mod 2^32
IP7 = $C5775851;
// P7^-1 mod 2^32
IP8 = $BD478E1D;
// P8^-1 mod 2^32
IP9 = $701FC485;
// P9^-1 mod 2^32
IP10 = $36BB9F79;
// P10^-1 mod 2^32
IP11 = $E9489799;
// P11^-1 mod 2^32
IP12 = $4D4F12C9;
// P12^-1 mod 2^32
{$ENDIF}
{prime count probability cum. prop.
256 44 0.17187500 0.171875
0.17187500000000000
193 97 0.50259067 0.08638277202072538860103626943
97 49 0.50515464 0.043636658031088082901554404145
17 9 0.52941176 0.023101760134105455653764096312
13 7 0.53846154 0.012439409302979860736642205707
7 4 0.57142857 0.007108233887417063278081260404
5 3 0.60000000 0.004264940332450237966848756242
0.02481419829789229
9 4 0.44444444 0.001895529036644550207488336108
673 337 0.50074294 0.000949172786551580118757160874
257 129 0.50194553 0.000476433032938341771671882306
241 121 0.50207469 0.000239204966744976574158911863
0.05608635715837826
41 21 0.51219512 0.00012251961711328068432529632
37 19 0.51351351 0.00006291547905817116222109811
31 16 0.51612903 0.000032472505320346406307663541
29 15 0.51724138 0.00001679612344155848602120528
23 12 0.52173913 0.000008763194839073992706715798
19 10 0.52631579 0.00000461220781003894352985042
3 2 0.66666667 0.00000307480520669262901990028
0.01285426991142190
61 31 0.50819672 0.000001562605924712647534703421
59 30 0.50847458 0.000000794545385447108915950892
53 27 0.50943396 0.00000040476840390701774963536
47 24 0.51063830 0.000000206690248803583531728695
43 22 0.51162791 0.000000105748499387879946465844
11 6 0.54545455 0.000000057680999666116334435915
0.01875923702111852
83 42 0.50602410 0.000000029187975734661277666366
79 40 0.50632911 0.000000014778721890967735527274
73 37 0.50684932 0.000000007490585068024742664509
71 36 0.50704225 0.000000003798043133082968111582
67 34 0.50746269 0.000000001927365172012252474534
0.03341421236054701
251 126 0.50199203 0.000000000967521958858740286021
239 120 0.50209205 0.00000000048578508394581102227
233 117 0.50214592 0.000000000243934999234591800882
229 115 0.50218341 0.000000000122500108785930380356
0.06355832852268207
227 114 0.50220264 0.000000000061519878421128032425
223 112 0.50224215 0.000000000030897876157696590276
211 106 0.50236967 0.000000000015522155794861794167
199 100 0.50251256 0.000000000007800078288875273451
0.06367405193497382
197 99 0.50253807 0.000000000003919836297455086658
191 96 0.50261780 0.00000000000197017950029156188
181 91 0.50276243 0.000000000000990532234953216194
179 90 0.50279330 0.00000000000049803296729491317
0.06384973956033545
173 87 0.50289017 0.000000000000250455885287037259
167 84 0.50299401 0.000000000000125977810563539699
163 82 0.50306748 0.000000000000063375340283498499
157 79 0.50318471 0.000000000000031889502435645742
0.06403090664631079
151 76 0.50331126 0.00000000000001605034559674885
149 75 0.50335570 0.000000000000008079033018497743
139 70 0.50359712 0.000000000000004068577779099583
137 69 0.50364964 0.000000000000002049137713561104
0.06425743762218754
131 66 0.50381679 0.000000000000001032389993091854
127 64 0.50393701 0.000000000000000520259524077785
113 57 0.50442478 0.000000000000000262431795331272
109 55 0.50458716 0.00000000000000013241971324055
0.06462216392983359
107 54 0.50467290 0.000000000000000066828640327007
103 52 0.50485437 0.000000000000000033738731038877
101 51 0.50495050 0.000000000000000017036388940423
89 45 0.50561798 0.000000000000000008613904520439
0.06505001641855890
1036 Bytes needed for all Lookuptable
252 Bytes needed for no expensive test
Follow code examine in avg. ~2.5 Cycles per Digit if A is NOT a perfect Sqr.
on 2048 Bit A ~ 3000 times faster as NSqrt() }
var
R,T: Cardinal;
S: Long96;
begin
Result := (A =
nil)
or (PInt(A).FCount = 0);
if not Result
and not PInt(A).FNeg
and NCheckSqrTable(A[0])
then
begin
// prop: 0.171875 pass above
// ~12 Cycles
S := NSum(A);
// S = A mod 2^96-1 // A.Count x Additions + 1 DIV
R := NMod(S, P1);
T := R
mod 193;
if Odd(T193[T
shr 5]
shr (T
and 31))
then Exit;
T := R
mod 97;
if Odd(T97 [T
shr 5]
shr (T
and 31))
then Exit;
if Odd(T17
shr (R
mod 17))
then Exit;
if Odd(T13
shr (R
mod 13))
then Exit;
if Odd(T7
shr (R
mod 7))
then Exit;
if Odd(T5
shr (R
mod 5))
then Exit;
NCalcCheck;
// ~ 4 Cycles
R := NMod(S, P2);
// 3 x 32Bit DIV's
if Odd(T9
shr (R
mod 9))
then Exit;
T := R
mod 673;
if Odd(T673[T
shr 5]
shr (T
and 31))
then Exit;
T := R
mod 257;
if Odd(T257[T
shr 5]
shr (T
and 31))
then Exit;
T := R
mod 241;
if Odd(T241[T
shr 5]
shr (T
and 31))
then Exit;
// prop: 0.0002392 pass above, 1 of 4181 Candidates pass
{$IFDEF ExpensiveSqrTest}
R := DModDInv2k32(PInt(A).FNum, P3, PInt(A).FCount, IP3);
T := DInvMul2k32(R, 41, I41);
if Odd(T41[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 37, I37);
if Odd(T37[T
shr 5]
shr (T
and 31))
then Exit;
if Odd(T31
shr DInvMul2k32(R, 31, I31))
then Exit;
if Odd(T29
shr DInvMul2k32(R, 29, I29))
then Exit;
if Odd(T23
shr DInvMul2k32(R, 23, I23))
then Exit;
if Odd(T19
shr DInvMul2k32(R, 19, I19))
then Exit;
if Odd(T3
shr DInvMul2k32(R, 3, I3))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P4, PInt(A).FCount, IP4);
T := DInvMul2k32(R, 61, I61);
if Odd(T61[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 59, I59);
if Odd(T59[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 53, I53);
if Odd(T53[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 47, I47);
if Odd(T47[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 43, I43);
if Odd(T43[T
shr 5]
shr (T
and 31))
then Exit;
if Odd(T11
shr DInvMul2k32(R, 11, I11))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P5, PInt(A).FCount, IP5);
T := DInvMul2k32(R, 83, I83);
if Odd(T83[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 79, I79);
if Odd(T79[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 73, I73);
if Odd(T73[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 71, I71);
if Odd(T71[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 67, I67);
if Odd(T67[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P6, PInt(A).FCount, IP6);
T := DInvMul2k32(R, 251, I251);
if Odd(T251[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 239, I239);
if Odd(T239[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 233, I233);
if Odd(T233[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 229, I229);
if Odd(T229[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P7, PInt(A).FCount, IP7);
T := DInvMul2k32(R, 227, I227);
if Odd(T227[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 223, I223);
if Odd(T223[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 211, I211);
if Odd(T211[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 199, I199);
if Odd(T199[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P8, PInt(A).FCount, IP8);
T := DInvMul2k32(R, 197, I197);
if Odd(T197[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 191, I191);
if Odd(T191[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 181, I181);
if Odd(T181[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 179, I179);
if Odd(T179[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P9, PInt(A).FCount, IP9);
T := DInvMul2k32(R, 173, I173);
if Odd(T173[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 167, I167);
if Odd(T167[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 163, I163);
if Odd(T163[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 157, I157);
if Odd(T157[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P10, PInt(A).FCount, IP10);
T := DInvMul2k32(R, 151, I151);
if Odd(T151[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 149, I149);
if Odd(T149[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 139, I139);
if Odd(T139[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 137, I137);
if Odd(T137[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P11, PInt(A).FCount, IP11);
T := DInvMul2k32(R, 131, I131);
if Odd(T131[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 127, I127);
if Odd(T127[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 113, I113);
if Odd(T113[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 109, I109);
if Odd(T109[T
shr 5]
shr (T
and 31))
then Exit;
R := DModDInv2k32(PInt(A).FNum, P12, PInt(A).FCount, IP12);
T := DInvMul2k32(R, 107, I107);
if Odd(T107[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 103, I103);
if Odd(T103[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 101, I101);
if Odd(T101[T
shr 5]
shr (T
and 31))
then Exit;
T := DInvMul2k32(R, 89, I89 );
if Odd(T89 [T
shr 5]
shr (T
and 31))
then Exit;
// 1 of 116091372690195275 candidates come here
{$ENDIF}
Result := True;
end;
end;
function NCheckSqrt(
const A: IInteger): Boolean;
asm // check if A is a perfect Square, by doing a NSqrt(nil, A)
MOV EDX,EAX
XOR EAX,EAX
JMP NSqrt
// Result := NSqrt(nil, A);
end;
function NIsPerfectSqr(
const A: IInteger; FullTest: Boolean): Boolean;
// test if A is a perfect square
begin
Result := NCheckSqr(A)
and (
not FullTest
or NCheckSqrt(A));
end;
function NIsPerfectSqr(
const N: Int64): Boolean;
var
R: Cardinal;
begin
Result := N = 0;
if not Result
and (N > 0)
and NCheckSqrTable(N)
then
Result := (DSqrt64(@R, @R, @N, 2) = 0)
and (R = 0);
end;
function NIsPerfectPower(
var B: IInteger;
const N: IInteger; Bound: Cardinal = 0): Cardinal;
// check if N is a perfect power, if true returns Base in B and Exponent (smallprime) in Result
// implementation of
// "Detecting perfect powers in essentially linear time"
// Daniel J. Bernstein, file: powers.ps
var
X,Y: IInteger;
I,K,L: Integer;
P: TSmallPrimeSieve;
label
Error;
begin
Result := 0;
if Bound = 0
then Bound := TSmallPrimeSieve.MaxPrime;
I := NSgn(N, True);
if I < 0
then goto Error;
if I < 3
then
begin
if @B <>
nil then NSet(B, I);
Result := 1;
Exit;
end;
if NCheckSqr(N)
then
if @B <>
nil then
begin
if NSqrt(X, N)
then
begin
NSwp(B, X);
Result := 2;
Exit;
end;
end else
if NCheckSqrt(N)
then
begin
Result := 2;
Exit;
end;
if not NOdd(N)
then
begin
// to slow,
{ I := 2;
L := Primes[I];
while L < K do
begin
if NRoot(X, N, L) then
begin
if @B <> nil then NSwp(B, X);
Result := L;
Exit;
end;
if L >= Bound then goto Error;
Inc(I);
L := Primes[I];
end; }
goto Error;
end else
if Bound > 2
then
begin
K := NSize(N);
L := (K +1)
shr 1;
if NInvMod2k(Y, N, L +1)
then // Y = N^-1 mod 2^(Round(K/2) +1)
begin
// first check to exponent 2,
// suppressed: NIsPerfectSqr() is always faster
{ if NSqrtMod2k_2adic(X, Y, L) then // X = Y^1/2 mod 2^(Round(K/2)
begin
if NIsPower(N, X, 2) then
begin
if @B <> nil then NSwp(B, X);
Result := 2;
Exit;
end;
NCpl(X, L); // X = 2^Round(K/2) - X
if NIsPower(N, X, 2) then
begin
if @B <> nil then NSwp(B, X);
Result := 2;
Exit;
end;
end;}
// now check to smallprime powers
P := Primes;
NAutoRelease(P);
I := 2;
L := P[I];
while L < K
do
begin
if NRootMod2k_2adic(X, Y, NInt(L), (K + L -1)
div L)
and NIsPower(N, X, L)
then
begin
if @B <>
nil then NSwp(B, X);
Result := L;
Exit;
end;
if Cardinal(L) >= Bound
then goto Error;
Inc(I);
L := P[I];
end;
if @B <>
nil then NSet(B, N);
Result := 1;
end else
begin
Error:
if @B <>
nil then NSet(B, 0);
Result := 0;
end;
end;
end;