We have 4 fishes in each tank. We have 4 tanks in sector Gamma. We have 4 sharks in sector Gamma.

I immediately thought it would be cool to whip a quick Julia code to get the full enumeration of all left truncatable primes, count the number of branches and also get the largest left truncatable prime.

using Primes function get_left_primes(s::String) p_arr=Array{String,1}() for i=1:9 number_s="$i$s" if isprime(parse(BigInt, number_s)) push!(p_arr,number_s) end end p_arr end function get_all_left_primes(l) r_l= Array{String,1}() n_end_points=0 for i in l new_l=get_left_primes(i) isempty(new_l) && (n_end_points+=1) append!(r_l,new_l) next_new_l,new_n=get_all_left_primes(new_l) n_end_points+=new_n # counting the chains append!(r_l,next_new_l) end r_l, n end

The first function just prepends a number (expressed in String for convenience) and checks for it possible primes that can emerge from a single digit prepending. For example:

julia> get_left_primes("17") 2-element Array{String,1}: "317" "617"

The second function, just makes extensive use of the first to get all left truncatable primes and also count the number of branches.

julia> all_left_primes, n_branches=get_all_left_primes([""]) (String["2", "3", "5", "7", "13", "23", "43", "53", "73", "83" … "6435616333396997", "6633396997", "76633396997", "963396997", "16396997", "96396997", "616396997", "916396997", "396396997", "4396396997"], 1442) julia> n_branches 1442 julia> all_left_primes 4260-element Array{String,1}: "2" "3" "5" "7" "13" "23" ⋮ "96396997" "616396997" "916396997" "396396997" "4396396997"

So we the full list of possible left truncatable primes with a length 4260. Also the total number of branches came to 1442.

We now get the largest left truncatable primes with the following one liner:

julia> largest_left_prime=length.(all_left_primes)|>indmax|> x->all_left_primes[x] "357686312646216567629137"

After this fun exploration, I found an implementation in Julia for just getting the largest left truncatable prime for any base in Rosseta Code.

]]>`Julia`

has good documentation on dealing with Dates and Time, however that is often in the context constructing and Date and Time objects. In this post, I am focus on the ability to iterate over Dates and Times. This is very useful in countless application.
We start of by capturing this moment and moving ahead into the future

julia> this_moment=now() 2018-04-01T23:13:33.437

In one hour that will be

julia> this_moment+Dates.Hour(1) 2018-04-02T00:13:33.437

Notice that `Julia`

was clever enough properly interpret that we will be on the in another day after exactly one hour. Thanks to it multiple dispatch of the `DateTime`

type to be able to do TimeType period arithmatic.

You can then write a nice `for`

loop that does something every four hours for the next two days.

julia> for t=this_moment:Dates.Hour(4):this_moment+Dates.Day(2) println(t) #or somethings special with that time end 2018-04-01T23:13:33.437 2018-04-02T03:13:33.437 2018-04-02T07:13:33.437 2018-04-02T11:13:33.437 2018-04-02T15:13:33.437 2018-04-02T19:13:33.437 2018-04-02T23:13:33.437 2018-04-03T03:13:33.437 2018-04-03T07:13:33.437 2018-04-03T11:13:33.437 2018-04-03T15:13:33.437 2018-04-03T19:13:33.437 2018-04-03T23:13:33.437

Often we are not so interested in the full dates. For example if we are reading a video file and we want to get a frame every 5 seconds while using `VideoIO.jl`

. We can deal here with the simpler `Time`

type.

julia> video_start=Dates.Time(0,5,20) 00:05:20

Here we are interested in starting 5 minutes and 20 seconds into the video.

Now we can make a nice loop from the start to finish

for t=video_start:Dates.Second(5):video_start+Dates.Hour(2) h=Dates.Hour(t).value m=Dates.Minute(t).value s=Dates.Second(t).value ms=Dates.Millisecond(t).value # Do something interesting with ffmpeg seek on the video end

`Julia`

is significantly faster than `c`

. I got this insightful response
Random number generation is not super cheap, so you're essentially benchmarking the performance for the built-in RNG.

— Keno Fischer (@KenoFischer) February 8, 2018

So I decided to dig deeper. Basically the standard `c`

`rand()`

is not that good. So instead I searched for the fastest Mersenne Twister there is. I downloaded the latest code and compiled it in the fastest way for my architecture.

/* eurler2.c */ #include <stdio.h> /* printf, NULL */ #include <stdlib.h> /* srand, rand */ #include "SFMT.h" /* fast Mersenne Twister */ sfmt_t sfmt; double r2() { return sfmt_genrand_res53(&sfmt); } double euler(long int n) { long int m=0; long int i; for(i=0; i<n; i++){ double the_sum=0; while(1) { m++; the_sum+=r2(); if(the_sum>1.0) break; } } return (double)m/(double)n; } int main () { sfmt_init_gen_rand(&sfmt,123456); printf ("Euler : %2.5f\n", euler(1000000000)); return 0; }

I had to compile with a whole bunch of flags which I induced from `SFMT`

‘s `Makefile`

to get faster performance.

gcc -O3 -finline-functions -fomit-frame-pointer -DNDEBUG -fno-strict-aliasing --param max-inline-insns-single=1800 -Wall -std=c99 -msse2 -DHAVE_SSE2 -DSFMT_MEXP=1279 -ISFMT-src-1.5.1 -o eulerfast SFMT.c euler2.c

And after all that trouble we got the performance down to 18 seconds. Still slower that `Julia`

‘s 16 seconds.

$ time ./eulerfast Euler : 2.71824 real 0m18.075s user 0m18.085s sys 0m0.001s

Probably, we could do a bit better with more tweaks, and probably exceed `Julia`

‘s performance with some effort. But at that point, I got tired of pushing this further. The thing I love about `Julia`

is how well it is engineered and hassle free. It is quite phenomenal the performance you get out of it, with so little effort. And for basic technical computing things, like random number generation, you don’t have to dig hard for a better library. The “batteries included” choices in the `Julia`

‘s standard library are pretty good.

]]>it is quite nice though that you can rely on the quality of Base for numerics.

— Chris Rackauckas (@ChrisRackauckas) February 8, 2018

Happy e-day of the century! 2 7 18

Curious fact about Euler's Number: Pick a uniformly random number in [0,1] and repeat until the sum of the numbers picked is >1. You'll on average pick e≈2.718… numbers! pic.twitter.com/P0wx9hQu79

— Fermat's Library (@fermatslibrary) February 7, 2018

So I spend a few minutes coding this into `Julia`

function euler(n) m=0 for i=1:n the_sum=0.0 while true m+=1 the_sum+=rand() (the_sum>1.0) && break; end end m/n end

Timing this on my machine, I got

julia> @time euler(1000000000) 15.959913 seconds (5 allocations: 176 bytes) 2.718219862

Gave a little under 16 seconds.

Tried a `c`

implementation

#include <stdio.h> /* printf, NULL */ #include <stdlib.h> /* srand, rand */ #include <time.h> /* time */ double r2() { return (double)rand() / (double)((unsigned)RAND_MAX + 1); } double euler(long int n) { long int m=0; long int i; for(i=0; i<n; i++){ double the_sum=0; while(1) { m++; the_sum+=r2(); if(the_sum>1.0) break; } } return (double)m/(double)n; } int main () { printf ("Euler : %2.5f\n", euler(1000000000)); return 0; }

and compiling with either `gcc`

gcc -Ofast euler.c

or `clang`

`clang -Ofast euler.c`

gave a timing twice as long

$ time ./a.out Euler : 2.71829 real 0m36.213s user 0m36.238s sys 0m0.004s

For the curios, I am using this version of `Julia`

julia> versioninfo() Julia Version 0.6.3-pre.0 Commit 93168a6 (2017-12-18 07:11 UTC) Platform Info: OS: Linux (x86_64-linux-gnu) CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz WORD_SIZE: 64 BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell) LAPACK: libopenblas64_ LIBM: libopenlibm LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

Now one should not put too much emphasis on such micro benchmarks. However, I found this a very curious examples when a high level language like `Julia`

could be twice as fast a `c`

. The `Julia`

language authors must be doing some amazing mojo.

Give the strength to continue

In this land full of pain

How can I continue

to walk among the lifeless…the cowards and the ill-bred

while knowing that the noble ones are dying…or dead

In this reign of fear, you bend…or disappear

Ozymandias now has pathetic pretender of silly sneer and cold command

His colossal wrecks stretch everywhere, further spreading despair

The servile and hungry rise in thunderous applause to Works of their new god

The smell of decay stifling, but the fouls reckon it ambrosia

Cowardice begetting dystopia

A reckoning is overdue….

]]>`Julia`

to find that special number. Have a look at the video here before we continue`using Primes`

function getBigestFactors{T}(a::T,b::T)

c=T()

for k in keys(a)

c[k]=max(a[k],b[k])

`end`

for k in keys(b)

c[k]=max(a[k],b[k])

`end`

c

`end`

superfactors=reduce(getBigestFactors,map(factor,BigInt.(1:100)))

n=prodfactors(superfactors)

`Julia`

can be amazingly expressive. In line 14 are getting all the prime factors of all the numbers from 1 to 100 using the `map`

function. We are using the type `BigInt`

since will be generating pretty large numbers. The standard `Int`

will simply overflow. We then use `reduce`

to get all the common factors for the numbers from 1 to 100.

If your are slight confused lets take a bit slower. First we generate a one dimensional array of `BigInt`

s from 1 to 100.

BigInt.(1:100)

This makes use of `Julia`

‘s very powerful dot factorization syntax. We then `map`

that array into another array of `Primes.Factorization{BigInt}`

through the `factor`

function in the `Primes`

package. The `Primes.Factorization`

type is subtype of `Associative`

. I learned that from reading the implementation.

The cool thing about `Associative`

types that the can be very conveniently accessed. Lets have a closer look.

a=factor(24) println(a)

yields

Primes.Factorization(2=>3,3=>1)

That is

You can then do the following

print(a[2])

yielding

`3`

This means that the number 24 has a 2 raised to the third power as a factor.

Or, be more adventurous and try

print(a[5])

which will print a

`0`

This makes perfect sense as and `1`

is always a factor of any number. We will always get 0 exponent for any number that is not a prime factor of a given number. The `Primes.jl`

package authors have done great use of `Julia`

‘s very powerful type system.

To solve this problem we need to look at all the prime factors for all the numbers from 1 to 100 find the highest exponent of the of each of those prime factors. To do that we implement the `getBigestFactor`

function when does that for any two prime factorizations. We plug that into `reduce`

and *et voilà*!

The superfactors in line 14 are

2^6 ⋅ 3^4 ⋅ 5^2 ⋅ 7^2 ⋅ 11 ⋅ 13 ⋅ 17 ⋅ 19 ⋅ 23 ⋅ 29 ⋅ 31 ⋅ 37 ⋅ 41 ⋅ 43 ⋅ 47 ⋅ 53 ⋅ 59 ⋅ 61 ⋅ 67 ⋅ 71 ⋅ 73 ⋅ 79 ⋅ 83 ⋅ 89 ⋅ 97

Finally we multiply them to get our special number at line 16

`n=69720375229712477164533808935312303556800`

Are we done yet?

As the saying goes “there is more than one way to skin a cat”. Some are more elegant than others.

Here is a more elegant and computationally more efficient way.

Using Primes

function getAllFactors{T}(n::T)

p_factors=filter(isprime,T.(1:n))

exponents=[floor(T,log(p,n)) for p in p_factors]

zip(p_factors,exponents) |> Dict |> Primes.Factorization

`end`

Here we just get the prime numbers in the range from 1 to `n`

(line 3). Using those primes we then get the highest exponents that will not yield a number outside the range (line 4). Finally, we package everything as a `Primes.Factorization`

(line 5).

In `Julia`

REPL we get

julia> getAllFactors(big"100") 2^6 ⋅ 3^4 ⋅ 5^2 ⋅ 7^2 ⋅ 11 ⋅ 13 ⋅ 17 ⋅ 19 ⋅ 23 ⋅ 29 ⋅ 31 ⋅ 37 ⋅ 41 ⋅ 43 ⋅ 47 ⋅ 53 ⋅ 59 ⋅ 61 ⋅ 67 ⋅ 71 ⋅ 73 ⋅ 79 ⋅ 83 ⋅ 89 ⋅ 97 julia> prodfactors(ans) 69720375229712477164533808935312303556800

At that point, you might ask, why stop at a 100. `Julia`

is lightening fast and powerful. Just for fun we can do the same for numbers from 1 to 1000.

julia> getAllFactors(big"1000") 2^9 ⋅ 3^6 ⋅ 5^4 ⋅ 7^3 ⋅ 11^2 ⋅ 13^2 ⋅ 17^2 ⋅ 19^2 ⋅ 23^2 ⋅ 29^2 ⋅ 31^2 ⋅ 37 ⋅ 41 ⋅ 43 ⋅ 47 ⋅ 53 ⋅ 59 ⋅ 61 ⋅ 67 ⋅ 71 ⋅ 73 ⋅ 79 ⋅ 83 ⋅ 89 ⋅ 97 ⋅ 101 ⋅ 103 ⋅ 107 ⋅ 109 ⋅ 113 ⋅ 127 ⋅ 131 ⋅ 137 ⋅ 139 ⋅ 149 ⋅ 151 ⋅ 157 ⋅ 163 ⋅ 167 ⋅ 173 ⋅ 179 ⋅ 181 ⋅ 191 ⋅ 193 ⋅ 197 ⋅ 199 ⋅ 211 ⋅ 223 ⋅ 227 ⋅ 229 ⋅ 233 ⋅ 239 ⋅ 241 ⋅ 251 ⋅ 257 ⋅ 263 ⋅ 269 ⋅ 271 ⋅ 277 ⋅ 281 ⋅ 283 ⋅ 293 ⋅ 307 ⋅ 311 ⋅ 313 ⋅ 317 ⋅ 331 ⋅ 337 ⋅ 347 ⋅ 349 ⋅ 353 ⋅ 359 ⋅ 367 ⋅ 373 ⋅ 379 ⋅ 383 ⋅ 389 ⋅ 397 ⋅ 401 ⋅ 409 ⋅ 419 ⋅ 421 ⋅ 431 ⋅ 433 ⋅ 439 ⋅ 443 ⋅ 449 ⋅ 457 ⋅ 461 ⋅ 463 ⋅ 467 ⋅ 479 ⋅ 487 ⋅ 491 ⋅ 499 ⋅ 503 ⋅ 509 ⋅ 521 ⋅ 523 ⋅ 541 ⋅ 547 ⋅ 557 ⋅ 563 ⋅ 569 ⋅ 571 ⋅ 577 ⋅ 587 ⋅ 593 ⋅ 599 ⋅ 601 ⋅ 607 ⋅ 613 ⋅ 617 ⋅ 619 ⋅ 631 ⋅ 641 ⋅ 643 ⋅ 647 ⋅ 653 ⋅ 659 ⋅ 661 ⋅ 673 ⋅ 677 ⋅ 683 ⋅ 691 ⋅ 701 ⋅ 709 ⋅ 719 ⋅ 727 ⋅ 733 ⋅ 739 ⋅ 743 ⋅ 751 ⋅ 757 ⋅ 761 ⋅ 769 ⋅ 773 ⋅ 787 ⋅ 797 ⋅ 809 ⋅ 811 ⋅ 821 ⋅ 823 ⋅ 827 ⋅ 829 ⋅ 839 ⋅ 853 ⋅ 857 ⋅ 859 ⋅ 863 ⋅ 877 ⋅ 881 ⋅ 883 ⋅ 887 ⋅ 907 ⋅ 911 ⋅ 919 ⋅ 929 ⋅ 937 ⋅ 941 ⋅ 947 ⋅ 953 ⋅ 967 ⋅ 971 ⋅ 977 ⋅ 983 ⋅ 991 ⋅ 997 julia> prodfactors(ans) 7128865274665093053166384155714272920668358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000

