Factoring a Number

Unit FGInt.pas

Unit FGInt;

Interface

Uses SysUtils, Math;

Type
   TCompare = (Lt, St, Eq, Er);
   TSign = (negative, positive);
   TFGInt = Record
      Sign : TSign;
      Number : Array Of LongWord;
   End;

Procedure FGIntToBase2String(Const FGInt : TFGInt; Var S : String);
Procedure Base2StringToFGInt(S : String; Var FGInt : TFGInt);
Procedure Base10StringToFGInt(Base10 : String; Var FGInt : TFGInt);
Procedure FGIntToBase10String(Const FGInt : TFGInt; Var Base10 : String);
Procedure FGIntDestroy(Var FGInt : TFGInt);
Function FGIntCompareAbs(Const FGInt1, FGInt2 : TFGInt) : TCompare;
Procedure FGIntAdd(Const FGInt1, FGInt2 : TFGInt; Var Sum : TFGInt);
Procedure FGIntChangeSign(Var FGInt : TFGInt);
Procedure FGIntSub(Var FGInt1, FGInt2, dif : TFGInt);
Procedure FGIntMulByIntbis(Var FGInt : TFGInt; by : LongWord);
Procedure FGIntDivByIntBis(Var FGInt : TFGInt; by : LongWord; Var modres : LongWord);
Procedure FGIntModByInt(Const FGInt : TFGInt; by : LongWord; Var modres : LongWord);
Procedure FGIntAbs(Var FGInt : TFGInt);
Procedure FGIntCopy(Const FGInt1 : TFGInt; Var FGInt2 : TFGInt);
Procedure FGIntShiftRight(Var FGInt : TFGInt);
Procedure FGIntShiftRightBy31(Var FGInt : TFGInt);
Procedure FGIntSubBis(Var FGInt1 : TFGInt; Const FGInt2 : TFGInt);
Procedure FGIntMul(Const FGInt1, FGInt2 : TFGInt; Var Prod : TFGInt);
Procedure FGIntSquare(Const FGInt : TFGInt; Var Square : TFGInt);
Procedure FGIntShiftLeftBy31(Var FGInt : TFGInt);
Procedure FGIntDivMod(Var FGInt1, FGInt2, QFGInt, MFGInt : TFGInt);
Procedure FGIntDiv(Var FGInt1, FGInt2, QFGInt : TFGInt);
Procedure FGIntMod(Var FGInt1, FGInt2, MFGInt : TFGInt);
Procedure FGIntSquareMod(Var FGInt, Modb, FGIntSM : TFGInt);
Procedure FGIntMulMod(Var FGInt1, FGInt2, base, FGIntres : TFGInt);
Procedure FGIntModBis(Const FGInt : TFGInt; Var FGIntOut : TFGInt; b, head : LongWord);
Procedure FGIntMulModBis(Const FGInt1, FGInt2 : TFGInt; Var Prod : TFGInt; b, head : LongWord);
Procedure FGIntMontgomeryMod(Const GInt, base, baseInv : TFGInt; Var MGInt : TFGInt; b : Longword; head : LongWord);
Procedure FGIntMontgomeryModExp(Var FGInt, exp, modb, res : TFGInt);
Procedure FGIntGCD(Const FGInt1, FGInt2 : TFGInt; Var GCD : TFGInt);
Procedure FGIntRabinMiller(Var FGIntp : TFGInt; nrtest : Longword; Var ok : boolean);
Procedure FGIntModInv(Const FGInt1, base : TFGInt; Var Inverse : TFGInt);

Implementation

Var
   primes : Array[1..1228] Of integer =
      (3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
      131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
      257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
      397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
      547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
      683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
      853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009,
      1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123,
      1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
      1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429,
      1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553,
      1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
      1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847,
      1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
      1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131,
      2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
      2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,
      2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593,
      2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719,
      2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
      2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037,
      3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
      3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359,
      3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527,
      3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659,
      3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
      3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967,
      3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
      4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273,
      4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
      4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637,
      4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789,
      4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957,
      4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101,
      5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281,
      5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
      5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623,
      5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779,
      5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903,
      5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101,
      6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269,
      6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397,
      6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599,
      6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779,
      6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947,
      6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
      7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283,
      7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487,
      7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607,
      7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
      7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951,
      7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161,
      8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, 8293, 8297, 8311,
      8317, 8329, 8353, 8363, 8369, 8377, 8387, 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, 8513, 8521,
      8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, 8681,
      8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831,
      8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007,
      9011, 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181,
      9187, 9199, 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, 9293, 9311, 9319, 9323, 9337, 9341, 9343,
      9349, 9371, 9377, 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, 9461, 9463, 9467, 9473, 9479, 9491,
      9497, 9511, 9521, 9533, 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, 9643, 9649, 9661, 9677, 9679,
      9689, 9697, 9719, 9721, 9733, 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, 9817, 9829, 9833, 9839,
      9851, 9857, 9859, 9871, 9883, 9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973);

// Convert a FGInt to a binary string (base 2) & visa versa
Procedure FGIntToBase2String(Const FGInt : TFGInt; Var S : String);
Var
   i : LongWord;
   j : integer;
Begin
   S := '';
   For i := 1 To FGInt.Number[0] Do
   Begin
      For j := 0 To 30 Do
         If (1 And (FGInt.Number[i] Shr j)) = 1 Then
            S := '1' + S
         Else
            S := '0' + S;
   End;
   While (length(S) > 1) And (S[1] = '0') Do
      delete(S, 1, 1);
   If S = '' Then S := '0';
End;

Procedure Base2StringToFGInt(S : String; Var FGInt : TFGInt);
Var
   i, j, size : LongWord;
Begin
   While (S[1] = '0') And (length(S) > 1) Do
      delete(S, 1, 1);
   size := length(S) Div 31;
   If (length(S) Mod 31) <> 0 Then size := size + 1;
   SetLength(FGInt.Number, (size + 1));
   FGInt.Number[0] := size;
   j := 1;
   FGInt.Number[j] := 0;
   i := 0;
   While length(S) > 0 Do
   Begin
      If S[length(S)] = '1' Then
         FGInt.Number[j] := FGInt.Number[j] Or (1 Shl i);
      i := i + 1;
      If i = 31 Then
      Begin
         i := 0;
         j := j + 1;
         If j <= size Then FGInt.Number[j] := 0;
      End;
      delete(S, length(S), 1);
   End;
   FGInt.Sign := positive;
End;

// Convert a base 10 string to a FGInt
Procedure Base10StringToFGInt(Base10 : String; Var FGInt : TFGInt);
Var
   i, size : LongWord;
   j : word;
   S, x : String;
   sign : TSign;

   Procedure GIntDivByIntBis1(Var GInt : TFGInt; by : LongWord; Var modres : word);
   Var
      i, size, rest, temp : LongWord;
   Begin
      size := GInt.Number[0];
      temp := 0;
      For i := size Downto 1 Do
      Begin
         temp := temp * 10000;
         rest := temp + GInt.Number[i];
         GInt.Number[i] := rest Div by;
         temp := rest Mod by;
      End;
	modres := temp;
      While (GInt.Number[size] = 0) And (size > 1) Do
         size := size - 1;
      If size <> GInt.Number[0] Then
      Begin
         SetLength(GInt.Number, size + 1);
         GInt.Number[0] := size;
      End;
   End;

Begin
   While (length(Base10) > 1) and (Not (Base10[1] In ['-', '0'..'9'])) Do
      delete(Base10, 1, 1);
   If copy(Base10, 1, 1) = '-' Then
   Begin
      Sign := negative;
      delete(Base10, 1, 1);
   End
   Else
      Sign := positive;
   While (length(Base10) > 1) And (copy(Base10, 1, 1) = '0') Do
      delete(Base10, 1, 1);
   size := length(Base10) Div 4;
   If (length(Base10) Mod 4) <> 0 Then size := size + 1;
   SetLength(FGInt.Number, size + 1);
   FGInt.Number[0] := size;
   For i := 1 To (size - 1) Do
   Begin
      x := copy(Base10, length(Base10) - 3, 4);
      FGInt.Number[i] := StrToInt(x);
      delete(Base10, length(Base10) - 3, 4);
   End;
   FGInt.Number[size] := StrToInt(Base10);

   S := '';
   While (FGInt.Number[0] <> 1) Or (FGInt.Number[1] <> 0) Do
   Begin
      GIntDivByIntBis1(FGInt, 2, j);
      S := inttostr(j) + S;
   End;
   If S = '' Then S := '0';
   FGIntDestroy(FGInt);
   Base2StringToFGInt(S, FGInt);
   FGInt.Sign := sign;
End;

// Convert a FGInt to a base 10 string
Procedure FGIntToBase10String(Const FGInt : TFGInt; Var Base10 : String);
Var
   S : String;
   j : LongWord;
   temp : TFGInt;
Begin
   FGIntCopy(FGInt, temp);
   Base10 := '';
   While (temp.Number[0] > 1) Or (temp.Number[1] > 0) Do
   Begin
      FGIntDivByIntBis(temp, 10000, j);
      S := IntToStr(j);
      While Length(S) < 4 Do
         S := '0' + S;
      Base10 := S + Base10;
   End;
   Base10 := '0' + Base10;
   While (length(Base10) > 1) And (Base10[1] = '0') Do
      delete(Base10, 1, 1);
   If FGInt.Sign = negative Then Base10 := '-' + Base10;
End;


// Destroy a FGInt to free memory
Procedure FGIntDestroy(Var FGInt : TFGInt);
Begin
   FGInt.Number := Nil;
End;

// Compare 2 FGInts in absolute value, returns
// Lt if FGInt1 > FGInt2, St if FGInt1 < FGInt2, Eq if FGInt1 = FGInt2,
// Er otherwise
Function FGIntCompareAbs(Const FGInt1, FGInt2 : TFGInt) : TCompare;
Var
   size1, size2, i : LongWord;
Begin
   FGIntCompareAbs := Er;
   size1 := FGInt1.Number[0];
   size2 := FGInt2.Number[0];
   If size1 > size2 Then FGIntCompareAbs := Lt Else
      If size1 < size2 Then FGIntCompareAbs := St Else
      Begin
         i := size2;
         While (FGInt1.Number[i] = FGInt2.Number[i]) And (i > 1) Do i := i - 1;
         If FGInt1.Number[i] = FGInt2.Number[i] Then FGIntCompareAbs := Eq Else
            If FGInt1.Number[i] < FGInt2.Number[i] Then FGIntCompareAbs := St Else
               If FGInt1.Number[i] > FGInt2.Number[i] Then FGIntCompareAbs := Lt;
      End;
End;

// Add 2 FGInts, FGInt1 + FGInt2 = Sum
Procedure FGIntAdd(Const FGInt1, FGInt2 : TFGInt; Var Sum : TFGInt);
Var
   i, size1, size2, size, rest, Trest : LongWord;
Begin
   size1 := FGInt1.Number[0];
   size2 := FGInt2.Number[0];
   If size1 < size2 Then
      FGIntAdd(FGInt2, FGInt1, Sum)
   Else
   Begin
      If FGInt1.Sign = FGInt2.Sign Then
      Begin
         Sum.Sign := FGInt1.Sign;
         setlength(Sum.Number, (size1 + 2));
         rest := 0;
         For i := 1 To size2 Do
         Begin
            Trest := FGInt1.Number[i];
            Trest := Trest + FGInt2.Number[i];
            Trest := Trest + rest;
            Sum.Number[i] := Trest And 2147483647;
            rest := Trest Shr 31;
         End;
         For i := (size2 + 1) To size1 Do
         Begin
            Trest := FGInt1.Number[i] + rest;
            Sum.Number[i] := Trest And 2147483647;
            rest := Trest Shr 31;
         End;
         size := size1 + 1;
         Sum.Number[0] := size;
         Sum.Number[size] := rest;
         While (Sum.Number[size] = 0) And (size > 1) Do
            size := size - 1;
         If Sum.Number[0] <> size Then SetLength(Sum.Number, size + 1);
         Sum.Number[0] := size;
      End
      Else
      Begin
         If FGIntCompareAbs(FGInt2, FGInt1) = Lt Then
            FGIntAdd(FGInt2, FGInt1, Sum)
         Else
         Begin
            SetLength(Sum.Number, (size1 + 1));
            rest := 0;
            For i := 1 To size2 Do
            Begin
               Trest := 2147483648;
               TRest := Trest + FGInt1.Number[i];
               TRest := Trest - FGInt2.Number[i];
               TRest := Trest - rest;
               Sum.Number[i] := Trest And 2147483647;
               If (Trest > 2147483647) Then
                  rest := 0
               Else
                  rest := 1;
            End;
            For i := (size2 + 1) To size1 Do
            Begin
               Trest := 2147483648;
               TRest := Trest + FGInt1.Number[i];
               TRest := Trest - rest;
               Sum.Number[i] := Trest And 2147483647;
               If (Trest > 2147483647) Then
                  rest := 0
               Else
                  rest := 1;
            End;
            size := size1;
            While (Sum.Number[size] = 0) And (size > 1) Do
               size := size - 1;
            If size <> size1 Then SetLength(Sum.Number, size + 1);
            Sum.Number[0] := size;
            Sum.Sign := FGInt1.Sign;
         End;
      End;
   End;
End;


Procedure FGIntChangeSign(Var FGInt : TFGInt);
Begin
   If FGInt.Sign = negative Then FGInt.Sign := positive Else FGInt.Sign := negative;
End;

// Substract 2 FGInts, FGInt1 - FGInt2 = dif
Procedure FGIntSub(Var FGInt1, FGInt2, dif : TFGInt);
Begin
   FGIntChangeSign(FGInt2);
   FGIntAdd(FGInt1, FGInt2, dif);
   FGIntChangeSign(FGInt2);
End;

// multiply a FGInt by an integer, FGInt * by = res, by < 1000000000
Procedure FGIntMulByIntbis(Var FGInt : TFGInt; by : LongWord);
Var
   i, size, rest : LongWord;
   Trest : int64;
Begin
   size := FGInt.Number[0];
   Setlength(FGInt.Number, size + 2);
   rest := 0;
   For i := 1 To size Do
   Begin
      Trest := FGInt.Number[i];
      TRest := Trest * by;
      TRest := Trest + rest;
      FGInt.Number[i] := Trest And 2147483647;
      rest := Trest Shr 31;
   End;
   If rest <> 0 Then
   Begin
      size := size + 1;
      FGInt.Number[size] := rest;
   End
   Else
      SetLength(FGInt.Number, size + 1);
   FGInt.Number[0] := size;
End;

// divide a FGInt by an integer, FGInt = FGInt * by + modres
Procedure FGIntDivByIntBis(Var FGInt : TFGInt; by : LongWord; Var modres : LongWord);
Var
   i, size : LongWord;
   temp, rest : int64;
Begin
   size := FGInt.Number[0];
   temp := 0;
   For i := size Downto 1 Do
   Begin
      temp := temp Shl 31;
      rest := temp Or FGInt.Number[i];
      FGInt.Number[i] := rest Div by;
      temp := rest Mod by;
   End;
   modres := temp;
   While (FGInt.Number[size] = 0) And (size > 1) Do
      size := size - 1;
   If size <> FGInt.Number[0] Then
   Begin
      SetLength(FGInt.Number, size + 1);
      FGInt.Number[0] := size;
   End;
End;


// Reduce a FGInt modulo by (=an integer), FGInt mod by = modres
Procedure FGIntModByInt(Const FGInt : TFGInt; by : LongWord; Var modres : LongWord);
Var
   i, size : LongWord;
   temp, rest : int64;
Begin
   size := FGInt.Number[0];
   temp := 0;
   For i := size Downto 1 Do
   Begin
      temp := temp Shl 31;
      rest := temp Or FGInt.Number[i];
      temp := rest Mod by;
   End;
   modres := temp;
   If FGInt.sign = negative Then modres := by - modres;
End;


// Returns the FGInt in absolute value
Procedure FGIntAbs(Var FGInt : TFGInt);
Begin
   FGInt.Sign := positive;
End;


// Copy a FGInt1 into FGInt2
Procedure FGIntCopy(Const FGInt1 : TFGInt; Var FGInt2 : TFGInt);
Begin
   FGInt2.Sign := FGInt1.Sign;
   FGInt2.Number := Nil;
   FGInt2.Number := Copy(FGInt1.Number, 0, FGInt1.Number[0] + 1);
End;


// Shift the FGInt to the right in base 2 notation, ie FGInt = FGInt div 2
Procedure FGIntShiftRight(Var FGInt : TFGInt);
Var
   l, m, i, size : LongWord;
Begin
   size := FGInt.Number[0];
   l := 0;
   For i := size Downto 1 Do
   Begin
      m := FGInt.Number[i] And 1;
      FGInt.Number[i] := (FGInt.Number[i] Shr 1) Or l;
      l := m Shl 30;
   End;
   If (FGInt.Number[size] = 0) And (size > 1) Then
   Begin
      setlength(FGInt.Number, size);
      FGInt.Number[0] := size - 1;
   End;
End;


// FGInt = FGInt / 2147483648
Procedure FGIntShiftRightBy31(Var FGInt : TFGInt);
Var
   size, i : LongWord;
Begin
   size := FGInt.Number[0];
   If size > 1 Then
   Begin
      For i := 1 To size - 1 Do
      Begin
         FGInt.Number[i] := FGInt.Number[i + 1];
      End;
      SetLength(FGInt.Number, Size);
      FGInt.Number[0] := size - 1;
   End
   Else
      FGInt.Number[1] := 0;
End;


// FGInt1 = FGInt1 - FGInt2, use only when 0 < FGInt2 < FGInt1
Procedure FGIntSubBis(Var FGInt1 : TFGInt; Const FGInt2 : TFGInt);
Var
   i, size1, size2, rest, Trest : LongWord;
Begin
   size1 := FGInt1.Number[0];
   size2 := FGInt2.Number[0];
   rest := 0;
   For i := 1 To size2 Do
   Begin
      Trest := (2147483648 Or FGInt1.Number[i]) - FGInt2.Number[i] - rest;
      If (Trest > 2147483647) Then
         rest := 0
      Else
         rest := 1;
      FGInt1.Number[i] := Trest And 2147483647;
   End;
   For i := size2 + 1 To size1 Do
   Begin
      Trest := (2147483648 Or FGInt1.Number[i]) - rest;
      If (Trest > 2147483647) Then
         rest := 0
      Else
         rest := 1;
      FGInt1.Number[i] := Trest And 2147483647;
   End;
   i := size1;
   While (FGInt1.Number[i] = 0) And (i > 1) Do
      i := i - 1;
   If i <> size1 Then
   Begin
      SetLength(FGInt1.Number, i + 1);
      FGInt1.Number[0] := i;
   End;
End;


// Multiply 2 FGInts, FGInt1 * FGInt2 = Prod
Procedure FGIntMul(Const FGInt1, FGInt2 : TFGInt; Var Prod : TFGInt);
Var
   i, j, size, size1, size2, rest : LongWord;
   Trest : int64;
Begin
   size1 := FGInt1.Number[0];
   size2 := FGInt2.Number[0];
   size := size1 + size2;
   SetLength(Prod.Number, (size + 1));
   For i := 1 To size Do
      Prod.Number[i] := 0;

   For i := 1 To size2 Do
   Begin
      rest := 0;
      For j := 1 To size1 Do
      Begin
         Trest := FGInt1.Number[j];
         Trest := Trest * FGInt2.Number[i];
         Trest := Trest + Prod.Number[j + i - 1];
         Trest := Trest + rest;
         Prod.Number[j + i - 1] := Trest And 2147483647;
         rest := Trest Shr 31;
      End;
      Prod.Number[i + size1] := rest;
   End;

   Prod.Number[0] := size;
   While (Prod.Number[size] = 0) And (size > 1) Do
      size := size - 1;
   If size <> Prod.Number[0] Then
   Begin
      SetLength(Prod.Number, size + 1);
      Prod.Number[0] := size;
   End;
   If FGInt1.Sign = FGInt2.Sign Then
      Prod.Sign := Positive
   Else
      prod.Sign := negative;
End;


// Square a FGInt, FGInt² = Square
Procedure FGIntSquare(Const FGInt : TFGInt; Var Square : TFGInt);
Var
   size, size1, i, j, rest : LongWord;
   Trest : int64;
Begin
   size1 := FGInt.Number[0];
   size := 2 * size1;
   SetLength(Square.Number, (size + 1));
   Square.Number[0] := size;
   For i := 1 To size Do
      Square.Number[i] := 0;
   For i := 1 To size1 Do
   Begin
      Trest := FGInt.Number[i];
      Trest := Trest * FGInt.Number[i];
      Trest := Trest + Square.Number[2 * i - 1];
      Square.Number[2 * i - 1] := Trest And 2147483647;
      rest := Trest Shr 31;
      For j := i + 1 To size1 Do
      Begin
         Trest := FGInt.Number[i] Shl 1;
         Trest := Trest * FGInt.Number[j];
         Trest := Trest + Square.Number[i + j - 1];
         Trest := Trest + rest;
         Square.Number[i + j - 1] := Trest And 2147483647;
         rest := Trest Shr 31;
      End;
      Square.Number[i + size1] := rest;
   End;
   Square.Sign := positive;
   While (Square.Number[size] = 0) And (size > 1) Do
      size := size - 1;
   If size <> (2 * size1) Then
   Begin
      SetLength(Square.Number, size + 1);
      Square.Number[0] := size;
   End;
End;

// FGInt = FGInt * 2147483648
Procedure FGIntShiftLeftBy31(Var FGInt : TFGInt);
Var
   f1, f2 : LongWord;
   i, size : longint;
Begin
   size := FGInt.Number[0];
   SetLength(FGInt.Number, size + 2);
   f1 := 0;
   For i := 1 To (size + 1) Do
   Begin
      f2 := FGInt.Number[i];
      FGInt.Number[i] := f1;
      f1 := f2;
   End;
   FGInt.Number[0] := size + 1;
End;


// Divide 2 FGInts, FGInt1 = FGInt2 * QFGInt + MFGInt, MFGInt is always positive
Procedure FGIntDivMod(Var FGInt1, FGInt2, QFGInt, MFGInt : TFGInt);
Var
   one, zero, temp1, temp2 : TFGInt;
   s1, s2 : TSign;
   j, s, t : LongWord;
   i : int64;
Begin
   s1 := FGInt1.Sign;
   s2 := FGInt2.Sign;
   FGIntAbs(FGInt1);
   FGIntAbs(FGInt2);
   FGIntCopy(FGInt1, MFGInt);
   FGIntCopy(FGInt2, temp1);

   If FGIntCompareAbs(FGInt1, FGInt2) <> St Then
   Begin
      s := FGInt1.Number[0] - FGInt2.Number[0];
      SetLength(QFGInt.Number, (s + 2));
      QFGInt.Number[0] := s + 1;
      For t := 1 To s Do
      Begin
         FGIntShiftLeftBy31(temp1);
         QFGInt.Number[t] := 0;
      End;
      j := s + 1;
      QFGInt.Number[j] := 0;
      While FGIntCompareAbs(MFGInt, FGInt2) <> St Do
      Begin
         While FGIntCompareAbs(MFGInt, temp1) <> St Do
         Begin
            If MFGInt.Number[0] > temp1.Number[0] Then
            Begin
               i := MFGInt.Number[MFGInt.Number[0]];
               i := i Shl 31;
               i := i + MFGInt.Number[MFGInt.Number[0] - 1];
               i := i Div (temp1.Number[temp1.Number[0]] + 1);
            End
            Else
               i := MFGInt.Number[MFGInt.Number[0]] Div (temp1.Number[temp1.Number[0]] + 1);
            If (i <> 0) Then
            Begin
               FGIntCopy(temp1, temp2);
               FGIntMulByIntBis(temp2, i);
               FGIntSubBis(MFGInt, temp2);
               QFGInt.Number[j] := QFGInt.Number[j] + i;
               If FGIntCompareAbs(MFGInt, temp2) <> St Then
               Begin
                  QFGInt.Number[j] := QFGInt.Number[j] + i;
                  FGIntSubBis(MFGInt, temp2);
               End;
               FGIntDestroy(temp2);
            End
            Else
            Begin
               QFGInt.Number[j] := QFGInt.Number[j] + 1;
               FGIntSubBis(MFGInt, temp1);
            End;
         End;
         If MFGInt.Number[0] <= temp1.Number[0] Then
            If FGIntCompareAbs(temp1, FGInt2) <> Eq Then
            Begin
               FGIntShiftRightBy31(temp1);
               j := j - 1;
            End;
      End;
   End
   Else
      Base10StringToFGInt('0', QFGInt);
   s := QFGInt.Number[0];
   While (s > 1) And (QFGInt.Number[s] = 0) Do
      s := s - 1;
   If s < QFGInt.Number[0] Then
   Begin
      setlength(QFGInt.Number, s + 1);
      QFGInt.Number[0] := s;
   End;
   QFGInt.Sign := positive;

   FGIntDestroy(temp1);
   Base10StringToFGInt('0', zero);
   Base10StringToFGInt('1', one);
   If s1 = negative Then
   Begin
      If FGIntCompareAbs(MFGInt, zero) <> Eq Then
      Begin
         FGIntadd(QFGInt, one, temp1);
         FGIntDestroy(QFGInt);
         FGIntCopy(temp1, QFGInt);
         FGIntDestroy(temp1);
         FGIntsub(FGInt2, MFGInt, temp1);
         FGIntDestroy(MFGInt);
         FGIntCopy(temp1, MFGInt);
         FGIntDestroy(temp1);
      End;
      If s2 = positive Then QFGInt.Sign := negative;
   End
   Else
      QFGInt.Sign := s2;
   FGIntDestroy(one);
   FGIntDestroy(zero);

   FGInt1.Sign := s1;
   FGInt2.Sign := s2;
End;

// Same as above but doesn 't compute MFGInt
Procedure FGIntDiv(Var FGInt1, FGInt2, QFGInt : TFGInt);
Var
   one, zero, temp1, temp2, MFGInt : TFGInt;
   s1, s2 : TSign;
   j, s, t : LongWord;
   i : int64;
Begin
   s1 := FGInt1.Sign;
   s2 := FGInt2.Sign;
   FGIntAbs(FGInt1);
   FGIntAbs(FGInt2);
   FGIntCopy(FGInt1, MFGInt);
   FGIntCopy(FGInt2, temp1);

   If FGIntCompareAbs(FGInt1, FGInt2) <> St Then
   Begin
      s := FGInt1.Number[0] - FGInt2.Number[0];
      SetLength(QFGInt.Number, (s + 2));
      QFGInt.Number[0] := s + 1;
      For t := 1 To s Do
      Begin
         FGIntShiftLeftBy31(temp1);
         QFGInt.Number[t] := 0;
      End;
      j := s + 1;
      QFGInt.Number[j] := 0;
      While FGIntCompareAbs(MFGInt, FGInt2) <> St Do
      Begin
         While FGIntCompareAbs(MFGInt, temp1) <> St Do
         Begin
            If MFGInt.Number[0] > temp1.Number[0] Then
            Begin
               i := MFGInt.Number[MFGInt.Number[0]];
               i := i Shl 31;
               i := i + MFGInt.Number[MFGInt.Number[0] - 1];
               i := i Div (temp1.Number[temp1.Number[0]] + 1);
            End
            Else
               i := MFGInt.Number[MFGInt.Number[0]] Div (temp1.Number[temp1.Number[0]] + 1);
            If (i <> 0) Then
            Begin
               FGIntCopy(temp1, temp2);
               FGIntMulByIntBis(temp2, i);
               FGIntSubBis(MFGInt, temp2);
               QFGInt.Number[j] := QFGInt.Number[j] + i;
               If FGIntCompareAbs(MFGInt, temp2) <> St Then
               Begin
                  QFGInt.Number[j] := QFGInt.Number[j] + i;
                  FGIntSubBis(MFGInt, temp2);
               End;
               FGIntDestroy(temp2);
            End
            Else
            Begin
               QFGInt.Number[j] := QFGInt.Number[j] + 1;
               FGIntSubBis(MFGInt, temp1);
            End;
         End;
         If MFGInt.Number[0] <= temp1.Number[0] Then
            If FGIntCompareAbs(temp1, FGInt2) <> Eq Then
            Begin
               FGIntShiftRightBy31(temp1);
               j := j - 1;
            End;
      End;
   End
   Else
      Base10StringToFGInt('0', QFGInt);
   s := QFGInt.Number[0];
   While (s > 1) And (QFGInt.Number[s] = 0) Do
      s := s - 1;
   If s < QFGInt.Number[0] Then
   Begin
      setlength(QFGInt.Number, s + 1);
      QFGInt.Number[0] := s;
   End;
   QFGInt.Sign := positive;

   FGIntDestroy(temp1);
   Base10StringToFGInt('0', zero);
   Base10StringToFGInt('1', one);
   If s1 = negative Then
   Begin
      If FGIntCompareAbs(MFGInt, zero) <> Eq Then
      Begin
         FGIntadd(QFGInt, one, temp1);
         FGIntDestroy(QFGInt);
         FGIntCopy(temp1, QFGInt);
         FGIntDestroy(temp1);
         FGIntsub(FGInt2, MFGInt, temp1);
         FGIntDestroy(MFGInt);
         FGIntCopy(temp1, MFGInt);
         FGIntDestroy(temp1);
      End;
      If s2 = positive Then QFGInt.Sign := negative;
   End
   Else
      QFGInt.Sign := s2;
   FGIntDestroy(one);
   FGIntDestroy(zero);
   FGIntDestroy(MFGInt);

   FGInt1.Sign := s1;
   FGInt2.Sign := s2;
End;


// Same as above but this computes MFGInt in stead of QFGInt
// MFGInt = FGInt1 mod FGInt2
Procedure FGIntMod(Var FGInt1, FGInt2, MFGInt : TFGInt);
Var
   one, zero, temp1, temp2 : TFGInt;
   s1, s2 : TSign;
   s, t : LongWord;
   i : int64;
Begin
   s1 := FGInt1.Sign;
   s2 := FGInt2.Sign;
   FGIntAbs(FGInt1);
   FGIntAbs(FGInt2);
   FGIntCopy(FGInt1, MFGInt);
   FGIntCopy(FGInt2, temp1);

   If FGIntCompareAbs(FGInt1, FGInt2) <> St Then
   Begin
      s := FGInt1.Number[0] - FGInt2.Number[0];
      For t := 1 To s Do
         FGIntShiftLeftBy31(temp1);
      While FGIntCompareAbs(MFGInt, FGInt2) <> St Do
      Begin
         While FGIntCompareAbs(MFGInt, temp1) <> St Do
         Begin
            If MFGInt.Number[0] > temp1.Number[0] Then
            Begin
               i := MFGInt.Number[MFGInt.Number[0]];
               i := i Shl 31;
               i := i + MFGInt.Number[MFGInt.Number[0] - 1];
               i := i Div (temp1.Number[temp1.Number[0]] + 1);
            End
            Else
               i := MFGInt.Number[MFGInt.Number[0]] Div (temp1.Number[temp1.Number[0]] + 1);
            If (i <> 0) Then
            Begin
               FGIntCopy(temp1, temp2);
               FGIntMulByIntBis(temp2, i);
               FGIntSubBis(MFGInt, temp2);
               If FGIntCompareAbs(MFGInt, temp2) <> St Then FGIntSubBis(MFGInt, temp2);
               FGIntDestroy(temp2);
            End
            Else
               FGIntSubBis(MFGInt, temp1);
//         If FGIntCompareAbs(MFGInt, temp1) <> St Then FGIntSubBis(MFGInt,temp1);
         End;
         If MFGInt.Number[0] <= temp1.Number[0] Then
            If FGIntCompareAbs(temp1, FGInt2) <> Eq Then FGIntShiftRightBy31(temp1);
      End;
   End;

   FGIntDestroy(temp1);
   Base10StringToFGInt('0', zero);
   Base10StringToFGInt('1', one);
   If s1 = negative Then
   Begin
      If FGIntCompareAbs(MFGInt, zero) <> Eq Then
      Begin
         FGIntSub(FGInt2, MFGInt, temp1);
         FGIntDestroy(MFGInt);
         FGIntCopy(temp1, MFGInt);
         FGIntDestroy(temp1);
      End;
   End;
   FGIntDestroy(one);
   FGIntDestroy(zero);

   FGInt1.Sign := s1;
   FGInt2.Sign := s2;
End;


// Square a FGInt modulo Modb, FGInt^2 mod Modb = FGIntSM
Procedure FGIntSquareMod(Var FGInt, Modb, FGIntSM : TFGInt);
Var
   temp : TFGInt;
Begin
   FGIntSquare(FGInt, temp);
   FGIntMod(temp, Modb, FGIntSM);
   FGIntDestroy(temp);
End;

// Multiply 2 FGInts modulo base, (FGInt1 * FGInt2) mod base = FGIntres
Procedure FGIntMulMod(Var FGInt1, FGInt2, base, FGIntres : TFGInt);
Var
   temp : TFGInt;
Begin
   FGIntMul(FGInt1, FGInt2, temp);
   FGIntMod(temp, base, FGIntres);
   FGIntDestroy(temp);
End;

// Procedures for Montgomery Exponentiation
Procedure FGIntModBis(Const FGInt : TFGInt; Var FGIntOut : TFGInt; b, head : LongWord);
Var
   i : LongWord;
Begin
   If b <= FGInt.Number[0] Then
   Begin
      SetLength(FGIntOut.Number, (b + 1));
      For i := 0 To b Do
         FGIntOut.Number[i] := FGInt.Number[i];
      FGIntOut.Number[b] := FGIntOut.Number[b] And head;
      i := b;
      While (FGIntOut.Number[i] = 0) And (i > 1) Do
         i := i - 1;
      If i < b Then SetLength(FGIntOut.Number, i + 1);
      FGIntOut.Number[0] := i;
      FGIntOut.Sign := positive;
   End
   Else
      FGIntCopy(FGInt, FGIntOut);
End;

Procedure FGIntMulModBis(Const FGInt1, FGInt2 : TFGInt; Var Prod : TFGInt; b, head : LongWord);
Var
   i, j, size, size1, size2, t, rest : LongWord;
   Trest : int64;
Begin
   size1 := FGInt1.Number[0];
   size2 := FGInt2.Number[0];
   size := min(b, size1 + size2);
   SetLength(Prod.Number, (size + 1));
   For i := 1 To size Do
      Prod.Number[i] := 0;

   For i := 1 To size2 Do
   Begin
      rest := 0;
      t := min(size1, b - i + 1);
      For j := 1 To t Do
      Begin
         Trest := FGInt1.Number[j];
	   Trest := Trest * FGInt2.Number[i];
	   Trest := Trest + Prod.Number[j + i - 1];
	   Trest := Trest + rest;
         Prod.Number[j + i - 1] := Trest And 2147483647;
         rest := Trest Shr 31;
      End;
      If (i + size1) <= b Then Prod.Number[i + size1] := rest;
   End;

   Prod.Number[0] := size;
   If size = b Then Prod.Number[b] := Prod.Number[b] And head;
   While (Prod.Number[size] = 0) And (size > 1) Do
      size := size - 1;
   If size < Prod.Number[0] Then
   Begin
      SetLength(Prod.Number, size + 1);
      Prod.Number[0] := size;
   End;
   If FGInt1.Sign = FGInt2.Sign Then
      Prod.Sign := Positive
   Else
      prod.Sign := negative;
End;

Procedure FGIntMontgomeryMod(Const GInt, base, baseInv : TFGInt; Var MGInt : TFGInt; b : Longword; head : LongWord);
Var
   m, temp, temp1 : TFGInt;
   r : LongWord;
Begin
   FGIntModBis(GInt, temp, b, head);
   FGIntMulModBis(temp, baseInv, m, b, head);
   FGIntMul(m, base, temp1);
   FGIntDestroy(temp);
   FGIntAdd(temp1, GInt, temp);
   FGIntDestroy(temp1);
   MGInt.Number := copy(temp.Number, b - 1, temp.Number[0] - b + 2);
   MGInt.Sign := positive;
   MGInt.Number[0] := temp.Number[0] - b + 1;
   FGIntDestroy(temp);
   If (head Shr 30) = 0 Then FGIntDivByIntBis(MGInt, head + 1, r)
   Else FGIntShiftRightBy31(MGInt);
   If FGIntCompareAbs(MGInt, base) <> St Then FGIntSubBis(MGInt, base);
   FGIntDestroy(temp);
   FGIntDestroy(m);
End;

Procedure FGIntMontgomeryModExp(Var FGInt, exp, modb, res : TFGInt);
Var
   temp2, temp3, baseInv, r, zero : TFGInt;
   i, j, t, b, head : LongWord;
   S : String;
Begin
   Base2StringToFGInt('0', zero);
   FGIntMod(FGInt, modb, res);
   If FGIntCompareAbs(res, zero)=Eq then
	Begin
	  FGIntDestroy(zero);
	  Exit;
	End else FGIntDestroy(res);
   FGIntDestroy(zero);

   FGIntToBase2String(exp, S);
   t := modb.Number[0];
   b := t;

   If (modb.Number[t] Shr 30) = 1 Then t := t + 1;
   SetLength(r.Number, (t + 1));
   r.Number[0] := t;
   r.Sign := positive;
   For i := 1 To t Do
      r.Number[i] := 0;
   If t = modb.Number[0] Then
   Begin
      head := 2147483647;
      For j := 29 Downto 0 Do
      Begin
         head := head Shr 1;
         If (modb.Number[t] Shr j) = 1 Then
         Begin
            r.Number[t] := 1 Shl (j + 1);
            break;
         End;
      End;
   End
   Else
   Begin
      r.Number[t] := 1;
      head := 2147483647;
   End;

   FGIntModInv(modb, r, temp2);
   If temp2.Sign = negative Then
      FGIntCopy(temp2, BaseInv)
   Else
   Begin
      FGIntCopy(r, BaseInv);
      FGIntSubBis(BaseInv, temp2);
   End;
//   FGIntBezoutBachet(r, modb, temp2, BaseInv);
   FGIntAbs(BaseInv);
   FGIntDestroy(temp2);
   FGIntMod(r, modb, res);
   FGIntMulMod(FGInt, res, modb, temp2);
   FGIntDestroy(r);

   For i := length(S) Downto 1 Do
   Begin
      If S[i] = '1' Then
      Begin
         FGIntmul(res, temp2, temp3);
         FGIntDestroy(res);
         FGIntMontgomeryMod(temp3, modb, baseinv, res, b, head);
         FGIntDestroy(temp3);
      End;
      FGIntSquare(temp2, temp3);
      FGIntDestroy(temp2);
      FGIntMontgomeryMod(temp3, modb, baseinv, temp2, b, head);
      FGIntDestroy(temp3);
   End;
   FGIntDestroy(temp2);
   FGIntMontgomeryMod(res, modb, baseinv, temp3, b, head);
   FGIntCopy(temp3, res);
   FGIntDestroy(temp3);
   FGIntDestroy(baseinv);
End;


// Compute the Greatest Common Divisor of 2 FGInts
Procedure FGIntGCD(Const FGInt1, FGInt2 : TFGInt; Var GCD : TFGInt);
Var
   k : TCompare;
   zero, temp1, temp2, temp3 : TFGInt;
Begin
   k := FGIntCompareAbs(FGInt1, FGInt2);
   If (k = Eq) Then FGIntCopy(FGInt1, GCD) Else
      If (k = St) Then FGIntGCD(FGInt2, FGInt1, GCD) Else
      Begin
         Base10StringToFGInt('0', zero);
         FGIntCopy(FGInt1, temp1);
         FGIntCopy(FGInt2, temp2);
         While (temp2.Number[0] <> 1) Or (temp2.Number[1] <> 0) Do
         Begin
            FGIntMod(temp1, temp2, temp3);
            FGIntCopy(temp2, temp1);
            FGIntCopy(temp3, temp2);
            FGIntDestroy(temp3);
         End;
         FGIntCopy(temp1, GCD);
         FGIntDestroy(temp2);
         FGIntDestroy(zero);
      End;
End;

// Perform a Rabin Miller Primality Test nrtest times on FGIntp, returns ok=true when FGIntp passes the test
Procedure FGIntRabinMiller(Var FGIntp : TFGInt; nrtest : Longword; Var ok : boolean);
Var
   j, b, i : LongWord;
   m, z, temp1, temp2, temp3, zero, one, two, pmin1 : TFGInt;
   ok1, ok2 : boolean;
Begin
   randomize;
   j := 0;
   Base10StringToFGInt('0', zero);
   Base10StringToFGInt('1', one);
   Base10StringToFGInt('2', two);
   FGIntsub(FGIntp, one, temp1);
   FGIntsub(FGIntp, one, pmin1);

   b := 0;
   While ((temp1.Number[1] Mod 2) = 0)and(FGIntCompareAbs(temp1, zero) <> eq) Do 
   Begin
      b := b + 1;
      FGIntShiftRight(temp1);
   End;
   m := temp1;

   i := 0;
   ok := true;
   Randomize;
   While (i < nrtest) And ok Do
   Begin
      i := i + 1;
      Base10StringToFGInt(inttostr(Primes[Random(1227) + 1]), temp2);
      FGIntMontGomeryModExp(temp2, m, FGIntp, z);
      FGIntDestroy(temp2);
      ok1 := (FGIntCompareAbs(z, one) = Eq);
      ok2 := (FGIntCompareAbs(z, pmin1) = Eq);
      If Not (ok1 Or ok2) Then
      Begin

         While (ok And (j < b)) Do
         Begin
            If (j > 0) And ok1 Then ok := false
            Else
            Begin
               j := j + 1;
               If (j < b) And (Not ok2) Then
               Begin
                  FGIntSquaremod(z, FGIntp, temp3);
                  FGIntCopy(temp3, z);
                  ok1 := (FGIntCompareAbs(z, one) = Eq);
                  ok2 := (FGIntCompareAbs(z, pmin1) = Eq);
                  If ok2 Then j := b;
               End
               Else If (Not ok2) And (j >= b) Then ok := false;
            End;
         End;

      End
   End;

   FGIntDestroy(zero);
   FGIntDestroy(one);
   FGIntDestroy(two);
   FGIntDestroy(m);
   FGIntDestroy(z);
   FGIntDestroy(pmin1);
End;

// Find the (multiplicative) Modular inverse of a FGInt in a finite ring
// of additive order base
Procedure FGIntModInv(Const FGInt1, base : TFGInt; Var Inverse : TFGInt);
Var
   zero, one, r1, r2, r3, tb, gcd, temp, temp1, temp2 : TFGInt;
Begin
   Base10StringToFGInt('1', one);
   FGIntGCD(FGInt1, base, gcd);
   If FGIntCompareAbs(one, gcd) = Eq Then
   Begin
      FGIntcopy(base, r1);
      FGIntcopy(FGInt1, r2);
      Base10StringToFGInt('0', zero);
      Base10StringToFGInt('0', inverse);
      Base10StringToFGInt('1', tb);

      Repeat
         FGIntDestroy(r3);
         FGIntdivmod(r1, r2, temp, r3);
         FGIntCopy(r2, r1);
         FGIntCopy(r3, r2);

         FGIntmul(tb, temp, temp1);
         FGIntsub(inverse, temp1, temp2);
         FGIntDestroy(inverse);
         FGIntDestroy(temp1);
         FGIntCopy(tb, inverse);
         FGIntCopy(temp2, tb);

         FGIntDestroy(temp);
      Until FGIntCompareAbs(r3, zero) = Eq;

      If inverse.Sign = negative Then
      Begin
         FGIntadd(base, inverse, temp);
         FGIntCopy(temp, inverse);
      End;

      FGIntDestroy(tb);
      FGIntDestroy(r1);
      FGIntDestroy(r2);
   End;
   FGIntDestroy(gcd);
   FGIntDestroy(one);
End;

End.

Unit UBrentRho.pas

unit UBrentRho;

interface

uses FGInt;

  function BrentRho(var A: TFGInt; var d: integer; var Z: TFGInt): string;

implementation

uses Forms, Math, SysUtils, Windows, Unit1;

  function BrentRho(var A: TFGInt; var d: integer; var Z: TFGInt): string;
  var
    okay: boolean;
    s: string;
    r, i, j, k, m: integer;
    xx, yy, tt, tmp, temp1, temp2, factorGGT, one: TFGInt; 
 
  label neu; 
 
  begin 
    if(Length(A.Number) = 0)then 
      exit;
    Base10StringToFGInt('0', xx);
    Base10StringToFGInt('0', yy);
    Base10StringToFGInt('0', temp1);
    Base10StringToFGInt('1', temp2);
    Base10StringToFGInt('1', factorGGT);
    Base10StringToFGInt('1', one);
    Base10StringToFGInt('0', tt);
    Base10StringToFGInt('0', tmp);
    result:='';

    TRY
      FGIntRabinMiller(A, 3, okay);
      if(okay)then begin
        FGIntToBase10String(A, s);
        result:= format('%s', [s]);
        d:=0;
      end;
      WHILE not(okay)DO BEGIN
neu:
        d:=0;
        r:=1;
        m:=1024;

        while(FGIntCompareAbs(factorGGT, one) = eq)do begin
          FGIntCopy(yy, xx);        //xx = yy
          for i:=1 to r do begin
             FGIntSquare(yy, tt);
             FGIntAdd(tt, Z, tmp);
             FGIntMod(tmp, A, yy);  //yy = (yy²+1) mod A
          end;
          k:=0;
          j:= min(m, r);
          repeat
            for i:=j downto 1 do begin
               FGIntSquare(yy, tt);
               FGIntAdd(tt, Z, tmp);
               FGIntMod(tmp, A, yy);
               FGIntSub(xx, yy, temp1);    //temp1 = xx - yy
               FGIntAbs(temp1);
               FGIntMul(temp1, temp2, tt);
               FGIntMod(tt, A, temp2);     //temp2 = (temp1*temp2) mod A
            end;

            FGIntGCD(temp2, A, factorGGT);
            d:=d + 1;
            k:=k + m;

            if(d and 7 = 0)then begin
              Application.ProcessMessages;
              QueryPerformanceCounter (t2x);       // Zählerstand 2
              form1.label1.caption:=FormatFloat('0.0 ms',1000*(t2x-t1x)/frequenz);
//Form1.lblZeit.Caption:= Format('D: %d  time: %.1f sec', [d, (t//- time)/1000]);
              if abbruch then begin
                Form1.Memo1.Lines.Add(Format('abgebrochen bei D: %d', [d]));
                exit;
              end;
            end;
          until((k >= r) or (FGIntCompareAbs(factorGGT, one) <> eq));
          r:= r shl 1;  //r:= r + r;
        end;

        if(FGIntCompareAbs(factorGGT, one) <> eq)then begin
          FGIntToBase10String(factorGGT, s);
          FGIntDiv(A, factorGGT, A);

          if(FGIntCompareAbs(A, one) = eq)then begin
            FGIntCopy(factorGGT, A);
            FGIntAdd(Z, one, Z);      //anderes Z versuchen !!
            Base10StringToFGInt('1', temp2);
            Base10StringToFGInt('1', factorGGT);
            goto neu;
          end else
            form1.memo1.lines.add('Teiler = '+s);
//          result:= format('%s%s * ', [result, s]);

          FGIntRabinMiller(A, 3, okay);
          if(okay)then begin
            FGIntToBase10String(A, s);
            form1.memo1.lines.add('Teiler = '+s);
//          result:= format('%s%s', [result, s]);
          end else
            Base10StringToFGInt('1', factorGGT);

        end;
      END;
      FGIntToBase10String(Z, s);
//      Form1.lblZeit.Caption:= Format('D: %d    Z: %s    time: %.2f sec   fertig!',
//              [d, s, (GetTickCount - time)/1000]);

    FINALLY
      FGIntDestroy(xx);
      FGIntDestroy(yy);
      FGIntDestroy(temp1);
      FGIntDestroy(temp2);
      FGIntDestroy(factorGGT);
      FGIntDestroy(one);
      FGIntDestroy(tt); 
      FGIntDestroy(tmp); 
    END; 
  end;    //BrentRhoTFGInt 

end.

Unit1 :

uses fgint, ubrentrho

var
  Form1: TForm1;
  abbruch:boolean;
  t1x,t2x,frequenz : TLargeInteger;
  
//

function IstPrime(N: Cardinal): Boolean;
//Primzahltest nach Hagen Reddmann
const
  MaxPrime = 137;
  MaxTrial = MaxPrime * MaxPrime; // maximal trial div bounds
  MinPrime = 137;                 // trial div bounds if N >= MaxTrial, must be <= MaxPrime 
  Primes: array[0..31] of Byte =
   (   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,
      61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137);

  InvPrimes: array[0..31] of Cardinal =  // InvPrimes[i] = Primes[i]^-1 mod 2^32 
   ($AAAAAAAB,$CCCCCCCD,$B6DB6DB7,$BA2E8BA3,$C4EC4EC5,$F0F0F0F1,$286BCA1B,$E9BD37A7,
    $4F72C235,$BDEF7BDF,$914C1BAD,$C18F9C19,$2FA0BE83,$677D46CF,$8C13521D,$A08AD8F3,
    $C10C9715,$07A44C6B,$E327A977,$C7E3F1F9,$613716AF,$2B2E43DB,$FA3F47E9,$5F02A3A1, 
    $7C32B16D,$D3431B57,$8D28AC43,$DA6C0965,$0FDBC091,$EFDFBF7F,$C9484E2B,$077975B9);
// Montgomery Version of IsPrimeHR_IA32_1()
asm
       TEST   AL,1
       JZ     @@0                          // even ??
       CMP    EAX,7
       JA     @@1                          // > 7 ?? 
       DEC    EAX 
       SETNZ  AL 
       RET

@@0:   CMP    EAX,2                        // even numbers 
       SETZ   AL 
       RET 
 
@@1:   PUSH   EBP                          // do trial divsion to small primes 
       PUSH   EBX
       PUSH   ESI 
       PUSH   EDI 

       MOV    EBX,EAX 
       CMP    EAX,MaxTrial 
       MOV    EBP,MinPrime 
       JAE    @@2 
 
       PUSH   EAX
       FILD   DWord Ptr [ESP] 
       FSQRT 
       FISTP  DWord Ptr [ESP]
       POP    EBP                          // EBP = Sqrt(N) 
 
@@2:   MOV    EDI,Offset Primes
       MOV    ESI,Offset InvPrimes
       XOR    ECX,ECX
@@3:   MOVZX  EDX,Byte Ptr [EDI + ECX]     // take care, InvPrimes[] MUST 
                                           // be after Primes[] declared
       MOV    EAX,EBX 
       CMP    EDX,EBP 
       JA     @@5 
       IMUL   EAX,[ESI + ECX * 4] 
       INC    ECX 
       MUL    EDX 
       JC     @@3 
       TEST   EDX,EDX 
@@4:   POP    EDI 
       POP    ESI 
       POP    EBX 
       POP    EBP
       SETNZ  AL 
       RET 
 
@@5:   CMP    EBX,MaxTrial                 // N <= MaxPrime^2 ?? 
       MOV    EAX,EBX 
       JBE    @@4

       IMUL   EAX,EBX                      // compute domain param U = -N^-1 mod 2^32 
       SUB    EAX,2                        // Lookuptable can reduce from 72 donwto 32 cycles 
       IMUL   EAX,EBX 
       MOV    EDX,EAX 
       IMUL   EAX,EBX 
       ADD    EAX,2 
       IMUL   EAX,EDX
       MOV    EDX,EAX 
       IMUL   EAX,EBX 
       ADD    EAX,2 
       IMUL   EAX,EDX 
       MOV    EBP,EAX
       IMUL   EBP,EBX 
       ADD    EBP,2 
       IMUL   EBP,EAX                      // U = -N^-1 mod^2^32 in EBP
 
       MOV    EDI,EBX 
       MOV    EAX,EBX 
       DEC    EDI                          // N -1
       NEG    EAX 
       BSF    ECX,EDI
       MUL    EAX
       PUSH   ECX                          // bits remain         [ESP + 20] 
       MOV    ESI,EAX 
       BSR    ECX,EDI 
       MOV    EAX,EDX 
       XOR    EDX,EDX 
       NEG    ECX 
       DIV    EBX                          // div, can't be removed 
       MOV    EAX,ESI 
       ADD    ECX,32
       DIV    EBX                          // div 
       SHL    EDI,CL 
       MOV    EAX,EDX 
       IMUL   EAX,EBP
       MOV    ESI,EDX                      // C = -N^2 mod N, to fast convert into 
       MUL    EBX                          // montgomery domain
       PUSH   ESI                          // C                    [ESP + 16]
       ADD    EAX,ESI 
       ADC    EDX,0 
 
       PUSH   EDX                          // +1 in montgomery     [ESP + 12] 
       PUSH   EDI                          // bit mask exponent    [ESP +  8]
       NEG    EDX
       ADD    EDX,EBX 
 
       CMP    EBX,$08A8D7F                 // N < $08A8D7F, do SPP to bases 31,73 
       PUSH   EDX                          // -1 in montgomery     [ESP +  4] 
       JAE    @@6
       MOV    EAX,31
       CALL   @@9 
       MOV    EAX,73 
       PUSH   Offset @@7 
       JMP    @@9 
 
@@6:   MOV    EAX,2                        // do SPP to bases 2,7,61 
       CALL   @@9
       MOV    EAX,7 
       CALL   @@9
       MOV    EAX,61 
       CALL   @@9 
 
@@7:   INC    EAX
@@8:   LEA    ESP,[ESP + 4 * 5]            // don't change flags !! 
       JMP    @@4

@@9:   MUL    DWord Ptr [ESP + 16]         // convert base in montgomery 
       MOV    EDI,EAX                      // Base' = Base * C mod N 
       IMUL   EAX,EBP                      // montgomery REDC
       MOV    ESI,EDX 
       MUL    EBX 
       ADD    EAX,EDI 
       ADC    EDX,ESI 
       MOV    ECX,[ESP + 8]                // bit mask of exponent N -1 
       MOV    EAX,EDX 
       PUSH   EDX

@@A:   MUL    EAX                          // X = X^2 mod N 
       MOV    EDI,EAX 
       IMUL   EAX,EBP 
       MOV    ESI,EDX
       MUL    EBX 
       ADD    EAX,EDI 
       ADC    EDX,ESI 
       JNC    @@B 
       SUB    EDX,EBX 
@@B:   ADD    ECX,ECX
       MOV    EAX,EDX
       JNC    @@D 
       MUL    DWord Ptr [ESP]              // X = X * Base mod N 
       MOV    EDI,EAX
       IMUL   EAX,EBP 
       MOV    ESI,EDX 
       MUL    EBX 
       ADD    EAX,EDI
       ADC    EDX,ESI 
       JNC    @@C 
       SUB    EDX,EBX 
@@C:   TEST   ECX,ECX 
       MOV    EAX,EDX 
@@D:   JNZ    @@A 
       CMP    EAX,EBX 
       JB     @@E
       SUB    EAX,EBX 
@@E:   CMP    EAX,[ESP + 16]               // == +1 ??
       MOV    ECX,[ESP + 24]               // bits remain 
       JE     @@J
@@F:   CMP    EAX,[ESP + 8]                // == -1 ?? 
       JE     @@J
       DEC    ECX
       JNG    @@I 
       MUL    EAX 
       MOV    EDI,EAX 
       IMUL   EAX,EBP
       MOV    ESI,EDX 
       MUL    EBX
       ADD    EAX,EDI
       ADC    EDX,ESI
       JC     @@G
       CMP    EDX,EBX
       JB     @@H
@@G:   SUB    EDX,EBX
@@H:   CMP    EDX,[ESP + 16]               // <> +1 ??
       MOV    EAX,EDX
       JNE    @@F
@@I:   ADD    ESP,8
       XOR    EAX,EAX
       JMP    @@8
@@J:   POP    EDX
end;

procedure TForm1.Button1Click(Sender: TObject);
label 1;
var k:string;
    i,probeprim,d:integer;
    proberest:longword;
    pa,pb,pc,x,y,konstante,sum,a,zz:tfgint;
    prim,kontrolle:boolean;
begin
    if abbruch=false then begin abbruch:=true; exit end;
    button1.caption:='Abbruch';
    label3.caption:='';
    abbruch:=false;
    memo1.clear;
    application.processmessages;
    k:=edit1.text;
    i:=1;
    repeat
      if (length(k)>=i) and (not (k[i] in ['0'..'9'])) then delete(k,i,1)
                                                       else inc(i);
    until i>length(k);
    if (k='0') then k:='1';
    memo1.lines.add('Zahl = '+k);

    QueryPerformanceFrequency (frequenz);      // Frequenz des Zählers
    QueryPerformanceCounter (t1x);       // Zählerstand 1

    //Probedivisionen
    probeprim:=2;
    Base10StringToFGInt(k,pa);
    FGIntAbs(pa);
    FGIntModByInt(pa,probeprim,proberest);
    while proberest=0 do
    begin
      FGIntDivByIntBis(pa,probeprim,proberest);
      memo1.lines.add('Teiler = 2');
      FGIntToBase10String(pa,k);
      FGIntModByInt(pa,probeprim,proberest);
    end;

    probeprim:=3;
    Base10StringToFGInt(k,pa);
    repeat
      FGIntModByInt(pa,probeprim,proberest);
      if proberest=0 then
      begin
        while proberest=0 do
        begin
          FGIntDivByIntBis(pa,probeprim,proberest);
          memo1.lines.add('Teiler = '+inttostr(probeprim));
          FGIntToBase10String(pa,k);
          FGIntModByInt(pa,probeprim,proberest);
        end
      end
      else
      begin
        inc(probeprim,2);
        while not istprime(probeprim) do inc(probeprim)
      end;
    until probeprim>10000;
    QueryPerformanceCounter (t2x);       // Zählerstand 2
    label1.caption:=FormatFloat('0.0 ms',1000*(t2x-t1x)/frequenz);

//    memo1.lines.add('');
//    memo1.lines.add('Restzahl = '+k);

    Base10StringToFGInt(k,pa);
    Base10StringToFGInt('2',pb);
    Base10StringToFGInt('2',pc);
    Base10StringToFGInt('1',konstante);
    Base10StringToFGInt('1',sum);

  repeat
    //Test auf = 1
    kontrolle:=(FGIntCompareAbs(pa,konstante)=eq);
    if kontrolle=false then begin

    //Primzahltest
    FGIntRabinMiller(pa,1,prim);

    FGIntToBase10String(pa,k);
    if prim then memo1.lines.add('Teiler = '+k)
    else
    begin

      Base10StringToFGInt(k, A);
      Base10StringToFGInt('1', Zz);
      //Brent-Rho-Verfahren
      k:= BrentRho(A, d, Zz);

      if not abbruch then begin
        Base10StringToFGInt('1', pa);
      end else begin
        memo1.lines.add('Abbruch');
        FGIntToBase10String(pa,k);
        memo1.lines.add('Restzahl = '+k);
      end;
    end;
   end;
   until kontrolle or prim or abbruch;

1:  if not abbruch then
      memo1.lines.add('Zahl faktorisiert');
    QueryPerformanceCounter (t2x);       // Zählerstand 2
    form1.label1.caption:=FormatFloat('0.0 ms',1000*(t2x-t1x)/frequenz);
    FGIntDestroy(pa);
    FGIntDestroy(pb);
    FGIntDestroy(pc);
    FGIntDestroy(x);
    FGIntDestroy(y);
    FGIntDestroy(konstante);
    FGIntDestroy(sum);
    FGIntDestroy(a);
    FGIntDestroy(zz);
    button1.caption:='Berechnung';
    abbruch:=true;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
    abbruch:=true;
end;