On Tuesday, March 24, 2015 at 10:24:59 PM UTC+8, Ian wrote: > On Tue, Mar 24, 2015 at 12:20 AM, Rustom Mody <rustompm...@gmail.com> wrote: > > On Tuesday, March 24, 2015 at 10:51:11 AM UTC+5:30, Ian wrote: > >> Iteration in log space. On my desktop, this calculates fib(1000) in > >> about 9 us, fib(100000) in about 5 ms, and fib(10000000) in about 7 > >> seconds. > >> > >> def fib(n): > >> assert n >= 0 > >> if n == 0: > >> return 0 > >> > >> a = b = x = 1 > >> c = y = 0 > >> n -= 1 > >> > >> while True: > >> n, r = divmod(n, 2) > >> if r == 1: > >> x, y = x*a + y*b, x*b + y*c > >> if n == 0: > >> return x > >> b, c = a*b + b*c, b*b + c*c > >> a = b + c > > > > This is rather arcane! > > What are the identities used above? > > It's essentially the same matrix recurrence that Gregory Ewing's > solution uses, but without using numpy (which doesn't support > arbitrary precision AFAIK) and with a couple of optimizations. > > The Fibonacci recurrence can be expressed using linear algebra as: > > F_1 = [ 1 0 ] > > T = [ 1 1 ] > [ 1 0 ] > > F_(n+1) = F_n * T > > I.e., given that F_n is a vector containing fib(n) and fib(n-1), > multiplying by the transition matrix T results in a new vector > containing fib(n+1) and fib(n). Therefore: > > F_n = F_1 * T ** (n-1) > > The code above evaluates this expression by multiplying F_1 by powers > of two of T until n-1 is reached. x and y are the two elements of the > result vector, which at the end of the loop are fib(n) and fib(n-1). > a, b, and c are the three elements of the (symmetric) transition > matrix T ** p, where p is the current power of two. > > The last two lines of the loop updating a, b, and c could equivalently > be written as: > > a, b, c = a*a + b*b, a*b + b*c, b*b + c*c > > A little bit of algebra shows that if a = b + c before the assignment, > the equality is maintained after the assignment (in fact the elements > of T ** n are fib(n+1), fib(n), and fib(n-1)), so the two > multiplications needed to update a can be optimized away in favor of a > single addition.
Well, solving a homogeneous difference equation of 2 degrees and generating the solution sequence for a particular one like the Finbnaci series is a good programming practice. A more general programming practice is to generate the solution series of an arbitrary homogeneous difference euqation of integer coefficients when a real or complex solution does exist. -- https://mail.python.org/mailman/listinfo/python-list
