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

`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 \(2^3 3\)

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 \(5^0=1\) 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