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.

03/27/17

Using pipes while running external programs in Julia

Recently I was using Julia to run ffprobe to get the length of a video file. The trouble was the ffprobe was dumping its output to stderr and I wanted to take that output and run it through grep. From a bash shell one would typically run:

ffprobe somefile.mkv 2>&1 |grep Duration

This would result in an output like

 Duration: 00:04:44.94, start: 0.000000, bitrate: 128 kb/s

This works because we used 2>&1 to redirect stderr to stdout which would in be piped to grep.

If you were try to run this in Julia

julia> run(`ffprobe somefile.mkv 2>&1 |grep Duration`)

you will get errors. Julia does not like pipes | inside the backticks command (for very sensible reasons). Instead you should be using Julia’s pipeline command. Also the redirection 2>&1 will not work. So instead, the best thing to use is and instance of Pipe. This was not in the manual. I stumbled upon it in an issue discussion on GitHub. So a good why to do what I am after is to run.

julia> p=Pipe()
Pipe(uninit => uninit, 0 bytes waiting)
 
julia> run(pipeline(`ffprobe -i  somefile.mkv`,stderr=p))

This would create a pipe object p that is then used to capture stderr after the execution of the command. Next we need to close the input end of the pipe.

julia> close(p.in)

Finally we can use the pipe with grep to filter the output.

julia> readstring(pipeline(p,`grep Duration`))
"  Duration: 00:04:44.94, start: 0.000000, bitrate: 128 kb/s\n"

We can then do a little regex magic to get the duration we are after.

julia> matchall(r"(\d{2}:\d{2}:\d{2}.\d{2})",ans)[1]
"00:04:44.94"
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.

09/25/16

Julia calling C: A minimal example

This blog is a “Hello World” example of Julia calling C.

We start of by at bit of C code we want to call from Julia. We write the following in calc_mean.c

double mean(double a, double b) {
  return (a+b) / 2;
}

To build the library, we need to create a Makefile

CC=gcc 
 
CFLAGS=-c -Wall -fPIC
 
SOURCES=calc_mean.c 
OBJECTS=$(SOURCES:.c=.o)
 
.c.o:
    $(CC) $(CFLAGS) $< -o $@ 
 
lib: $(OBJECTS)
    $(CC) -shared -fPIC -o libmean.so $(OBJECTS)
 
clean:
    rm *.o *.so

The option fPIC and -shared are essential for Julia to be able to resolve the function in our library. Now we are almost ready to build our library. From the bash terminal we invoke:

make lib

This will generate a libmean.so file.

In Julia we call the function in our c library by

x=ccall((:mean,"libmean"),Float64,(Float64,Float64),2.0,5.0)
println(x)
3.5

For this to work,

  • Julia must be running either on the same path where libmean.so resides,
  • the path to libmean.so is in LD_LIBRARY_PATH, or
  • the path to the library is pushed to Libdl.DL_LOAD_PATH via

push!(Libdl.DL_LOAD_PATH,"path_to_libmean.so")

P.S. Thanks to Christopher Rackauckas for tips on Julia highlighting.

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 .