You're right, I hadn't, but it sounds to me like it worked pretty well in your tests. It works better in Python than in Julia, apparently:

>>> [int(round(((1+sqrt(5))**n - (1-sqrt(5))**n)/(2**n*sqrt(5)))) for n in range(128)]
[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, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073L, 4807526976L, 7778742049L, 12586269025L, 20365011074L, 32951280099L, 53316291173L, 86267571272L, 139583862445L, 225851433717L, 365435296162L, 591286729879L, 956722026041L, 1548008755920L, 2504730781961L, 4052739537881L, 6557470319842L, 10610209857723L, 17167680177565L, 27777890035288L, 44945570212853L, 72723460248141L, 117669030460994L, 190392490709135L, 308061521170130L, 498454011879265L, 806515533049395L, 1304969544928660L, 2111485077978055L, 3416454622906716L, 5527939700884771L, 8944394323791488L, 14472334024676260L, 23416728348467744L, 37889062373144008L, 61305790721611752L, 99194853094755776L, 160500643816367552L, 259695496911123328L, 420196140727490880L, 679891637638614272L, 1100087778366105088L, 1779979416004719360L, 2880067194370824704L, 4660046610375544832L, 7540113804746369024L, 12200160415121913856L, 19740274219868282880L, 31940434634990198784L, 51680708854858489856L, 83621143489848688640L, 135301852344707186688L, 218922995834555891712L, 354224848179263111168L, 573147844013818970112L, 927372692193082081280L, 1500520536206901248000L, 2427893228399983329280L, 3928413764606884839424L, 6356306993006868692992L, 10284720757613753532416L, 16641027750620622225408L, 26925748508234379952128L, 43566776258855008468992L, 70492524767089384226816L, 114059301025944392695808L, 184551825793033785311232L, 298611126818978194784256L, 483162952612011980095488L, 781774079430990309097472L, 1264937032043002356301824L, 2046711111473992665399296L, 3311648143516995021701120L, 5358359254990987687100416L, 8670007398507983245672448L, 14028366653498970932772864L, 22698374052006954178445312L, 36726740705505935848636416L, 59425114757512894322049024L, 96151855463018834465652736L, 155576970220531703017897984L]

But, although it gives close answers up to cfib(604) (after which it throws an exception) it still gives slightly inexact answers starting at 308061521170129 (cfib(71)), which is between 2⁴⁸ and 2⁴⁹. I'd expect it to make it a few items further before going astray, since it's presumably doing floating-point calculations with 53 bits of mantissa that are supposed to be rounded to less than half an ulp.

I'm comparing with

>>> fibs = range(2)
>>> while len(fibs) < 128: fibs.append(fibs[-2] + fibs[-1])
...
>>> fibs
[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, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073L, 4807526976L, 7778742049L, 12586269025L, 20365011074L, 32951280099L, 53316291173L, 86267571272L, 139583862445L, 225851433717L, 365435296162L, 591286729879L, 956722026041L, 1548008755920L, 2504730781961L, 4052739537881L, 6557470319842L, 10610209857723L, 17167680177565L, 27777890035288L, 44945570212853L, 72723460248141L, 117669030460994L, 190392490709135L, 308061521170129L, 498454011879264L, 806515533049393L, 1304969544928657L, 2111485077978050L, 3416454622906707L, 5527939700884757L, 8944394323791464L, 14472334024676221L, 23416728348467685L, 37889062373143906L, 61305790721611591L, 99194853094755497L, 160500643816367088L, 259695496911122585L, 420196140727489673L, 679891637638612258L, 1100087778366101931L, 1779979416004714189L, 2880067194370816120L, 4660046610375530309L, 7540113804746346429L, 12200160415121876738L, 19740274219868223167L, 31940434634990099905L, 51680708854858323072L, 83621143489848422977L, 135301852344706746049L, 218922995834555169026L, 354224848179261915075L, 573147844013817084101L, 927372692193078999176L, 1500520536206896083277L, 2427893228399975082453L, 3928413764606871165730L, 6356306993006846248183L, 10284720757613717413913L, 16641027750620563662096L, 26925748508234281076009L, 43566776258854844738105L, 70492524767089125814114L, 114059301025943970552219L, 184551825793033096366333L, 298611126818977066918552L, 483162952612010163284885L, 781774079430987230203437L, 1264937032042997393488322L, 2046711111473984623691759L, 3311648143516982017180081L, 5358359254990966640871840L, 8670007398507948658051921L, 14028366653498915298923761L, 22698374052006863956975682L, 36726740705505779255899443L, 59425114757512643212875125L, 96151855463018422468774568L, 155576970220531065681649693L]

For exact answers, you clearly can't do better than O(N) (where N is the magnitude of the input parameter, not its length in bits) because the output is of size O(N). But for approximate answers, you can.

I agree about the uses for doubly recursive Fibonacci!