Now we get super special number that is 443 digits long. Lets finish by a `Julia`

one liner to get the number of digits

julia> getAllFactors(big"1000") |> prodfactors|> digits |>length 433

There was a cool ratio between the outer and inner circle radii that is expressed as

.

I used `Compose.jl`

to rapidly do that.

using Compose set_default_graphic_size(20cm, 20cm) ϕ=sqrt((pi -2)/(4-pi)) R=10 r=R/ϕ ctx=context(units=UnitBox(-10, -10, 20, 20)) composition = compose(ctx, (ctx, rectangle(-r/√2,-r/√2,r*√2,r*√2),fill("white")), (ctx,circle(0,0,r),fill("blue")), (ctx,circle(0,0,R),fill("white")), (ctx,rectangle(-10,-10,20,20),fill("red"))) composition |> SVG("inscribed.svg")

The rest of the video describes two related approaches for carrying out the counting. Often when I run into complex counting problems, I like to do a sanity check using brute force computation to make sure I have not missed anything. `Julia`

is fantastic choice for doing such computation. It has `C`

like speed, and with an expressiveness that rivals many other high level languages.

Without further ado, here is the Julia code I used to verify my solution the problem.

function unique_combs(n=4)

pat_lookup=Dict{String,Bool}()

for i=0:10^n-1

d=digits(i,10,n) # The digits on an integer in an array with padding

ds=d |> sort |> join # putting the digits in a string after sorting

get(pat_lookup,ds,false) || (pat_lookup[ds]=true)

`end`

println("The number of unique digits is $(length(pat_lookup))")

`end`

In `line 2`

we create a dictionary that we will be using to check if the number fits a previously seen pattern. The loop starting in `line 3`

, examines all possible ordered combinations. The `digits`

function in `line 4`

takes any integer and generate an array of its constituent digits. We generate the unique digit string in `line 5`

using pipes, by first sorting the integer array of digits and then combining them in a string. In `line 6`

we check if the pattern of digits was seen before and make use of quick short short-circuit evaluation to avoid an `if-then`

statement.

`Julia`

calling `C`

. It mimics how one would go about writing `C`

code, wrapping it a library and then calling it from Julia. Today I came across and even more minimal way of doing that while reading an excellent blog on Julia’s syntactic loop fusion. Associated with the blog was notebook that explores the matter further.
Basically, you an write you `C`

in a string and pass it directly to the compiler. It goes something like

C_code= """ double mean(double a, double b) { return (a+b) / 2; } """ const Clib=tempname() open(`gcc -fPIC -O3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f print(f, C_code) end

The `tempname`

function generate a unique temporary file path. On my Linux system `Clib`

will be string like `"/tmp/juliaivzRkT"`

. That path is used to generate a library name `"/tmp/juliaivzRkT.so"`

which will then used in the `ccall`

:

julia> x=ccall((:mean,Clib),Float64,(Float64,Float64),2.0,5.0) 3.5

This approach would be be recommended if are writing anything sophisticated in `C`

. However, it fun to experiment with for short bits of `C`

code that you might like to call from Julia. Saves you the hassle of creating a `Makefile`

, compiling, etc…

I thought it would be fun to try solving it using Julia’s award winning JuMP package. Before we get started, please watch the above video-you might want to pause at 2:24 if you want to solve it yourself.

To attempt this problem in Julia, you will have to install the `JuMP`

package.

julia> Pkg.add("JuMP")

`JuMP`

provides an algebraic modeling language for dealing with mathematical optimization problems. Basically, that allows you to focus on describing your problem in a simple syntax and it would then take care of transforming that description in a form that can be handled by any number of *solvers*. Those solvers can deal with several types of optimization problems, and some solvers are more generic than others. It is important to pick the right solver for the problem that you are attempting.

The problem premises are:

1. There are 50 creatures in total. That includes sharks outside the tanks and fish

2. Each SECTOR has anywhere from 1 to 7 sharks, with no two sectors having the same number of sharks.

3. Each tank has an equal number of fish

4. In total, there are 13 or fewer tanks

5. SECTOR ALPHA has 2 sharks and 4 tanks

6. SECTOR BETA has 4 sharsk and 2 tanks

We want to find the number of tanks in sector GAMMA!

Here we identify the problem as mixed integer non-linear program (MINLP). We know that because the problem involves an integer number of fish tanks, sharks, and number of fish inside each tank. It also non-linear (quadratic to be exact) because it involves multiplying two two of the problem variables to get the total number or creatures. Looking at the table of solvers in the JuMP manual. pick the `Bonmin`

solver from `AmplNLWriter`

package. This is an open source solver, so installation should be hassle free.

julia> Pkg.add("AmplNLWriter")

We are now ready to write some code.

using JuMP, AmplNLWriter # Solve model m = Model(solver=BonminNLSolver()) # Number of fish in each tank @variable(m, n>=1, Int) # Number of sharks in each sector @variable(m, s[i=1:3], Int) # Number of tanks in each sector @variable(m, nt[i=1:3]>=0, Int) @constraints m begin # Constraint 2 sharks[i=1:3], 1 <= s[i] <= 7 numfish[i=1:3], 1 <= nt[i] # Missing uniqueness in restriction # Constraint 4 sum(nt) <= 13 # Constraint 5 s[1] == 2 nt[1] == 4 # Constraint 6 s[2] == 4 nt[2] == 2 end # Constraints 1 & 3 @NLconstraint(m, s[1]+s[2]+s[3]+n*(nt[1]+nt[2]+nt[3]) == 50) # Solve it status = solve(m) sharks_in_each_sector=getvalue(s) fish_in_each_tank=getvalue(n) tanks_in_each_sector=getvalue(nt) @printf("We have %d fishes in each tank.\n", fish_in_each_tank) @printf("We have %d tanks in sector Gamma.\n",tanks_in_each_sector[3]) @printf("We have %d sharks in sector Gamma.\n",sharks_in_each_sector[3])

In that representation we could not capture the restriction that “no two sectors having the same number of sharks”. We end up with the following output:

We have 4 fishes in each tank. We have 4 tanks in sector Gamma. We have 4 sharks in sector Gamma.

Since the problem domain is limited, we can possible fix that by adding a constrain that force the number of sharks in sector Gamma to be greater than 4.

@constraint(m,s[3]>=5)

This will result in an answer that that does not violate any of the stated constraints.

We have 3 fishes in each tank. We have 7 tanks in sector Gamma. We have 5 sharks in sector Gamma.

However, this seems like a bit of kludge. The proper way go about it is represent the number of sharks in the each sector as binary array, with only one value set to 1.

# Number of sharks in each sector @variable(m, s[i=1:3,j=1:7], Bin)

We will have to modify our constraint block accordingly

@constraints m begin # Constraint 2 sharks[i=1:3], sum(s[i,:]) == 1 u_sharks[j=1:7], sum(s[:,j]) <=1 # uniquness # Constraint 4 sum(nt) <= 13 # Constraint 5 s[1,2] == 1 nt[1] == 4 # Constraint 6 s[2,4] == 1 nt[2] == 2 end

We invent a new variable array `st`

to capture the number of sharks in each sector. This simply obtained by multiplying the binary array by the vector

@variable(m,st[i=1:3],Int) @constraint(m, st.==s*collect(1:7))

We rewrite our last constraint as

# Constraints 1 & 3 @NLconstraint(m, st[1]+st[2]+st[3]+n*(nt[1]+nt[2]+nt[3]) == 50)

After the model has been *solved*, we extract our output for the number of sharks.

sharks_in_each_sector=getvalue(st)

…and we get the correct output.

This problem might have been an overkill for using a full blown mixed integer non-linear optimizer. It can be solved by a simple table as shown in the video. However, we might not alway find ourselves in such a fortunate position. We could have also use mixed integer quadratic programming solver such as Gurobi which would be more efficient for that sort of problem. Given the small problem size, efficiency hardly matters here.

]]>