rationalize() for BigFloats has an upper limit?

Question

I have the following code:

function recursion(i::BigFloat)
    r = BigFloat(0)
    if i >= 1.0
        i = i - 1.0
        r = 1.0/(2.0+recursion(i))
    end
    return r
end

function main()
    solution = 0

    i = BigFloat(1)
    while i < 1000
        sum = 1 + recursion(i)
        print("\nsum: ",rationalize(sum))
        if length(digits(numerator(rationalize(sum)))) > length(digits(denominator(rationalize(sum))))
            solution = solution + 1
            print("     <------- found")
        end
        i = i + 1.0
    end
    
    return solution
end

solution = main()

The goal is to find expansions of the sqrt(two) using the following equation:

(image taken from "Square Root of 2" wikipedia page)

where the rational representation has a numerator with a number of digits greater than the denominator's number of digits.

My code's logic seems to be working, however there is an upper limit where it seems to be breaking the rationalize() function when using a BigFloat type.

sum: 3//2
sum: 7//5
sum: 17//12
sum: 41//29
sum: 99//70
sum: 239//169
sum: 577//408
sum: 1393//985     <------- found
sum: 3363//2378
sum: 8119//5741
sum: 19601//13860
sum: 47321//33461
sum: 114243//80782     <------- found
sum: 275807//195025
sum: 665857//470832
sum: 1607521//1136689
sum: 3880899//2744210
sum: 9369319//6625109
sum: 22619537//15994428
sum: 54608393//38613965
sum: 131836323//93222358     <------- found
sum: 318281039//225058681
sum: 768398401//543339720
sum: 1855077841//1311738121
sum: 4478554083//3166815962
sum: 10812186007//7645370045     <------- found
sum: 26102926097//18457556052
sum: 63018038201//44560482149
sum: 152139002499//107578520350
sum: 367296043199//259717522849
sum: 886731088897//627013566048
sum: 2140758220993//1513744654945
sum: 5168247530883//3654502875938
sum: 12477253282759//8822750406821     <------- found
sum: 30122754096401//21300003689580
sum: 72722761475561//51422757785981
sum: 175568277047523//124145519261542
sum: 423859315570607//299713796309065
sum: 1023286908188737//723573111879672     <------- found
sum: 2470433131948081//1746860020068409
sum: 5964153172084899//4217293152016490
sum: 14398739476117879//10181446324101389
sum: 34761632124320657//24580185800219268
sum: 83922003724759193//59341817924539925
sum: 202605639573839043//143263821649299118
sum: 489133282872437279//345869461223138161
sum: 1180872205318713601//835002744095575440     <------- found
sum: 2850877693509864481//2015874949414289041
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522
sum: 6882627592338442563//4866752642924153522

It gets "stuck" at the fraction 6882627592338442563//4866752642924153522. Presumably because there is some sort of upper limit within the rationalize function.

Why is it getting stuck here and how can I circumvent the problem?

Answer 1

There are a few issues here, all dealing with making sure that you're using types that can get as accurate or as large as you want.

First, you need to break the habit of putting decimals in your exact numbers — as in 1.0 and 2.0, which are interpreted as literal Float64s. That's fine when you're only ever dealing with 64-bit floats, but you don't want to do algebra with a Float64 and a BigFloat. In this case, since your floats are exact integers, just use the integers 1 and 2, and Julia will promote properly to BigFloat. For things that aren't exact integers, you can explicitly give literal Rationals, and the compiler will convert them to the appropriate type efficiently. Or you can explicitly convert to the number type you need, as in sqrt(big(2)) — though usually it's better to detect the input type and then do something like sqrt(T(2)). Unlike in — say, C++ — if you enter 1/2, you'll get 0.5 as a Float64. On the other hand, 1/(2+recursion(i)) will first convert 2+recursion(i) to the appropriate BigFloat, and then take the true inverse.

Second, if you look at the docstring for rationalize, you see that it has an optional argument of the type, which defaults to Int. This means the result must have numerator and denominator of this type. But Int64 can only represent numbers up to typemax(Int64), which is 9223372036854775807. And your "stuck" numerator is ~75% as large, so you probably just can't get much larger. Instead of using the default, just call it as rationalize(BigInt, sum), and you'll be able to get much larger numbers on top and bottom.

Third, note the third argument of rationalize, which says that it will stop trying to find a better approximation once the result is within tol, which defaults to the machine precision of the input. For sqrt(big(2)), that's about 1.7e-77. So rationalize won't give you anything more once the result is within that distance of sqrt(big(2)). So even with Big*, you're results will get "stuck" eventually.

Anyway, putting that all together:

function recursion(i::BigFloat)
    r = BigFloat(0)
    if i >= 1
        i = i - 1
        r = 1/(2+recursion(i))
    end
    return r
end

function main()
    solution = 0

    i = BigFloat(1)
    while i < 1000
        sum = 1 + recursion(i)
        ratio = rationalize(BigInt, sum)
        print("\nsum: ",ratio)
        if length(digits(numerator(ratio))) > length(digits(denominator(ratio)))
            solution = solution + 1
            print("     <------- found")
        end
        i = i + 1
    end
    
    return solution
end

solution = main()

And the result looks like

sum: 3//2
sum: 7//5
sum: 17//12
sum: 41//29
sum: 99//70
sum: 239//169
sum: 577//408
sum: 1393//985     <------- found
sum: 3363//2378
sum: 8119//5741
sum: 19601//13860
sum: 47321//33461
sum: 114243//80782     <------- found
sum: 275807//195025
sum: 665857//470832
sum: 1607521//1136689
sum: 3880899//2744210
sum: 9369319//6625109
sum: 22619537//15994428
sum: 54608393//38613965
sum: 131836323//93222358     <------- found
sum: 318281039//225058681
sum: 768398401//543339720
sum: 1855077841//1311738121
sum: 4478554083//3166815962
sum: 10812186007//7645370045     <------- found
sum: 26102926097//18457556052
sum: 63018038201//44560482149
sum: 152139002499//107578520350
sum: 367296043199//259717522849
sum: 886731088897//627013566048
sum: 2140758220993//1513744654945
sum: 5168247530883//3654502875938
sum: 12477253282759//8822750406821     <------- found
sum: 30122754096401//21300003689580
sum: 72722761475561//51422757785981
sum: 175568277047523//124145519261542
sum: 423859315570607//299713796309065
sum: 1023286908188737//723573111879672     <------- found
sum: 2470433131948081//1746860020068409
sum: 5964153172084899//4217293152016490
sum: 14398739476117879//10181446324101389
sum: 34761632124320657//24580185800219268
sum: 83922003724759193//59341817924539925
sum: 202605639573839043//143263821649299118
sum: 489133282872437279//345869461223138161
sum: 1180872205318713601//835002744095575440     <------- found
sum: 2850877693509864481//2015874949414289041
sum: 6882627592338442563//4866752642924153522
sum: 16616132878186749607//11749380235262596085
sum: 40114893348711941777//28365513113449345692
sum: 96845919575610633161//68480406462161287469
sum: 233806732499933208099//165326326037771920630
sum: 564459384575477049359//399133058537705128729
sum: 1362725501650887306817//963592443113182178088     <------- found
sum: 3289910387877251662993//2326317944764069484905
sum: 7942546277405390632803//5616228332641321147898
sum: 19175002942688032928599//13558774610046711780701
sum: 46292552162781456490001//32733777552734744709300
sum: 111760107268250945908601//79026329715516201199301     <------- found
sum: 269812766699283348307203//190786436983767147107902
sum: 651385640666817642523007//460599203683050495415105
sum: 1572584048032918633353217//1111984844349868137938112
sum: 3796553736732654909229441//2684568892382786771291329
sum: 9165691521498228451812099//6481122629115441680520770
sum: 22127936779729111812853639//15646814150613670132332869
sum: 53421565080956452077519377//37774750930342781945186508
sum: 128971066941642015967892393//91196316011299234022705885     <------- found
sum: 311363698964240484013304163//220167382952941249990598278
sum: 751698464870122983994500719//531531081917181734003902441
sum: 1814760628704486452002305601//1283229546787304717998403160
sum: 4381219722279095887999111921//3097990175491791170000708761
sum: 10577200073262678228000529443//7479209897770887057999820682     <------- found
sum: 25535619868804452344000170807//18056409971033565286000350125
sum: 61648439810871582916000871057//43592029839838017630000520932
sum: 148832499490547618176001912921//105240469650709600546001391989
sum: 359313438791966819268004696899//254072969141257218722003304910
sum: 867459377074481256712011306719//613386407933224037990008001809
sum: 2094232192940929332692027310337//1480845785007705294702019308528
sum: 5055923762956339922096065927393//3575077977948634627394046618865
sum: 12206079718853609176884159165123//8631001740904974549490112546258     <------- found
sum: 29468083200663558275864384257639//20837081459758583726374271711381
sum: 71142246120180725728612927680401//50305164660422142002238655969020
sum: 171752575441025009733090239618441//121447410780602867730851583649421
sum: 414647397002230745194793406917283//293199986221627877463941823267862
sum: 1001047369445486500122677053453007//707847383223858622658735230185145     <------- found
sum: 2416742135893203745440147513823297//1708894752669345122781412283638152
sum: 5834531641231893991002972081099601//4125636888562548868221559797461449
sum: 14085805418356991727446091676022499//9960168529794442859224531878561050     <------- found
sum: 34006142477945877445895155433144599//24045973948151434586670623554583549
sum: 82098090374248746619236402542311697//58052116426097312032565778987728148
sum: 198202323226443370684367960517767993//140150206800346058651802181530039845
sum: 478502736827135487987972323577847683//338352530026789429336170142047807838
sum: 1155207796880714346660312607673463359//816855266853924917324142465625655521     <------- found
sum: 2788918330588564181308597538924774401//1972063063734639263984455073299118880
sum: 6733044458057842709277507685523012161//4760981394323203445293052612223893281
sum: 16255007246704249599863612909970798723//11494025852381046154570560297746905442
sum: 39243058951466341909004733505464609607//27749033099085295754434173207717704165
sum: 94741125149636933417873079920900017937//66992092050551637663438906713182313772
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 323466434400377142162623973268164663418//228725309250740208744750893347264645481
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709
sum: 228725309250740208744750893347264645481//161733217200188571081311986634082331709

So it does get "stuck", but just evaluate that final result compared to sqrt(big(2)), and you should be pretty happy.

Answer 2

Mike asked in another answer what the convergence of this algorithm looks like when precision is set to a high value.

The answer is that it looks just about perfectly steady. The residual decreases by roughly a factor of sqrt(2) each iteration. Here is the log of the residual versus iteration: