01/11/18

# Factorizations to get a very special number

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

1. using Primes
2. 
3. function getBigestFactors{T}(a::T,b::T)
4.     c=T()
5.     for k in keys(a)
6.         c[k]=max(a[k],b[k])
7.     end
8.     for k in keys(b)
9.         c[k]=max(a[k],b[k])
10.     end
11.     c
12. end
13. 
14. superfactors=reduce(getBigestFactors,map(factor,BigInt.(1:100)))
15. 
16. 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 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 $$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.

1. Using Primes
2. function getAllFactors{T}(n::T)
3.     p_factors=filter(isprime,T.(1:n))
4.     exponents=[floor(T,log(p,n)) for p in p_factors]
5.     zip(p_factors,exponents) |> Dict |> Primes.Factorization
6. 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