Posts Tagged ‘f#’

Project Euler Problem #231!

November 3rd, 2011

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

Tags: , , ,
Posted in Project Euler | Comments (0)

Find The Divisors Of A Number From Its Prime Factorization.

September 15th, 2011

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

Tags: , , , ,
Posted in Math | Comments (0)

Project Euler Problem #108!

January 12th, 2011

Problem #108 says:

In the following equation x, y, and n are positive integers.

\[ \frac {1}{x} + \frac {1}{y} = \frac{1}{n} \]

What is the least value of n for which the number of distinct solutions exceeds one-thousand?

NOTE: This problem is an easier version of problem 110; it is strongly advised that you solve this one first.

So two things help us with this solution. First, the number of solutions can be calculated as Solutions = ((NumberOfDivisors(n*n) + 1) / 2), Second we can maximize the number of divisors and keep N small by making N a product of small prime numbers. So what this solution does is create a step value out of small primes and then steps through a range looking for solutions. Once it has the list it simply sorts it by the smallest N value and prints it out.

This solution calculates the number of divisors by factoring the number and using one plus the exponents of the factors times one plus the other factors exponents to get the number of divisors. E.g. 2^4 = 16 and 16 has the divisors 1,2,4,8,16 or (4+1). We also get to keep the factorization simple by using our step value this keeps the number of primes we need to factor N down by simply repeating primes that have already been used. However if we search too large a space we will eventually need to expand our list of prime numbers so our factor function doesn’t stop working because it cannot factor part of the number.

This solution works for problem #108 but extending the range it searches shows that it will be too slow for problem #110. I do however have a lead on a solution for #110 that should be much faster than this solution could ever be even though the numbers are far larger.

Solution in F# written in monodevelop. Runs in well under a second.

open System
open System.Collections.Generic
open System.Numerics
open System.Diagnostics

let generatePrimeNumbers max =
    let rec generate number numberSequence =
        if number * number > max then numberSequence else
        let filteredNumbers = numberSequence
                              |> Seq.filter (fun v -> v = number || v % number <> 0I)
                              |> Seq.toArray |> Array.toSeq
        let newNumber = filteredNumbers |> Seq.find (fun v -> v > number)
        generate newNumber filteredNumbers
    generate 2I (seq { for i in 2I..max -> i })

let primes = generatePrimeNumbers 50I |> Seq.toArray

let Factor x =
    let factors = new Dictionary<BigInteger,int64>()
    let factoree = ref x
    for f in primes do
        if (!factoree = 1I) then
            ignore
        else
            while (!factoree % f = 0I) do
                if factors.ContainsKey(f) then
                    factors.[f] <- factors.[f] + 1L
                    factoree := !factoree / f
                else
                   factors.Add (f,1L)
                   factoree := !factoree / f
            ignore
    factors

let DivisorsFromFactoring x =
    let divisors = ref 1L
    let factors = Factor x
    for k in factors.Keys do
        divisors := !divisors * (factors.[k] + 1L)
    !divisors

let SolutionsForGivenN n =
    ((DivisorsFromFactoring (n*n)) + 1L) / 2L

let step = (2L*3L*5L*7L*11L)*2L
let lower = step
let upper = 500000L

let Answer =
    seq[
        for n in [lower .. step .. upper] do
            yield (n , SolutionsForGivenN (BigInteger n))
    ]
    |> Seq.sortBy (fun x -> (fst x))
    |> Seq.filter (fun x -> (snd x) > 1000L)
    |> Seq.take 1

for i in Answer do
    Console.WriteLine i  

Tags: , , ,
Posted in Project Euler | Comments (1)

Point in triangle testing.

January 6th, 2011

So according to my website stats people are ending up on this site after looking for the keywords “check if triangle contain the origin” (yeah I know that’s not good English but it’s what people enter into the search engine), no doubt looking for help with Project Euler Problem #102! (warning that link contains a solution to the problem; don’t click if you want to solve it yourself). However maybe they just want something to check a point and a triangle and see if the triangle encompasses the point. So I went a head and extracted my function for testing whether a triangle encompasses a point from my Project Euler solution and translated it to f# and cleaned it up a bit.

This function uses the half-plane test to check if a triangle contains a point (said point could be the origin if one wants.) this method could be extended for n-polygon shapes but as provided it’s limited to triangles.

Code provided in f#.

#light
open System

type Point(x:float, y) =
  member p.X = x
  member p.Y = y

let Origin = new Point(0.0,0.0)

let Sign (p1:Point) (p2:Point) (p3:Point) =
    (p1.X - p3.X) * (p2.Y - p3.Y) - (p2.X - p3.X) * (p1.Y - p3.Y)

let CheckPoint pt p1 p2 p3 =
    let test1 = (Sign pt p1 p2) < 0.0
    let test2 = (Sign pt p2 p3) < 0.0
    let test3 = (Sign pt p3 p1) < 0.0
    ((test1 = test2) && (test2 = test3))

let A = new Point(-340.0,495.0)
let B = new Point(-153.0,-910.0)
let C = new Point(835.0,-947.0)

let X = new Point(-175.0,41.0)
let Y = new Point(-421.0,-714.0)
let Z = new Point(574.0,-645.0)

Console.WriteLine (CheckPoint Origin A B C)
Console.WriteLine (CheckPoint Origin X Y Z)

Tags: ,
Posted in Math, Software | Comments (0)

Project Euler Problem #78!

January 6th, 2011

Problem #78 says:

Let p(n) represent the number of different ways in which n coins can be separated into piles.
For example, five coins can separated into piles in exactly seven different ways, so p(5)=7.
OOOOO
OOOO O
OOO OO
OOO O O
OO OO O
OO O O O
O O O O O

Find the least value of n for which p(n) is divisible by one million.

This question is actually a slightly more complicated version of #76; and once again we’re dealing with integer partitions. This time however we include the number itself in the list of partitions. Also, the method for generating integer partitions provided in my solution for #76 is too slow for this problem witch deals with many more numbers and much larger numbers.

So the solution is to use Euler’s Pentagonal Theorm and this will let us generate the number of partitions are a good clip. To speed things up we’ll use a Dictionary to store all of the previously computed Pentagonal numbers. We do this mostly because BigInteger operations are slow and duplicating those operations adds up quickly. Another possible speedup could be found in Parallelizing the whole mess and checking multiple numbers in our testing range at the same time. As this is basically an embarrisingly parallel program doing this in parallel should give close to the optimal speedup. If one were to use the f# powerpack the answer function could be rewritten to use PSeq with minimal changes. Another thing that would speedup this program is modifiying it to only tracks the last 7 digits of the number of paritions (7 digits is the minimum number of digits needed to check if a number divides evenly by one million) which would allow the program to use float or decimal values instead of BigInt; which would make the math faster. Lastly one could also add logic to terminate the program when a working number is found provided the program checks the numbers in order; as provided below the program will continue to run untill it has generated the entire list. However this solution runs in about one minute on my machine and that works for me.

Solution provided in f# written in monodevelop.

open System
open System.Diagnostics
open System.Collections.Generic
open System.Numerics

let Pentagonal x = x*(3I*x-1I)/2I

let GeneralPentagonal x =
    if x<0I then
        0I
    else
        if x%2I = 0I then
            Pentagonal ((x/2I)+1I)
        else
            Pentagonal (-((x/2I)+1I))

let mem = new System.Collections.Generic.Dictionary<BigInteger,BigInteger>()
let MemPentagonal x =
    if mem.ContainsKey(x) then
        mem.[x]
    else
        mem.Add(x,GeneralPentagonal x)
        mem.[x]     

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

let answer =
    List [
        let pt = new System.Collections.Generic.List<BigInteger>()
        pt.Add 1I
        for n = 1 to 60000 do
            let r = ref 0I
            let f = ref -1I
            let i = ref 0I
            let cont = ref true
            let k = ref 0I
            while !cont do
                k := MemPentagonal !i
                if !k > (BigInteger n) then
                    cont := false
                else
                    if !i % 2I = 0I then
                        f := 0I - !f
                    r := !r + (!f * pt.[n - (int32 !k)])
                    i := !i + 1I

            let consumeandignore = pt.Add !r
            yield (!r,n)
    ]
    |> Seq.filter (fun n -> (fst n) % 1000000I = 0I)
    |> Seq.head

[<EntryPoint>]
let main args =
    Console.WriteLine answer
    sw.Stop()
    Console.WriteLine sw.Elapsed
    0

Tags: , , ,
Posted in Project Euler | Comments (0)