// X=Fib(n) and n=invFib(X) have the same speed, // If (pseudo) random numbers with the bitLength of X are used, // invFib will be very fast. // X=Fib(23.000.000) takes 7540 mS (X has 15.967.636 bits ~ 16M bits) // invFib(X) takes 7640 mS, random numbers // with X's bitLength take 78 mS. // Another way to test Fibonacci numbers is: // X is a Fibo if 5*X^2 + 4 or 5*X^2 - 4 is a square number. // squaring a 16M bits number takes 8460 mS, // taking the square root takes 25000 mS, // Let's guess it's possible to take the other square root fast // then it seems to be about 400 times slower when X is not a Fibo, // 4 times slower when X is (still n will be unknown). // Did a few tests, // Fib(500.000.000) takes 10 minutes (347.120.956 bits) // Fib(600.000.000) crashes (2GB memory?) // so n_max is somewhere in between 500M and 600M // and in the line: // if (m > 9.07086593448401E-07) return -1; // 9.07086593448401E-07 should be: // Math.Pow(Math.E, c2 * Math.Log(500000000) - c3) // and for my machine I might change: // if (m > int.MaxValue - 1) return -2; // into: // if (m > 500000000) return -2; using Xint = System.Numerics.BigInteger; using System.Threading.Tasks; using System.Diagnostics; using System; class inverseFibonacci { 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 }; 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; // 1? private const double c3 = 35.73932166425341; private static int invFib(Xint X) { if (X <= fibSmall[33]) { int x = (int)X; if (x < 2) return x == 0 ? 0 : -1; int n = 3; for (int m = 16; m > 0; ) { for (; (m > 0) & (x > fibSmall[n]); m /= 2) n += m; for (; (m > 0) & (x < fibSmall[n]); m /= 2) n -= m; if (x == fibSmall[n]) return n; } return -1; } else { 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++) // Fnk is a multiple of Fn if (n % i == 0) if (Xint.Remainder(X, fibSmall[i]) > 0) return false; return X == Fib(n) ? true : false; } private static Xint Fib(int 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! 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() { sw.Restart(); Xint X = Fib(1000000); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); sw.Restart(); int n = invFib(X); sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds); Console.ReadLine(); } }
C# , System.Numerics, Multiplication, Karatsuba, Toom Cook, Division, Burnikel, Ziegler, Factorial, Luschny, Square Root, Zimmermann, Choose, Binomial Coefficient, Permutation, Combination, Eratosthenes, Primes, Fibonacci, Lucas, Pell, Catalan, Fast Random Number Generator, Overton
2011/08/06
Inverse Fibonacci, code
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment