05/17/17

Exploring Fibonacci Fractions with Julia

Recently, I came across a fascinating blog and video from Mind your Decisions. It is about how a fraction
\frac{1}{999,999,999,999,999,999,999,998,,999,999,999,999,999,999,999,999}
would show the Fibonacci numbers in order when looking at its decimal output.

On a spreadsheet and most standard programming languages, such output can not be attained due to the limited precision for floating point numbers. If you try this on R or Python, you will get an output of 1e-48.
Wolfram alpha,however,allows arbitrary precision.

In Julia by default we get a little better that R and Python

julia> 1/999999999999999999999998999999999999999999999999
1.000000000000000000000001000000000000000000000002000000000000000000000003000002e-48
 
julia> typeof(ans)
BigFloat

We observe here that we are getting the first few Fibonacci numbers 1, 1, 2, 3. We need more precision to get more numbers. Julia has arbitrary precision arithmetic baked into the language. We can crank up the precision of the BigFloat type on demand. Of course, the higher the precision, the slower the computation and the greater the memory we use. We do that by setprecision.

julia> setprecision(BigFloat,10000)
10000

Reevaluating, we get

julia> 1/999999999999999999999998999999999999999999999999
1.00000000000000000000000100000000000000000000000200000000000000000000000300000000000000000000000500000000000000000000000800000000000000000000001300000000000000000000002100000000000000000000003400000000000000000000005500000000000000000000008900000000000000000000014400000000000000000000023300000000000000000000037700000000000000000000061000000000000000000000098700000000000000000000159700000000000000000000258400000000000000000000418100000000000000000000676500000000000000000001094600000000000000000001771100000000000000000002865700000000000000000004636800000000000000000007502500000000000000000012139300000000000000000019641800000000000000000031781100000000000000000051422900000000000000000083204000000000000000000134626900000000000000000217830900000000000000000352457800000000000000000570288700000000000000000922746500000000000000001493035200000000000000002415781700000000000000003908816900000000000000006324598600000000000000010233415500000000000000016558014100000000000000026791429600000000000000043349443700000000000000070140873300000000000000113490317000000000000000183631190300000000000000297121507300000000000000480752697600000000000000777874204900000000000001258626902500000000000002036501107400000000000003295128009900000000000005331629117300000000000008626757127200000000000013958386244500000000000022585143371700000000000036543529616200000000000059128672987900000000000095672202604100000000000154800875592000000000000250473078196100000000000405273953788100000000000655747031984200000000001061020985772300000000001716768017756500000000002777789003528800000000004494557021285300000000007272346024814100000000011766903046099400000000019039249070913500000000030806152117012900000000049845401187926400000000080651553304939300000000130496954492865700000000211148507797805000000000341645462290670700000000552793970088475700000000894439432379146400000001447233402467622100000002341672834846768500000003788906237314390600000006130579072161159100000009919485309475549700000016050064381636708800000025969549691112258500000042019614072748967300000067989163763861225800000110008777836610193100000177997941600471418900000288006719437081612000000466004661037553030900000754011380474634642900001220016041512187673800001974027421986822316700003194043463499009990500005168070885485832307200008362114348984842297700013530185234470674604900021892299583455516902600035422484817926191507500057314784401381708410100092737269219307899917600150052053620689608327700242789322839997508245300392841376460687116573000635630699300684624818301028472075761371741391301664102775062056366209602692574850823428107600904356677625885484473810507049252476708912581411411405930102594397055221918455182579303309636633329861112681897706691855248316295261201016328488578177407943098723020343826493703204299739348832404671111147398462369176231164814351698201718008635835925499096664087184867000739850794865805193502836665349891529892378369837405200686395697571872674070550577925589950242511475751264321287522115185546301842246877472357697022053e-48

That is looking much better. However it we be nice if we could extract the Fibonacci numbers that are buried in that long decimal. Using the approach in the original blog. We define a function

y(x)=one(x)-x-x^2

and calculate the decimal

a=1/y(big"1e-24")

Here we use the non-standard string literal big"..." to insure proper interpretation of our input. Using BigFloat(1e-24)) would first construct at floating point with limited precision and then do the conversion. The initial loss of precision will not be recovered in the conversion, and hence the use of big. Now we extract our Fibonacci numbers by this function

function extract_fib(a)
   x=string(a)
   l=2
   fi=BigInt[]
   push!(fi,1)
   for i=1:div(length(x)-24,24)
        j=parse(BigInt,x[l+1+(i-1)*24:l+i*24])
        push!(fi,j)
   end
   fi
end

Here we first convert our very long decimal number of a string and they we exploit the fact the Fibonacci numbers occur in blocks that 24 digits in length. We get out output in an array of BigInt. We want to compare the output with exact Fibonacci numbers, we just do a quick and non-recursive implementation.

function fib(n)
    f=Vector{typeof(n)}(n+1)
    f[1]=f[2]=1;
    for i=3:n+1
       f[i]=f[i-1]+f[i-2]
    end
    f
end

Now we compare…

fib_exact=fib(200);
fib_frac=extract_fib(a);
for i in eachindex(fib_frac)
     println(fib_exact[i], " ", fib_exact[i]-fib_frac[i])
end

