Problem #231 says:

The binomial coefficient \( ^{10}C_{3} = 120\).
\(120 = 2^{3} \ 3 \ 5 = 2\ 2\ 2\ 3\ 5\), and \(2 + 2 + 2 + 3 + 5 = 14\).
So the sum of the terms in the prime factorisation of \(^{10}C_{3}\) is 14.

Find the sum of the terms in the prime factorisation of \(^{20000000}C_{15000000}\).

There’s a simple formula to calculate this. I don’t remember where I found it I wrote this solution a while ago and never posted it. Uses a bunch of prime numbers and the sum of factors to get the answer without brute force.The other way I suppose one could do this is there is a definition for the binomial coefficient as a product of a bunch of stuff. If we could eliminate all rounding error that approach could lead to a working solution.

Solution provided in F# and uses a Sieve of Atkins that seems to be slightly faster than other prime number generators I’ve been using.

#light

open System
open System.Collections
open System.Collections.Generic
open System.Diagnostics
open System.Threading
open Microsoft.FSharp.Collections
open Microsoft.FSharp.Math

let Atkins limit =

    let is_prime = new Dictionary<uint64,bool>()
    let keys = [5UL..limit]
    for i in keys do
        is_prime.Add(i,false)

    let sqrt_limit = uint64 (Math.Sqrt (float limit))
    let xys = [1UL..sqrt_limit]
    for x in xys do
        for y in xys do
            let mutable n = 4UL*(x*x) + (y*y)
            if (n <= limit && (n % 12UL = 1UL || n % 12UL = 5UL)) then
                is_prime.[n] <- not is_prime.[n]
            n <- 3UL*(x*x) + (y*y)
            if (n <= limit && (n % 12UL = 7UL)) then
                is_prime.[n] <- not is_prime.[n]
            n <- 3UL*(x*x) - (y*y)
            if ((x>y) && n <= limit && (n % 12UL = 11UL)) then
                is_prime.[n] <- not is_prime.[n]

    for i in keys do
        if is_prime.[i] then
            let max = limit / (i*i)
            for k in [1UL..max] do
                is_prime.[k*(i*i)] <- false
    [
    yield 2UL;yield 3UL;
    for i in keys do
        if is_prime.[i] then
            yield i
    ]

let primes = Atkins 20000000UL |> Seq.toArray

let SumOfFactors n p  =
    let nt = ref n
    [
        while (!nt > 0UL) do
            nt := !nt / p
            yield !nt
    ] |> Seq.sum

let Problem231 n k =
    primes
    |> PSeq.map (fun i -> ((SumOfFactors n i) - (SumOfFactors k i) - (SumOfFactors (n-k) i) )* i)
    |> Seq.sum

Console.WriteLine "Primes Finished"
let sw = new Stopwatch()
sw.Start()
Console.WriteLine  (Problem231 20000000UL 15000000UL)
sw.Stop()
Console.WriteLine(sw.Elapsed)
Console.ReadLine() |> ignore

The following code takes a number k and its distinct factors f and creates a SortedSet containing all of the divisors of that number. It does this by expanding the set of divisors but multiplying each member of the already known divisors by each prime factor and checking to make sure that the input number k mod the resulting number is 0. Doing this a few times gives a complete list of all of the divisors. There are other methods including recursion that could be used but for what I wrote this to do it should be efficient enough. This program was written to handle numbers of the form \( 10^{n} \) which have only two distinct factors (2,5). While it might be able to handle other numbers with a small set of distinct factors I believe that it would quickly become unwieldy when given a large number with a large prime factor list.

Also updates will be sporadic for a while because it’s my first semester at my new university and I’ve got a heavy course load.

#light

open System
open System.Collections
open System.Collections.Generic

let ExpandSet s f max=
    let s2 = new SortedSet<int64>()
    for i in s do
        for j in f do
            let n = i*j
            if (n <= max && max % n = 0L) then
                s2.Add(n) |> ignore
    s2.UnionWith(s)
    s2.UnionWith(f)
    s2

let fac =
    let t = new SortedSet<int64>()
    t.Add(2L)
    t.Add(5L)
    t

let DivisorsOf k f =
    let mutable divs = new SortedSet<int64>()
    divs.Add(1L) |> ignore
    let tries = int64(Math.Ceiling(Math.Log(float k)/Math.Log(2.)))
    for i in [1L .. tries] do
        Console.WriteLine(i) |> ignore
        divs <- ExpandSet divs f k
    divs

let PrintSet (s:SortedSet<int64>) =
    for i in s do
        Console.WriteLine(i)

let t = DivisorsOf 1000000000L fac

PrintSet t

So I found myself working on a small personal project in F# to gain more familiarity in the language when I decided the best way to solve my data storage needs would be a small local database file. So I went looking for ways to get sqlite into F# projects and couldn’t find a decent example. So I’ve started going through an example project in c# and trying to make it work in F#. I’ve got database creation and item insertion working so I figured what the heck I’ll post it here to make up for all of the posts I haven’t been writing lately.

I’m using the sqlite-net class from here compiled as a dll in a c# project. I’m then using that dll in my F# project. So if you’re looking for a pure f# implementation you at the wrong place. However its all il code at that level anyway. So armed with this dll I then create a f# type that looks a little something like this:

type Site (url:string) =
    let mutable id = new int()
    let mutable BD = ""
    let mutable site = url
    let mutable shown = false
    let mutable exemplarcontributor = false
    [<PrimaryKey>] [<AutoIncrement>]
    member x.ID with get() = id
                and set v = id <- v
    member x.ExemplarContributor with get() = exemplarcontributor
                                    and set v = exemplarcontributor <- v
    member x.Shown with get() = shown
                        and set v = shown <- v
    member x.BreakDown with get() = BD
                        and set v = BD <- v
    [<Indexed>]
    member x.Site with get() = site
                   and set v = site <- v
    member x.Url = url
    new() = Site "www.google.com"

now I should mention that this is directly from the test applet I’ve been using to get sqlite to work. as such it can probably be cleaned up quite a bit. Basically we have a URL and a few other details of a site and a few attributes. The ID is a primary key and auto incremented as shown above, note that this NEEDS to be put on the member declaration not the actual id. Also note that we need getter and setters on anything that we want to show up in the database, as such we need to make everything mutable. I’m still looking to see if there is another way to do that but this works for now. Anyways, armed with this type we can go to the next section, creating a database type that we’ll be adding custom stuff to in addition to inheriting from SQLiteConnection. And that looks like:

type Database (path) =
    inherit SQLiteConnection(path)
    member this.Setup = base.CreateTable<Site>()

Pretty basic right? All we’ve done is give the Database type a setup method that can be used to add the Site table to the database. If the database already has a Site table it’s not replaced but a safety check could easily be added to this method to make sure we don’t overwrite anything; however this is a throw away app and I’m not going to add it. Next we start doing stuff with these types. In fact we do this:

let D = new Database("test.sqlitedb")
D.Setup |> ignore

let s = new Site "www.google.com"

D.Insert(s) |> ignore

D.Commit|>ignore

Here we create a Database that we call D, call its setup method and pipe that to ignore (it returns an int to indicate errors, but I don’t want to catch it). Then we create a new instance of the Site type, do an insert again piping it to ignore and finally we commit everything and pipe that to ignore. And that’s it. We have a simple program that builds a useless database and adds a single entry to it. If you run this program multiple times it’ll just keep adding instances of Site to the database because it doesn’t overwrite the file.

Problem #99 says:

Comparing two numbers written in index form like 2^(11) and 3^(7) is not difficult,
as any calculator would confirm that 2^(11) = 2048 < 3^(7) = 2187.

However, confirming that 632382^(518061) > 519432^(525806) would be much more difficult,
as both numbers contain over three million digits.

Using base_exp.txt (right click and ‘Save Link/Target As…’), a 22K text file containing one
thousand lines with a base/exponent pair on each line, determine which line number has the greatest numerical value.

This is a problem that seems kinda complicated at first. Each of the base exponent sets creates a number that easily dwarfs anything a computer can currently handle. The solution is to do something else that still lets us compare them to one another. So instead of expanding the exponent base pair what we do is multiply the exponent by the natural log of the base and use that value. Using this method reading the file in is almost harder than solving the task.

Solution in F# requires .net 4, runs in under one second.

open System
open System.IO
open System.Diagnostics

let sw = new Stopwatch()
sw.Start()

let test b e =
    e * Math.Log(b)    

let CheckSets =
    let f = new StreamReader("base_exp.txt")
    let mutable lines = []
    let mutable line = f.ReadLine()
    while line <> null do
        lines <- lines @ [line]
        line <- f.ReadLine()
    lines

let mutable mv = 0.0
let mutable ml = 0
let mutable ln = 0

let FindAnswer =
    for s in CheckSets do
        ln <- ln + 1
        let ss = s.Split(",".ToCharArray())
        let b = float ss.[0]
        let e = float ss.[1]
        if (test b e) > mv then
            mv <- (test b e)
            ml <- ln
FindAnswer
Console.WriteLine(ml)
sw.Stop()
Console.WriteLine(sw.Elapsed)
Console.ReadLine() |> ignore

Problem #58 says:

Starting with 1 and spiraling anticlockwise in the following way, a square spiral with side length 7 is formed.

It is interesting to note that the odd squares lie along the bottom right diagonal,
but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime;
that is, a ratio of 8/13 ≈ 62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed.
If this process is continued, what is the side length of the square spiral for which the ratio of
primes along both diagonals first falls below 10%?

The below solution simply counts the offset between corners and checks that number to see if its prime or composite then the next cycle after its check it updates the offset amount and continues this way until its collected a list of prime and composite numbers that have less than 10% prime numbers. very simple, no optimizations applied. The solution could be made faster by using a simpler test and also by not keeping a list of composite numbers but a counter instead.

Solution in f# and requires the .net 4.0 run-time. runs in ~50 seconds.

#light
module ProjectEuler

open System
open System.Diagnostics
open Microsoft.FSharp.Linq
open Microsoft.FSharp.Collections
open Microsoft.FSharp.Math
open System.Numerics

let toBinary (n:BigInteger) =
    let mutable m = n
    let mutable r = []
    while m > BigInteger.Zero do
        r <- r @ [(m % (BigInteger 2))]
        m <- m / BigInteger 2
    r

let test (a:BigInteger) (n:BigInteger) =
    let (b:List<BigInteger>) = toBinary (n - BigInteger.One)
    let mutable d = BigInteger.One
    let mutable Prime = false
    let CheckList = [0 .. b.Length-1 ] |> List.rev
    for i in CheckList do
        let x = d
        d <- (d*d) % n
        if (d = BigInteger.One && x <> BigInteger.One && x <> n-BigInteger.One) then
            Prime <- true // complex
        if b.[i] = BigInteger.One then
            d <- (d*a) % n
    if d <> BigInteger.One then
        Prime <- true //complex
    Prime //if its still false then prime

let MillerRabin (n:uint64) (s:int) =
    let Rand = new System.Random()
    let mutable Prime = true
    if n < Convert.ToUInt64(Int32.MaxValue) then
        for j in [1 .. s+1] do
            let a = Rand.Next(1, Convert.ToInt32(n)-1)
            if (test (BigInteger a) (BigInteger n)) then
                Prime <- false
    else
        for j in [1 .. s+1] do
            let a = Rand.Next(1, Int32.MaxValue)
            if (test (BigInteger a) (BigInteger n)) then
                Prime <- false
    Prime

let watch = new Stopwatch()
watch.Start()

let rec Answer (currentnum:uint64) (step:uint64) (prime:List<uint64>) (composite:List<uint64>) (tick:bool) =
    let check = (float prime.Length)/((float composite.Length)+(float prime.Length))
    if check < 0.1 && currentnum > 100UL  then
        step-1UL
    else
        if (tick) then
            if (MillerRabin currentnum 10) then
                Answer (currentnum+step) step (currentnum::prime) composite (false)
            else
                Answer (currentnum+step) step (prime) (currentnum::composite) (false)
        else
            if (MillerRabin currentnum 10) then
                Answer (currentnum+step) (step+1UL) (currentnum::prime) composite (true)
            else
                Answer (currentnum+step) (step+1UL) (prime) (currentnum::composite) (true)

Console.WriteLine(Answer 3UL 2UL [] [1UL] true )
watch.Stop()
Console.WriteLine(watch.Elapsed)
Console.Beep() |> ignore
Console.ReadKey() |> ignore