// Catalan(1.000.000) in 810 ms (Athlon X4 640, XP, 2GB). // Amsterdam University Library thanks, there I read // "Computing Binomial Coefficients" (April 1987!) by P. Goetgheluck. // The blog about binomials (december 2010) has to be updated. // "BNM(uint n, uint k)" chooses the fastest algorithm, BNM1 or BNM3. // If k is small compared to n, BNM1 will be used. // Example: BNM1(12,4)=12!/8!/4!=9*10*11*12/(2*3*4)=9*5*11*6/(3*2) // Otherwise BNM3 will be used, it uses prime power factorization. // A sieve of Eratosthenes "getComposites" returns a BitArray, // and an approximation for the square root of n. // The exact square root is found by "root". // Goetgheluck's algorithm is used in "getBnmPrimes", // to get the power "exp(n,k,p)" of a prime. // A signed integer version "exp(int n, int k, int p)" isn't used, // but ..... is there for it's own sake. // Some care has been taken to work with arrays of uints, ulongs etc as long as possible. // Result: BNM(6.400.000 , 2.133.333) in 3510 ms, that's nearly 10 % faster. using Xint = System.Numerics.BigInteger; using System.Collections.Generic; using System.Threading.Tasks; using System.Collections; using System.Diagnostics; using System; class Binomial { public static Xint CAT(uint n) { return n < 3 ? n / 2 + 1 : n < 11 ? BNM1(n * 2, n++) / n : BNM3(n * 2, n++) / n; } public static Xint BNM(uint n, uint k) { if (k > n) return 0; if (k > n / 2) k = n - k; if (k == 0) return 1; if (k == 1) return n; if (k == 2) return n <= ushort.MaxValue ? n-- * n / 2 : (ulong)(n--) * n / 2; if (k < 11) return BNM1(n, k); if (n <= 64) return BNM3(n, k); if (n <= 128) return k < 12 ? BNM1(n, k) : BNM3(n, k); if (n <= 256) return k < 12 + 5 * (n - 128) / 128 ? BNM1(n, k) : BNM3(n, k); if (n <= 512) return k < 17 + 2 * (n - 256) / 256 ? BNM1(n, k) : BNM3(n, k); if (n <= 1024) return k < 19 + 9 * (n - 512) / 512 ? BNM1(n, k) : BNM3(n, k); if (n <= 2048) return k < 28 + 13 * (n - 1024) / 1024 ? BNM1(n, k) : BNM3(n, k); if (n <= 4096) return k < 41 + 17 * (n - 2048) / 2048 ? BNM1(n, k) : BNM3(n, k); if (n <= 8192) return k < 58 + 31 * (n - 4096) / 4096 ? BNM1(n, k) : BNM3(n, k); if (n <= 16384) return k < 89 + 45 * (n - 8192) / 8192 ? BNM1(n, k) : BNM3(n, k); if (n <= 32768) return k < 134 + 60 * (n - 16384) / 16384 ? BNM1(n, k) : BNM3(n, k); if (n <= 65536) return k < 194 + 88 * (n - 32768) / 32768 ? BNM1(n, k) : BNM3(n, k); if (n <= 131072) return k < 282 + 120 * (n - 65536) / 65536 ? BNM1(n, k) : BNM3(n, k); if (n <= 262144) return k < 402 + 162 * (n - 131072) / 131072 ? BNM1(n, k) : BNM3(n, k); if (n <= 524288) return k < 564 + 229 * (n - 262144) / 262144 ? BNM1(n, k) : BNM3(n, k); if (n <= 1048576) return k < 793 + 312 * (n - 524288) / 524288 ? BNM1(n, k) : BNM3(n, k); if (n <= 2097152) return k < 1105 + 430 * (n - 1048576) / 1048576 ? BNM1(n, k) : BNM3(n, k); if (n <= 4194304) return k < 1535 + 583 * (n - 2097152) / 2097152 ? BNM1(n, k) : BNM3(n, k); if (n <= 8388608) return k < 2118 + 806 * (n - 4194304) / 4194304 ? BNM1(n, k) : BNM3(n, k); if (n <= 16777216) return k < 2924 + (ulong)(1700) * (n - 8388608) / 8388608 ? BNM1(n, k) : BNM3(n, k); if (n <= 33554432) return k < 4624 + (ulong)(2328) * (n - 16777216) / 16777216 ? BNM1(n, k) : BNM3(n, k); if (n <= 67108864) return k < 6952 + (ulong)(3135) * (n - 33554432) / 33554432 ? BNM1(n, k) : BNM3(n, k); if (n <= 134217728) return k < 10087 + (ulong)(4282) * (n - 67108864) / 67108864 ? BNM1(n, k) : BNM3(n, k); return k < 10087 + 4282 * (n >> 13) / 16384 ? BNM1(n, k) : BNM3(n, k); } private static Xint BNM1(uint n, uint k) { ulong uu = 0; Xint U = 0; uint i = n - k + 1; uint u = i; if ((i++ & 1) == 0) { u /= 2; goto L1; } L0: if (i <= n && u <= ushort.MaxValue && i <= ushort.MaxValue) { u *= (i++ / 2); } else { if (i > n) goto M0; else { uu = (ulong)(u) * (i++ / 2); goto L2; } } L1: if (i <= n && u <= ushort.MaxValue && i <= ushort.MaxValue) { u *= i++; goto L0; } else { if (i > n) goto M0; else { uu = (ulong)(u) * i++; goto L3; } } L2: if (i <= n && uu <= uint.MaxValue) { uu *= i++; } else { if (i > n) goto M0; else { U = (Xint)(uu) * i++; goto L4; } } L3: if (i <= n && uu <= uint.MaxValue) { uu *= (i++ / 2); goto L2; } else { if (i > n) goto M0; else { U = (Xint)(uu) * (i++ / 2); goto L5; } } L4: if (i <= n) { U *= (i++ / 2); } else goto M0; L5: if (i <= n) { U *= i++; goto L4; } M0: switch (k) { case 3: if (uu == 0) return u / 3 << 1 - (int)(n & 1); if (U.IsZero) return uu / 3 << 1 - (int)(n & 1); return U / 3 << 1 - (int)(n & 1); case 4: if (uu == 0) return u / 6; if (U.IsZero) return uu / 6; return U / 6; case 5: if (uu == 0) return (u >> (int)(n & 1)) / 15; if (U.IsZero) return (uu >> (int)(n & 1)) / 15; return (U >> (int)(n & 1)) / 15; case 6: if (uu == 0) return u / 90; if (U.IsZero) return uu / 90; return U / 90; case 7: if (U.IsZero) return (uu >> (int)(n & 1)) / 315; return (U >> (int)(n & 1)) / 315; case 8: if (U.IsZero) return uu / 2520; return U / 2520; case 9: if (U.IsZero) return (uu >> (int)(n & 1)) / 11340; return (U >> (int)(n & 1)) / 11340; case 10: if (U.IsZero) return uu / 113400; return U / 113400; case 11: return (U >> (int)(n & 1)) / 623700; case 12: return U / 7484400; case 13: return (U >> (int)(n & 1)) / 48648600; case 14: return U / 681080400; case 15: return (U >> 3 + (int)(n & 1)) / 638512875; case 16: return (U >> 7) / 638512875; case 17: return (U >> 6 + (int)(n & 1)) / 10854718875; case 18: return (U >> 7) / 97692469875; case 19: return (U >> 6 + (int)(n & 1)) / 1856156927625; case 20: return (U >> 8) / 9280784638125; case 21: return (U >> 7 + (int)(n & 1)) / 194896477400625; case 22: return (U >> 8) / 2143861251406875; case 23: return (U >> 7 + (int)(n & 1)) / 49308808782358125; case 24: return (U >> 10) / 147926426347074375; case 25: return (U >> 9 + (int)(n & 1)) / 3698160658676859375; default: Xint V = 3698160658676859375; i = 26; M1: if (i <= k) V *= (i++ / 2); else return (U >> 9 + (int)(n & 1)) / V; if (i <= k) { V *= i++; goto M1; } return (U >> 10) / V; } } private static Xint BNM3(uint n, uint k) { if (n <= 65536) { uint[] f = getBnmPrimes(n, k); int i = f.Length; if (i == 1) return (Xint)(f[0]) << exp(n, k, 2); int j = 1 << fL2(i); uint max = 0; if (i != j) { f = uintSpecialProduct(f, i, j, out max); } for (; j > 1 && max <= ushort.MaxValue; j /= 2) f = uintProduct(f, j, out max); if (j == 1) return (Xint)(f[0]) << exp(n, k, 2); ulong mmax = 0; ulong[] ff = ulongProduct(f, j, out mmax); j /= 2; if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2); Array.Resize(ref f, 0); for (; j > 1 && mmax <= uint.MaxValue; j /= 2) ff = ulongProduct(ff, j, out mmax); if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2); Xint[] F = XintProduct(ff, j); j /= 2; if (j == 1) return F[0] << exp(n, k, 2); Array.Resize(ref ff, 0); for (; j > 1 && F[0].ToByteArray().Length < 376; j /= 2) F = SmallProduct(F, j); for (; j > 1; j /= 2) F = LargeProduct(F, j); return F[0] << exp(n, k, 2); } else { uint[] f = getBnmPrimes(n, k); int i = f.Length; if (i == 1) return (Xint)(f[0]) << exp(n, k, 2); int j = 1 << fL2(i); ulong mmax = 0; ulong[] ff; if (i != j) { ff = ulongSpecialProduct(f, i, j, out mmax); } else { ff = ulongProduct(f, j, out mmax); j /= 2; } if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2); Array.Resize(ref f, 0); for (; j > 1 && mmax <= uint.MaxValue; j /= 2) ff = ulongProduct(ff, j, out mmax); if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2); Xint[] F = XintProduct(ff, j); j /= 2; if (j == 1) return F[0] << exp(n, k, 2); Array.Resize(ref ff, 0); for (; j > 1 && F[0].ToByteArray().Length < 376; j /= 2) F = SmallProduct(F, j); for (; j > 1; j /= 2) F = LargeProduct(F, j); return F[0] << exp(n, k, 2); } } private static uint[] getBnmPrimes(uint n, uint k) { int i; BitArray composite = getComposites(n, out i); uint rt = root(n, i); List<uint> primes = new List<uint>(); int e = exp(n, k, 3); for (i = 0; i < e; i++) primes.Add(3); uint p = 5; i = 0; int j = 0; L0: if (p > rt) { e = (int)(n / 2); goto L1; } if (!composite[i++]) { e = exp(n, k, p); for (j = 0; j < e; j++) primes.Add(p); } p += 2; if (p > rt) { e = (int)(n / 2); goto L2; } if (!composite[i++]) { e = exp(n, k, p); for (j = 0; j < e; j++) primes.Add(p); } p += 4; goto L0; L1: if (p > e) goto L3; if (!composite[i++] && n % p < k % p) primes.Add(p); p += 2; L2: if (p > e) goto L3; if (!composite[i++] && n % p < k % p) primes.Add(p); p += 4; goto L1; L3: p = n - k + 1; p += 1 - (p & 1); p += p % 6 & 2; // if (p % 6 == 3) p += 2; i = (int)(p / 3 - 1); if ((i & 1) != 0) goto L5; L4: if (p > n) goto L6; if (!composite[i++]) primes.Add(p); p += 2; L5: if (p > n) goto L6; if (!composite[i++]) primes.Add(p); p += 4; goto L4; L6: return primes.ToArray(); } private static BitArray getComposites(uint n, out int i) { int len = (int)(n / 3); BitArray composite = new BitArray(len); i = 0; int d1 = 8, d2 = 0, p1 = 3, p2 = 1, s = 7, inc = 10; while (s < len) { if (!composite[i++]) { int j = s, m = s + p1; for (; m < len; j += inc, m += inc) { composite[j] = true; composite[m] = true; } if (j < len) composite[j] = true; } d2 += 8; s += d2; p2 += 8; inc += 4; if (!composite[i++]) { int j = s, m = s + p2; for (; m < len; j += inc, m += inc) { composite[j] = true; composite[m] = true; } if (j < len) composite[j] = true; } d1 += 16; s += d1; p1 += 4; inc += 8; } return composite; } private static uint root(uint n, int i) { uint rt = (uint)(i * 3 + 3); uint sq = rt * rt; for (; sq < n; sq += 2 * rt | 1, rt++) ; for (; sq > n; rt--, sq -= 2 * rt | 1) ; return rt; } private static int exp(uint n, uint k, uint p) { int e = 0; L0: uint qn = n / p; uint qk = k / p; if ((qn - qk) * p > n - k) { e++; n = qn; if (qk == 0) goto L2; k = qk; goto L1; } if (qk == 0) return e; n = qn; k = qk; goto L0; L1: qn = n / p; qk = k / p; if ((qn - qk) * p >= n - k) { e++; n = qn; if (qk == 0) goto L2; k = qk; goto L1; } if (qk == 0) return e; n = qn; k = qk; goto L0; L2: qn = n / p; if (n == qn * p) e++; else return e; n = qn; goto L2; } private static int exp(int n, int k, int p) { int e = 0, rn = 0, rk = 0; L0: n = Math.DivRem(n, p, out rn); k = Math.DivRem(k, p, out rk); if (rk > rn) { e++; if (k == 0) goto L2; goto L1; } if (k == 0) return e; goto L0; L1: n = Math.DivRem(n, p, out rn); k = Math.DivRem(k, p, out rk); if (rk >= rn) { e++; if (k == 0) goto L2; goto L1; } if (k == 0) return e; goto L0; L2: n = Math.DivRem(n, p, out rn); if (rn == 0) e++; else return e; goto L2; } private static uint[] uintSpecialProduct(uint[] f, int i, int j, out uint max) { max = 0; uint p = 0; uint[] newF = new uint[j]; int m = i - j; // ?!?! i = j - m; // ?!?! int n = 0; while (n < i) newF[n++] = f[m++]; i = 0; while (n < j) { p = f[m++] * f[i++]; if (max < p) max = p; newF[n++] = p; } return newF; } private static ulong[] ulongSpecialProduct(uint[] f, int i, int j, out ulong mmax) { mmax = 0; ulong pp = 0; ulong[] newF = new ulong[j]; int m = i - j; i = j - m; int n = 0; while (n < i) newF[n++] = f[m++]; i = 0; while (n < j) { pp = (ulong)(f[m++]) * f[i++]; if (mmax < pp) mmax = pp; newF[n++] = pp; } return newF; } private static uint[] uintProduct(uint[] f, int j, out uint max) { max = 0; uint p = 0; int k = j-- / 2; uint[] newF = new uint[k]; for (int i = 0; i < k; i++, j--) { p = f[i] * f[j]; if (max < p) max = p; newF[i] = p; } return newF; } private static ulong[] ulongProduct(uint[] f, int j, out ulong mmax) { mmax = 0; ulong pp = 0; int k = j-- / 2; ulong[] newF = new ulong[k]; for (int i = 0; i < k; i++, j--) { pp = (ulong)(f[i]) * f[j]; if (mmax < pp) mmax = pp; newF[i] = pp; } return newF; } private static ulong[] ulongProduct(ulong[] ff, int j, out ulong mmax) { mmax = 0; ulong pp = 0; int k = j-- / 2; ulong[] newF = new ulong[k]; for (int i = 0; i < k; i++, j--) { pp = ff[i] * ff[j]; if (mmax < pp) mmax = pp; newF[i] = pp; } return newF; } private static Xint[] XintProduct(ulong[] ff, int j) { int k = j-- / 2; Xint[] NewF = new Xint[k]; for (int i = 0; i < k; i++, j--) NewF[i] = (Xint)(ff[i]) * ff[j]; return NewF; } private static Xint[] SmallProduct(Xint[] F, int j) { int k = j-- / 2; Xint[] NewF = new Xint[k]; for (int i = 0; i < k; i++, j--) NewF[i] = F[i] * F[j]; return NewF; } private static Xint[] LargeProduct(Xint[] F, int j) { int k = j-- / 2; Xint[] NewF = new Xint[k]; for (int i = 0; i < k; i++, j--) NewF[i] = MTP(F[i], F[j]); return NewF; } public static int bL(Xint U) { byte[] bytes = (U.Sign * U).ToByteArray(); int i = bytes.Length - 1; return i * 8 | bitLengthMostSignificantByte(bytes[i]); } private static int bitLengthMostSignificantByte(byte b) { return b < 08 ? b < 02 ? b < 01 ? 0 : 1 : b < 04 ? 2 : 3 : b < 32 ? b < 16 ? 4 : 5 : b < 64 ? 6 : 7; } private static Xint MTP(Xint U, Xint V) { return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3); } private static Xint MTP(Xint U, Xint V, int n) { if (n <= 3000) return U * V; if (n <= 6000) return TC2(U, V, n); if (n <= 10000) return TC3(U, V, n); if (n <= 40000) return TC4(U, V, n); return TC2P(U, V, n); } private static Xint MTPr(Xint U, Xint V, int n) { if (n <= 3000) return U * V; if (n <= 6000) return TC2(U, V, n); if (n <= 10000) return TC3(U, V, n); return TC4(U, V, n); } private static Xint TC2(Xint U1, Xint V1, int n) { n >>= 1; Xint Mask = (Xint.One << n) - 1; Xint U0 = U1 & Mask; U1 >>= n; Xint V0 = V1 & Mask; V1 >>= n; Xint P0 = MTPr(U0, V0, n); Xint P2 = MTPr(U1, V1, n); return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0; } private static Xint TC3(Xint U2, Xint V2, int n) { n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3; Xint Mask = (Xint.One << n) - 1; Xint U0 = U2 & Mask; U2 >>= n; Xint U1 = U2 & Mask; U2 >>= n; Xint V0 = V2 & Mask; V2 >>= n; Xint V1 = V2 & Mask; V2 >>= n; Xint W0 = MTPr(U0, V0, n); Xint W4 = MTPr(U2, V2, n); Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n); U2 += U0; V2 += V0; Xint P2 = MTPr(U2 - U1, V2 - V1, n); Xint P1 = MTPr(U2 + U1, V2 + V1, n); Xint W2 = (P1 + P2 >> 1) - (W0 + W4); Xint W3 = W0 - P1; W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1); Xint W1 = P1 - (W4 + W3 + W2 + W0); return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint TC4(Xint U3, Xint V3, int n) { n >>= 2; Xint Mask = (Xint.One << n) - 1; Xint U0 = U3 & Mask; U3 >>= n; Xint U1 = U3 & Mask; U3 >>= n; Xint U2 = U3 & Mask; U3 >>= n; Xint V0 = V3 & Mask; V3 >>= n; Xint V1 = V3 & Mask; V3 >>= n; Xint V2 = V3 & Mask; V3 >>= n; Xint W0 = MTPr(U0, V0, n); // 0 U0 += U2; U1 += U3; V0 += V2; V1 += V3; Xint P1 = MTPr(U0 + U1, V0 + V1, n); // 1 Xint P2 = MTPr(U0 - U1, V0 - V1, n); // -1 U0 += 3 * U2; U1 += 3 * U3; V0 += 3 * V2; V1 += 3 * V3; Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n); // 2 Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n); // -2 Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); // 4 Xint W6 = MTPr(U3, V3, n); // inf Xint W1 = P1 + P2; Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6; Xint W2 = (W1 >> 1) - (W6 + W4 + W0); P1 = P1 - P2; P4 = P4 - P3; Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45; W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2); Xint W3 = (P1 >> 1) - (W1 + W5); return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint TC2P(Xint A, Xint B, int n) { n >>= 1; Xint Mask = (Xint.One << n) - 1; Xint[] U = new Xint[3]; U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A; Xint[] V = new Xint[3]; V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B; Xint[] P = new Xint[3]; Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n)); return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0]; } private static int fL2(int i) { return i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 : i < 1 << 02 ? 01 : 02 : i < 1 << 05 ? i < 1 << 04 ? 03 : 04 : i < 1 << 06 ? 05 : 06 : i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 : i < 1 << 10 ? 09 : 10 : i < 1 << 13 ? i < 1 << 12 ? 11 : 12 : i < 1 << 14 ? 13 : 14 : i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 : i < 1 << 18 ? 17 : 18 : i < 1 << 21 ? i < 1 << 20 ? 19 : 20 : i < 1 << 22 ? 21 : 22 : i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 : i < 1 << 26 ? 25 : 26 : i < 1 << 29 ? i < 1 << 28 ? 27 : 28 : i < 1 << 30 ? 29 : 30; } private static Stopwatch sw = new Stopwatch(); static void Main() { BNM(1000000, 500000); sw.Restart(); BNM(6400000, 2133333); sw.Stop(); Console.WriteLine("BNM(6.400.000 , 2.133.333)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(); sw.Restart(); Xint C = CAT(100); sw.Stop(); Console.WriteLine("Catalan(100)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(bL(C) + " bits"); Console.WriteLine(C); Console.WriteLine(); sw.Restart(); C = CAT(100); sw.Stop(); Console.WriteLine("Catalan(100)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(bL(C) + " bits"); Console.WriteLine(C); Console.WriteLine(); sw.Restart(); C = CAT(200); sw.Stop(); Console.WriteLine("Catalan(200)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(bL(C) + " bits"); Console.WriteLine(C); Console.WriteLine(); sw.Restart(); for (int i = 0; i < 1000; i++) { C = CAT(1000); } sw.Stop(); Console.WriteLine("Catalan(1000)"); Console.WriteLine(sw.ElapsedMilliseconds + " us"); Console.WriteLine(bL(C) + " bits"); Console.WriteLine(C); Console.WriteLine(); sw.Restart(); C = CAT(1000000); sw.Stop(); Console.WriteLine("Catalan(1.000.000)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(bL(C) + " bits"); Console.WriteLine(); sw.Restart(); C = CAT(2000000); sw.Stop(); Console.WriteLine("Catalan(2.000.000)"); Console.WriteLine(sw.ElapsedMilliseconds + " ms"); Console.WriteLine(bL(C) + " bits"); Console.ReadLine(); } }
2012/02/24
Binomial Coefficient part 2, Catalan Numbers
2011/12/29
Factorial by binary splitting, part 2
using Xint = System.Numerics.BigInteger; using System; class Factorial { public static Xint F(uint n) { uint a = 0; uint s = 0; Xint P = 1; Xint Q = 1; uint b = 1; for (int i = fL2((int)(n / 2)); i >= 0; i--) { a = n >> i; s = s + a / 2; a = a - 1 | 1; P = Q * P; Q = OddP(a, b) * Q; b = a + 2; } return Q * P << (int)s; } private static Xint OddP(uint a, uint b) { if (a == b) return a; uint m = (a + b) / 2; m += m & 1; return OddP(a, m + 1) * OddP(m - 1, b); } private static int fL2(int n) { int i = -1; for (; n > 0; n /= 2) i++; return i; } // 42!=42*41*.. // =41*39*...*3*1 * 42*40*...*4*2 // =41,1? * 42,2? // 42,2?=21! * 2^(42/2) // 21!=21*19*...*3*1 * 20*18*...*4*2 // =21,1? * 20,2? // 20,2?=10! * 2^(20/2) // 10!=9*7*5*3*1 * 10*8*6*4*2 // =9,1? * 10,2? // 10,2?=5! * 2^(10/2) // 5!=5*3*1 * 4*2 // =5,1? * 4,2? // 4,2?=2! * 2^( 4/2) // 2!=1*1 * 2 // =1,1? * 2,2? // 2,2?=1! * 2^( 2/2) // // = 41,1? * 21,1? * 9,1? * 5,1? * 2^(42/2+20/2+10/2+4/2+2/2) // = 5,1? * 9,1? * 21,1? * 41,1? * 2^(21+10+5+2+1) // = 5,1?^4 * 9,7?^3 * 21,11?^2 * 41,23? * 2^39 // = a^4 * b^3 * c^2 * d^1 << 39 // // = 1 * a // = a^1 * a*b // = a^2 * b^1 * a*b*c // = a^3 * b^2 * c^1 * a*b*c*d // = a^4 * b^3 * c^2* d^1 << 39 static void Main() { Console.WriteLine(F(0)); Console.WriteLine(F(100)); Console.ReadLine(); } }
2011/10/11
Half Companion Pell Numbers
// Pell(1000000) in 199 mS (Athlon X4 640, XP, 2GB) // Half Companion Pell Numbers are defined by: // H[0] = 1, H[1] = 1, H[n+2]= 2*H[n+1] + H[n] // The sequence is 1,1,3,7,17,41,99,239,... // Small values of n use a table. // Larger values use a binary powering approach. // HCPs(n) returns H[n] and H[n+1] // The most important formula used is: // // H[2n] = 2 * H[n]^2 - (-1)^n // // Each bit requires two squares. // HCP(n) returns H[n] // If n is odd the lowest 1 bit requires a multiply. // Trailing zero bits are done by a square each. // Pells(n) returns P[n] and P[n+1] // It derives a pair of Pell numbers // from a pair of HCP numbers: // // 2*P[n] = H[n+1] - H[n] // 2*P[n+1] = H[n+1] + H[n] // // Pell(n) returns P[n] // The least significant bit requires a multiply. using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class Half_Companion_Pell_Numbers { private static int[] smallHCP = { 1, 1, 3, 7, 17, 41, 99, 239, 577, 1393, 3363, 8119, 19601, 47321, 114243, 275807, 665857 }; public static Xint[] HCPs(int n) { if (n < 16) return new Xint[2] { smallHCP[n++], smallHCP[n] }; int i = fL2(n) - 4; int m = (n >> i) / 2 & 7 + 8; Xint H0 = smallHCP[m++], H1 = smallHCP[m]; for (; i >= 0; i--) { H0 = SQ(H0); H1 = SQ(H1); switch (((3 << i) & n) >> i) { case 0: H1 = H1 - H0 + 1; H0 = 2 * H0 - 1; break; case 1: H0 = H1 - H0 + 1; H1 = 2 * H1 + 1; break; case 2: H1 = H1 - H0 - 1; H0 = 2 * H0 + 1; break; case 3: H0 = H1 - H0 - 1; H1 = 2 * H1 - 1; break; } } return new Xint[2] { H0, H1 }; } public static Xint HCP(int n) { if (n < 17) return smallHCP[n]; if ((n & 1) == 1) return HCP_n_is_odd(n); int z = fL2(n & -n); int y = fL2(n) - z - 3; if (y < 0) { z += y; Xint H = smallHCP[n >> z]; for (; z > 0; z--) H = 2 * SQ(H) - 1; return H; } else { Xint H = 2 * SQ(HCP_n_is_odd(n >> z)) + 1; for (; z > 1; z--) H = 2 * SQ(H) - 1; return H; } } private static Xint HCP_n_is_odd(int n) { Xint[] H = HCPs(n / 2); return (n & 2) - 1 + 2 * MTP(H[1], H[0]); } private static int[] smallPell = { 0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461, 80782, 195025, 470832 }; public static Xint[] Pells(int n) { if (n < 16) return new Xint[2] { smallPell[n++], smallPell[n] }; Xint[] H = HCPs(n); return new Xint[2] { (H[1] - H[0]) / 2, (H[1] + H[0]) / 2 }; } public static Xint Pell(int n) { if (n < 17) return smallPell[n]; Xint[] H = HCPs(n / 2); return (n & 1) == 0 ? MTP(H[1] - H[0], H[0]) : MTP(H[1] + H[0], H[0]) + ((n & 2) - 1); } // use the easy to copy faster versions: // http://bigintegers.blogspot.com/2011/fibonacci-numbers-part-1.html private static Xint MTP(Xint U, Xint V) { return U * V; } private static Xint SQ(Xint U) { return U * U; } private static int fL2(int n) { int i = -1; for (; n > 0; n /= 2) i++; return i; } private static Stopwatch sw = new Stopwatch(); static void Main() { Pell(1000000); sw.Restart(); for (int n = 0; n < 10; n++) Pell(1000000); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); Console.ReadLine(); } }
2011/10/06
Pell Numbers
// Pell(1000000) in 207 mS (Athlon X4 640, XP, 2GB) // Pell Numbers are defined by: // P[0] = 0, P[1] = 1, P[n+2] = 2*P[n+1] + P[n] // The sequence is 0,1,2,5,12,29,70,169,... // Formula's used are: // (P[n])^2 = P[n+1]*P[n-1] - (-1)^n // P[m+n] = P[m]*P[n+1] + P[m-1]*P[n] // Small values of n use a table. // Larger values use a binary powering approach, // which scannes n from left to right, // from most to least significant bit. // Pells(n) returns P[n] and P[n+1] // Each bit requires one square and one multiply. // Pell(n) returns P[n] // The least significant bit requires one multiply. // Result: mS(Pells(n)) / mS(Pell(n)) ~ 1.41 // ~ square root of two!? using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class Pell_Numbers { private static int[] smallPell = { 0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461, 80782, 195025, 470832 }; public static Xint[] Pells(int n) { if (n < 16) return new Xint[2] { smallPell[n++], smallPell[n] }; int i = fL2(n) - 4; int m = (n >> i) / 2 & 7 + 8; Xint[] P = { smallPell[m++], smallPell[m] }; for (; i >= 0; i--) { switch (((3 << i) & n) >> i) { case 0: P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) }; P = new Xint[2] { P[1] - P[0], P[1] + P[0] + 1 }; break; case 2: P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) }; P = new Xint[2] { P[1] - P[0], P[1] + P[0] - 1 }; break; case 1: P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) }; P = new Xint[2] { P[1] - P[0] - 1, P[1] + P[0] }; break; case 3: P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) }; P = new Xint[2] { P[1] - P[0] + 1, P[1] + P[0] }; break; } } return P; } public static Xint Pell(int n) { if (n < 17) return smallPell[n]; Xint[] P = Pells(n / 2); return (n & 1) == 0 ? 2 * MTP(P[1] - P[0], P[0]) : 2 * MTP(P[1] + P[0], P[0]) + (1 - (n & 2)); } // use the easy to copy faster versions: // http://bigintegers.blogspot.com/2011/fibonacci-numbers-part-1.html private static Xint MTP(Xint U, Xint V) { return U * V; } private static Xint SQ(Xint U) { return U * U; } private static int fL2(int n) { int i = -1; for (; n > 0; n /= 2) i++; return i; } private static Stopwatch sw = new Stopwatch(); static void Main() { for (int n = 0; n < 10; n++) { sw.Restart(); Pells(1000000); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); } Console.ReadLine(); } }
2011/09/18
Lucas Numbers
// Lucas(1000000) in 55 mS (Athlon X4 640, XP, 2GB) // http://gmplib.org/manual/Lucas-Numbers-Algorithm.html#Lucas-Numbers-Algorithm // // mpz_lucnum2_ui derives a pair of Lucas numbers from a pair // of Fibonacci numbers with the following simple formulas. // // L[k] = F[k] + 2*F[k-1] // L[k-1] = 2*F[k] - F[k-1] // // mpz_lucnum_ui is only interested in L[n], and some work can be saved. // Trailing zero bits on n can be handled with a single square each. // // L[2k] = L[k]^2 - 2*(-1)^k // // And the lowest 1 bit can be handled with one multiply of a pair of // Fibonacci numbers, similar to what mpz_fib_ui does. // // L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k // Results: // // mS(Lucs(n)) / mS(Fibs(n)) ~ 1.00 // mS(Luc(2n+1)) / mS(Fib(2n+1)) ~ 1.00 // mS(Luc(2n)) / mS(Fib(2n)) ~ 0.76 // mS(Luc(n)) / ms(Fib(n)) ~ 0.88 // mS(Lucs(n)) / ms(Luc(n)) ~ 1.55 using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class Lucas { private static int[] lucSmall = { 2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843, 1364, 2207, 3571, 5778, 9349, 15127, 24476, 39603, 64079, 103682, 167761, 271443, 439204, 710647, 1149851, 1860498, 3010349 }; public static Xint[] Lucs(int n) // { L(n), L(n + 1) } { if (n < 31) return new Xint[2] { lucSmall[n++], lucSmall[n] }; Xint[] F = Fibs(n); return new Xint[2] { 2 * F[1] - F[0], F[1] + 2 * F[0] }; } public static Xint Luc(int n) // L(n) { if (n < 32) return lucSmall[n]; if ((n & 1) == 1) return Luc_n_is_odd(n); int z = fL2(n & -n); // trailing zeros int y = fL2(n) - z - 4; if (y < 0) { z += y; Xint L = lucSmall[n >> z]; for (; z > 0; z--) L = SQ(L) - 2; return L; } else { Xint L = SQ(Luc_n_is_odd(n >> z)) + 2; for (; z > 1; z--) L = SQ(L) - 2; return L; } } private static Xint Luc_n_is_odd(int n) { n = n / 2 - 1; Xint[] F = Fibs(n); return 4 - 8 * (n & 1) + MTP(5 * F[0], 2 * F[1] + F[0]); } private static int[] fibSmall = { 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578 }; public static Xint[] Fibs(int n) // { F(n), F(n + 1) } { if (n < 33) return new Xint[2] { fibSmall[n++], fibSmall[n] }; n++; int i = fL2(n) - 4; int m = n >> i; Xint[] F = { fibSmall[m - 1], fibSmall[m] }; m = ((1 << i) & n) >> i; for (--i; i >= 0; i--) { F = new Xint[2] { SQ(F[0]), SQ(F[1]) }; F = new Xint[2] { F[1] + F[0], 2 - 4 * m + 4 * F[1] - F[0] }; m = (1 << i & n) >> i; F[1 - m] = F[1] - F[0]; } return F; } public static Xint Fib(int n) // F(n) { if (n < 34) return fibSmall[n]; n++; int i = fL2(n) - 4; int m = n >> i; Xint[] F = { fibSmall[m - 1], fibSmall[m] }; m = ((1 << i) & n) >> i; for (--i; i > 0; i--) { F = new Xint[2] { SQ(F[0]), SQ(F[1]) }; F = new Xint[2] { F[1] + F[0], 2 - 4 * m + 4 * F[1] - F[0] }; m = ((1 << i) & n) >> i; F[1 - m] = F[1] - F[0]; } return (n & 1) == 0 ? 2 - 4 * m + MTP(3 * F[1] + F[0], F[1] - F[0]) : MTP(F[1] + 2 * F[0], F[1]); } private static Stopwatch sw = new Stopwatch(); static void Main() { Luc(1000000); for (int i = 0; i < 10; i++) { sw.Restart(); Luc(1000000); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); } Console.ReadLine(); } // use the easy to copy faster versions: // http://bigintegers.blogspot.com/2011/fibonacci-numbers-part-1.html private static Xint MTP(Xint U, Xint V) { return U * V; } private static Xint SQ(Xint U) { return U * U; } private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; } }
2011/09/05
Inverse Fibonacci filter(1/30.000.000)
// Another simple test for Fibonacci numbers: // Take the repeated digit sum of a Fibonacci number. // F11=89 8+9=17 1+7=8 // The repeated digit sum for odd Fibonacci numbers never is 2 or 7. // Why? I asked Dr Math: // http://mathforum.org/dr.math/ // The repeated digit sum is the same as taking the remainder after division by 9. // Next "proof": // F%9 // F0 0 0------<---- // F1 1 1--<-- | // F2 1 1 | | // F3 2 2 | | // F4 3 3 | | // F5 5 5 | | // F6 8 8 | | // F7 13 4 | | // F8 21 3 | | // F9 34 7 | | // F10 55 1 | | // F11 89 8 | | (F(n+2))%9 = ( (F(n+1))%9 + (F(n))%9 ) % 9 // F12 144 0 | | // F13 233 8 | | // F14 377 8 | | 377%9 = 8 // F15 610 7 | | 610%9 = 7 // F16 987 6 | | 987%9 = 377%9 + 610%9 = (8+7)%9=6 // F17 1597 4 | | // F18 2584 1 | | // F19 4181 5 | | // F20 6765 6 | | // F21 10946 2 | | // F22 17711 8 | | // F23 28657 1-->-- | // F24 46368 0------>---- // Division by 8 gives 6 different remainders, // 0%8=0, 1%8=1, 1%8=1, 2%8=2, 3%8=3, 5%8=5, 8%8=0, 5%8=5, 5%8=5, // 10%8=2,7%8=7, 9%8=1, 8%8=0....6 different remainders(0,1,2,3,5,7) // Division by 11 gives 7 different remainders. // Division by 18 gives 11 different remainders. // Division by F46 gives 66 different remainders. // Which means the chance a number is a Fibonacci number // becomes about 1 to 30.000.000(~F46/66). // +-----------------------------------+ // | | different | // | devisor | remainders | // |----------------------|------------| // | 8 | F6 | 6 | // | 11 | L5 | 7 | // | 18 | L6 | 11 | // | 21 | F8 | 9 | // | 29 | L7 | 10 | 10/29 < 9/21 // | 38 | F18/68 | 13 | // | 47 | L8 | 15 | // | 48 | F12/3 | 15 | // | 55 | F10 | 12 | 12*48 < 15*55 // | ..... | .... | .. | // | ........ | ..... | .. | // | 1201881744 | L45/2 | 69 | // | 1268860318 | F56/147 | 67 | // | 1536404311 | L55/199 | 81 | // | 1602508992 | F48/3 | 69 | // | 1836311903 | F46 | 66 | // +-----------------------------------+ // The devisors found are closely related to the closely related // Fibonacci and Lucas numbers, see the table below. // F12(1,2,3) is F12/1, F12/2, F12/3 is 144, 72, 48 // // +--------------------------------------------------------------------+ // | F6(1) | F32(1,3,7,21,47) | // | F8(1) | F34(1) | // | F10(1) | F36(1,2,3,4,6,8,9,18,24,48,51) | // | F12(1,2,3) | F38(1) | // | F14(1) | F39(2) | // | F16(1,3,7) | F40(3,5,7,11,15,21,35,55) | // | F18(1,2) | F42(1,2,4,8,26) | // | F20(1,3) | F44(1,3) | // | F22(1) | F46(1,139) | // | F24(1,2,3,4,7,8,9,48) | F48(3,4,6,7,8,9,14,18,21,23,24,28,36,47, | // | F26(1) | 48,56,92,94,96,139,144,161,168,329) | // | F27(2) | F50(11,25,55) | // | F28(1,3) | F54(436,872) | // | F30(1,2,4,5,8,10,20,22) | F56(147) | // |-------------------------|------------------------------------------| // | L5(1) | L24(1) | // | L6(1) | L25(1,7,11) | // | L7(1) | L27(1,2,4) | // | L8(1) | L28(1) | // | L9(1,2) | L29(1) | // | L10(1) | L31(1) | // | L11(1) | L33(1,2,4) | // | L12(1) | L34(1) | // | L13(1) | L35(1,29,71) | // | L14(1) | L36(1) | // | L15(1) | L37(1) | // | L16(1) | L39(1,2,4) | // | L17(1) | L41(1) | // | L18(1) | L43(1) | // | L19(1) | L45(2,4,11,19,31,38,76,124,181,209) | // | L21(1,2,4) | L49(29) | // | L22(1) | L55(199) | // | L23(1) | | // +--------------------------------------------------------------------+ // // F24/48=L12*3, F32/47=L16*21, probably there are more. // For X being a (pseudo) random number with 16.000.000 bits, // invFib(X) took 78 mS, now it takes less than 6 mS. using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class inverseFibonacci { private static double c0 = Math.Log(Math.Sqrt(5)); private static double c1 = Math.Log(Math.Sqrt(5) + 1) - Math.Log(2); private const double c2 = 1.0157640563849244; private const double c3 = 35.73932166425341; private static int invFib(Xint X) { if (X <= int.MaxValue) switch ((int)X) { case 0: return 0; case 2: return 3; case 3: return 4; case 5: return 5; case 8: return 6; case 13: return 7; case 21: return 8; case 34: return 9; case 55: return 10; case 89: return 11; case 144: return 12; case 233: return 13; case 377: return 14; case 610: return 15; case 987: return 16; case 1597: return 17; case 2584: return 18; case 4181: return 19; case 6765: return 20; case 10946: return 21; case 17711: return 22; case 28657: return 23; case 46368: return 24; case 75025: return 25; case 121393: return 26; case 196418: return 27; case 317811: return 28; case 514229: return 29; case 832040: return 30; case 1346269: return 31; case 2178309: return 32; case 3524578: return 33; case 5702887: return 34; case 9227465: return 35; case 14930352: return 36; case 24157817: return 37; case 39088169: return 38; case 63245986: return 39; case 102334155: return 40; case 165580141: return 41; case 267914296: return 42; case 433494437: return 43; case 701408733: return 44; case 1134903170: return 45; case 1836311903: return 46; default: return -1; } if (X.IsEven) { int x = (int)(X & 31); if ((x != 2) && ((x & 7) != 0)) return -1; } int r = (int)Xint.Remainder(X, 1836311903); if ((r != 0) && (r != 1) && (r != 1134903170) && (r != 2) && (r != 5) && (r != 13) && (r != 34) && (r != 89) && (r != 233) && (r != 610) && (r != 1597) && (r != 4181) && (r != 10946) && (r != 28657) && (r != 75025) && (r != 196418) && (r != 514229) && (r != 1346269) && (r != 3524578) && (r != 9227465) && (r != 24157817) && (r != 63245986) && (r != 165580141) && (r != 433494437) && (r != 3) && (r != 8) && (r != 21) && (r != 55) && (r != 144) && (r != 377) && (r != 987) && (r != 2584) && (r != 6765) && (r != 17711) && (r != 46368) && (r != 121393) && (r != 317811) && (r != 832040) && (r != 2178309) && (r != 5702887) && (r != 14930352) && (r != 39088169) && (r != 102334155) && (r != 267914296) && (r != 701408733) && (r != 1568397607) && (r != 1733977748) && (r != 1797223734) && (r != 1821381551) && (r != 1830609016) && (r != 1834133594) && (r != 1835479863) && (r != 1835994092) && (r != 1836190510) && (r != 1836265535) && (r != 1836294192) && (r != 1836305138) && (r != 1836309319) && (r != 1836310916) && (r != 1836311526) && (r != 1836311759) && (r != 1836311848) && (r != 1836311882) && (r != 1836311895) && (r != 1836311900) && (r != 1836311902)) return -1; double m = (Xint.Log(X) + c0) / c1; if (m > int.MaxValue - 1) return -2; int n = (int)Math.Round(m); m = Math.Abs(n - m); if (m > 9.07086593448401E-07) return -1; if (n < 69) return m < 7.105427357602E-15 ? n : -1; if (n < 129) return m < 1.4210854715203E-14 && isFib(X, n) ? n : -1; if (n < 169) return m < 2.8421709430405E-14 && isFib(X, n) ? n : -1; if (n < 338) return m < 5.6843418860809E-14 && isFib(X, n) ? n : -1; if (n < 1026) return m < 1.13686837721617E-13 && isFib(X, n) ? n : -1; if (n < 1087) return m < 2.27373675443233E-13 && isFib(X, n) ? n : -1; if (n < 4120) return m < 4.54747350886465E-13 && isFib(X, n) ? n : -1; if (n < 7593) return m < 9.09494701772929E-13 && isFib(X, n) ? n : -1; if (n < 8522) return m < 1.81898940354587E-12 && isFib(X, n) ? n : -1; if (n < 17136) return m < 3.63797880709172E-12 && isFib(X, n) ? n : -1; if (n < 35166) return m < 7.27595761418344E-12 && isFib(X, n) ? n : -1; if (n < 131083) return m < 1.45519152283670E-11 && isFib(X, n) ? n : -1; if (n < 259737) return m < 2.91038304567338E-11 && isFib(X, n) ? n : -1; if (n < 272394) return m < 5.82076609134675E-11 && isFib(X, n) ? n : -1; if (n < 1048589) return m < 1.16415321826936E-10 && isFib(X, n) ? n : -1; if (n < 1823140) return m < 2.32830643653871E-10 && isFib(X, n) ? n : -1; if (n < 2179038) return m < 4.65661287307740E-10 && isFib(X, n) ? n : -1; if (n < 8388610) return m < 9.31322574615480E-10 && isFib(X, n) ? n : -1; if (n < 14011307) return m < 1.86264514923097E-09 && isFib(X, n) ? n : -1; if (n < 24678101) return m < 3.72529029846192E-09 && isFib(X, n) ? n : -1; return m < Math.Pow(Math.E, c2 * Math.Log(n) - c3) && isFib(X, n) ? n : -1; } private static bool isFib(Xint X, int n) { for (int i = 3; i < 34; i++) if (n % i == 0) if (Xint.Remainder(X, fibSmall[i]) > 0) return false; return X == Fib(n) ? true : false; } private static int[] fibSmall = { 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578 }; public static Xint Fib(int n) // returns F(n) { if (n < 34) return fibSmall[n]; n++; int i = fL2(n) - 4; int m = n >> i; Xint[] F = { fibSmall[m - 1], fibSmall[m] }; m = ((1 << i) & n) >> i; for (--i; i > 0; i--) { F = new Xint[2] { SQ(F[0]), SQ(F[1]) }; F = new Xint[2] { F[1] + F[0], 4 * F[1] - F[0] + (2 - 4 * m) }; m = ((1 << i) & n) >> i; F[1 - m] = F[1] - F[0]; } if ((n & 1) == 1) return MTP(F[1] + 2 * F[0], F[1]); return MTP(3 * F[1] + F[0], F[1] - F[0]) + (2 - 4 * m); } // use the easy to copy faster versions: // http://bigintegers.blogspot.com/2011/fibonacci-numbers-part-1.html private static Xint MTP(Xint U, Xint V) { return U * V; } private static Xint SQ(Xint U) { return U * U; } private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; } private static Stopwatch sw = new Stopwatch(); static void Main() { Xint X = RND(16000000); invFib(X); sw.Restart(); for (int i = 0; i < 100; i++) { if (invFib(X) != -1) Console.WriteLine(i); } sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); Console.ReadLine(); } private static int seed; public static Xint RND(int n) { if (n < 2) return n; if (seed == int.MaxValue) seed = 0; else seed++; Random rand = new Random(seed); byte[] bytes = new byte[(n + 15) / 8]; rand.NextBytes(bytes); int i = bytes.Length - 1; bytes[i] = 0; n = i * 8 - n; i--; bytes[i] >>= n; bytes[i] |= (byte)(128 >> n); return new Xint(bytes); } }
Subscribe to:
Posts (Atom)