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

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

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```