We get a long sequence, we just focused here on when the discrepancy happens.

...
184551825793033096366333 0
298611126818977066918552 0
483162952612010163284885 0
781774079430987230203437 -1
1264937032042997393488322 999999999999999999999998
2046711111473984623691759 1999999999999999999999997
...

The output shows that just before the extracted Fibonacci number exceeds 24 digits, a discrepancy occurs. I am not quite sure why, but this was a fun exploration. Julia allows me to do mathematical explorations that would take one or even two orders of magnitude of effort to do in any other language.

01/13/17

Kaperkar’s Constant

I was recently introduced to Kaperkar’s Constant.
It is quite magical. You take any four digit number A, sort the digits from highest to lowest to create a new number A^{\text{high}}, sort the digits from lowest to highest to get A^{\text{low}}, and calculate and new number A= A^{\text{high}}- A^{\text{low}}. You repeat this procedure enough times and you end up with A=6174.

I made a nifty implementation of that in Julia below.

05/2/15

Monkeys and Coconuts

Here is my attempt to solve the monkeys and coconuts problem. The story goes
that five sailors were stranded on an island and had decided to gather some
coconuts for their provisions. They put all the coconuts in one large pile and
went to sleep. One sailor got up and fearing that there could be problems when the
time to came to divide the pile, he divided the pile five ways and noticing that
he has an extra coconut, he gave it to a monkey, and then hid his stash. The other
four sailor repeated the same procedure. When they woke up they noticed that
they had a smaller pile and proceeded to divide into five equal piles, this time
around there were no extra coconut left for the monkey. So the question is: what was
the size of the original pile?

We will denote by x_i as the size of the
pile after the ith sailor has carried out his procedure. In this system x_0
is original size of the pile. So following this procedure we then proceed as
follows:

x_1=\frac{4(x_0-1)}{5}
\cdots
x_i=\frac{4(x_{i-1}-1)}{5}
\cdots
x_5=\frac{4(x_4-1)}{5}

It is important to note that x_i \in \mathbb{N}_0 for i=0\ldots5, also \frac{x_5}{5}\in \mathbb{N}_0. Alternatively, we think of the reverse procedure and express the above as
x_4=\frac{5x_5}{4}+1
\cdots
x_i=\frac{5x_{i+1}}{4}+1
\cdots
x_0=\frac{5x_1}{4}+1

Observing that x_5 has to be divisible by 4 and 5 (last equation in the first system and first equation in the second), one can brute force the solution(s) by the following Julia code:

Which results in

[2496,1996,1596,1276,1020]
m=51, x₀= 3121
[14996,11996,9596,7676,6140]
m=307, x₀= 18746
[27496,21996,17596,14076,11260]
m=563, x₀= 34371
[39996,31996,25596,20476,16380]
m=819, x₀= 49996

This corresponds nicely to the answers that were obtained by rigorous derivation in the video, however it shows how programming can easily find such solutions by brute force.

If one would like to avoid the negative or blue concocts suggested in the video and also preserve the monkey. Below is an alternative derivation. Working through the first system, one gets:

x_{5}={{4\,\left({{4\,\left({{4\,\left({{4\,\left({{4\,\left(x_{0}-  1\right)}\over{5}}-1\right)}\over{5}}-1\right)}\over{5}}-1\right)  }\over{5}}-1\right)}\over{5}}=\frac{4^5x_0}{5^5}-\frac{8404}{5^5}

Hence,
5^5x_5=4^5x_0-8404.
Realizing the x_5 has to be necessarily divisible by 5, we denote the final share that each sailor gets in the last division by s. So our Diophantine equation becomes
5^6s=4^5x_0-8404.
It will have solutions at x_0=3121, 18746, 34371 \ldots= 3121+n5^6 \text{ for } n \in \mathbb{N}_0 .

03/3/10

Mathematical Distractions

While attending a conference and hanging out with a college later in the day he recounted to me a nice mathematical puzzle that he learned while growing up in Algeria.
Puzzle 1:
Three brothers were in a bind, they could not divide their father’s inheritance. Their father has left them with 17 camels and instructions that half of his camels should go to eldest son, a third to the second eldest, and a ninth to the youngest. Knowing that there is no point in having a fraction of camel, they wanted a solution that would be fair and yet that would result in no fractional division of camels. A passerby, seeing their predicament proposed to help.
In order to resolve their problem he donated his camel to the pool of 17 camels, bring thus the total to 18. One brother will hence get 9 camels, the other would get 6, and the youngest would get 2 camels. But 9+6+2=17, therefore the passerby concluded that he would walk away with the left over camel. Problem solved… but how?
Puzzle 2:
Two farmers decided to sit down and share their lunch food. One farmer had on him five loafs of bread and the other had seven loafs. Just as they were about to start, a passerby asked if he could join them. The three sat down and they equally divided the twelve loafs amongst themselves. When they were done, the passerby give them twelve dirhams for the meal. The farmer who contributed five loafs proposed that a fair division of that amount would be for him to take five dirhams and for the other farmer to take the remaining seven. A disagreement ensued. They sought the assistance of judge to resolve their disupte and his resolution was that the farmer who contributed five loafs should only walk away with three dirhams, while the other should have the remaining nine dirhams. How is this fair?