// Fibonacci(1000000) in 105 mS (Athlon X4 640, XP, 2GB) // Two formulas: // // F[2n-1] = F[n-1]^2 + F[n]^2 // // F[2n] = (2*F[n-1] + F[n]) * F[n] // // http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html?n=55 // Fibs(66) // 66 = 1000010 // 10000.. = 16 F0 = fibSmall(16) = 987 = F(16) // F1 = fibSmall(17) = 1597 = F(17) // F0 = (2*F1-F0)*F0 = (2*1597-987)*987 = 2178309 = F(32) // F1 = F1*F1+F0*F0 = 1597*1597+987*987 = 3524578 = F(33) // .....1. F0 = F1 = 3524578 = F(33) // F1 = F1+F0 = 3524578+2178309 = 5702887 = F(34) // F0 = (2*F1-F0)*F0 = 27777890035288 = F(66) // F1 = F1*F1+F0*F0 = 44945570212853 = F(67) // ......0 no change // // Fibs(66) = { Fib(66), Fib(67) } = { 27777890035288, 44945570212853 } // +--------------------------------+ // | Fibonacci.exe in milliSeconds | // |--------------------------------| // | n | F(n) | F(n+1)|Fibs(n)| // |--------------------------------| // | 156250 | 10 | 13 | 15 | // | 312500 | 23 | 28 | 37 | // | 625000 | 57 | 73 | 93 | // | 1250000 | 148 | 198 | 255 | // | 2500000 | 406 | 515 | 676 | // | 5000000 | 1110 | 1420 | 1860 | // |10000000 | 3010 | 3800 | 4930 | // |20000000 | 8080 | 10400 | 13500 | // +--------------------------------+ // !!! Project >> Add Reference >> System.Numerics !!! using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class program { private static int fibSmall(int n) { int[] f = { 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 }; return f[n]; } public static Xint Fib(int n) { if (n < 34) return fibSmall(n); int i = fL2(n) - 4; int m = n >> i; Xint[] F = { fibSmall(m++), fibSmall(m) }; for (--i; i > 0; i--) { F = new Xint[2] { MTP(2 * F[1] - F[0], F[0]), SQ(F[1]) + SQ(F[0]) }; if (((1 << i) & n) != 0) F = new Xint[2] { F[1], F[1] + F[0] }; } if ((n & 1) == 0) return MTP(2 * F[1] - F[0], F[0]); return SQ(F[1]) + SQ(F[0]); } public static Xint[] Fibs(int n) { if (n < 33) return new Xint[2] { fibSmall(n++), fibSmall(n) }; int i = fL2(n) - 4; int m = n >> i; Xint[] F = { fibSmall(m++), fibSmall(m) }; for (--i; i >= 0; i--) { F = new Xint[2] { MTP(2 * F[1] - F[0], F[0]), SQ(F[1]) + SQ(F[0]) }; if (((1 << i) & n) != 0) F = new Xint[2] { F[1], F[1] + F[0] }; } return F; } private static Stopwatch sw = new Stopwatch(); static void Main() { for (int i = 0; i < 20; i++) { sw.Restart(); Fib(1000000); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); } Console.ReadLine(); } 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 Xint SQ(Xint U) { return SQ(U, U.Sign * U.ToByteArray().Length << 3); } private static Xint SQ(Xint U, int n) { if (n <= 700) return U * U; if (n <= 3000) return Xint.Pow(U, 2); if (n <= 6000) return SQ2(U, n); if (n <= 10000) return SQ3(U, n); if (n <= 40000) return SQ4(U, n); return SQ2P(U, n); } private static Xint SQr(Xint U, int n) { if (n <= 3000) return Xint.Pow(U, 2); if (n <= 6000) return SQ2(U, n); if (n <= 10000) return SQ3(U, n); return SQ4(U, n); } private static Xint SQ2(Xint U1, int n) { n >>= 1; Xint U0 = U1 & ((Xint.One << n) - 1); U1 >>= n; Xint P0 = SQr(U0, n); Xint P2 = SQr(U1, n); return ((P2 << n) + (SQr(U0 + U1, n) - (P0 + P2)) << n) + P0; } private static Xint SQ3(Xint U2, int n) { n = (int)((long)(n) * 0x55555556 >> 32); Xint Mask = (Xint.One << n) - 1; Xint U0 = U2 & Mask; U2 >>= n; Xint U1 = U2 & Mask; U2 >>= n; Xint W0 = SQr(U0, n); Xint W4 = SQr(U2, n); Xint P3 = SQr((((U2 << 1) + U1) << 1) + U0, n); U2 += U0; Xint P2 = SQr(U2 - U1, n); Xint P1 = SQr(U2 + U1, 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 SQ4(Xint U3, 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 W0 = SQr(U0, n); // 0 U0 += U2; U1 += U3; Xint P1 = SQr(U0 + U1, n); // 1 Xint P2 = SQr(U0 - U1, n); // -1 U0 += 3 * U2; U1 += 3 * U3; Xint P3 = SQr(U0 + (U1 << 1), n); // 2 Xint P4 = SQr(U0 - (U1 << 1), n); // -2 Xint P5 = SQr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), n); // 4 Xint W6 = SQr(U3, 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 SQ2P(Xint A, int n) { n >>= 1; Xint[] U = new Xint[3]; U[0] = A & ((Xint.One << n) - 1); A >>= n; U[2] = A; U[1] = U[0] + A; Xint[] P = new Xint[3]; Parallel.For(0, 3, (int i) => P[i] = SQr(U[i], n)); return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0]; } public static int bL(Xint U) { byte[] bytes = (U.Sign * U).ToByteArray(); int i = bytes.Length - 1; return (i << 3) + 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 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; } }
2011/07/05
Fibonacci numbers part 1
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment