So I decided to dig deeper. Basically the standard crand() 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.

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.

On e-day, I came across this cool tweet from Fermat’s library

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

$ 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 me a reason to hang on
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

Mathematical recreation is fun. I came across a nice Mind your Decisions video and thought it would be cool to use Julia to find that special number. Have a look at the video here before we continue

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 BigInts 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à!

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.

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).

Recently, I watched a cool mind your decsions video on an inscribed circle and rectangle puzzle. In the video they showed a diagram that was not scale. I wanted to get a sense of how these differently shaped areas will match.

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

I came across a neat math puzzle involving counting the number of unique combinations in a hypothetical lock where digit order does not count. Before you continue, please watch at least the first minute of following video:

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.

Earlier I presented a minimal example of 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:

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…

Recently I came across a nice Ted-Ed video presenting a Fish Riddle.

I thought it would be fun to try solving it using Julia’s award winningJuMP 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] == 2end# 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

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.

We humans have always seemed to have sort of ineptitude of not being able to notice the great diversity of life on our planet. However, any feeling of emptiness is merely an illusion to a blind mind. Through the course of human evolution, our minds and our vast intellect were the most essential tools for the survival of our species.

But if you were to condense Earth’s entire existence into 24hrs, then we would have only been on this planet for a mere… 3 seconds !

Our minds have enabled us to think deeply about our purpose in life and to do incredible things. But we have become so blinded by the power we have that we have not looked at what insidious and irrevocable damage we have caused to the balance nature.

Do you know what is the most dangerous flaw our human minds have kept for the last 200,000 of our existence?…

Not knowing where to look or whom to blame for our mistakes. We have just become so ignorant, so blinded, and so utterly self-centered about only our well-being that we have not noticed issues like increasing famine cases around Nigeria and Somalia, or the rapid occurrences of some of the largest storms ever seen. So pernicious we have become to ourselves. How have we not noticed this catastrophic conundrum at all ?

Sure we have accomplished some of the most incredible feats of innovation in the natural world but we have not at all been wise enough to contemplate on the consequences our triumph cost. It’s good to be smart but not too smart for you own good.

We have been under the highly foolish illusion that what our minds and ideas have enabled us to do would not exponentially changing many highly essential factors for the greater good of our species and all life on Earth.

We have believed for so long that supposedly “small” problems that don’t fall into our best interests at heart should not be look at seriously for even the smallest or seemingly “simple” could bring the down fall of us all.

Do you want what is the deadliest disease ever?

Tuberculosis?

Measles?

Ebola ?

AIDs? *It’s a single, influential idea which makes us blind to what heinous and atrocious things happen by every passing day.

Is it highly ethical of us as a species to fall into the hands of such a vicious free-thinking paralytic? Would it help to just think that with a flick of a finger all of the problems that this crisis that we have fallen into would be solved ?

An idea can be very deadly by it self but there is a much more vicious motivator for simply ignoring the worsening problems we have caused.

We are silenced by payment of money. A bribery which many us might lightly take. Our minds can be easily persuaded into doing any thing by the motivation of payment or wealth. If we are given everything we could possibly need it’s still not enough to content us.

Our greedy and cunning simply want more of this pleasurable wealth. Some people may see a hundred pound not in there hands from there work. I think we have just become slaves to a printed peace of paper.