Archive for May, 2010

Project Euler Problem #132!

May 23rd, 2010

Problem #132 says:

A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k.

For example, R(10) = 1111111111 = 11×41×271×9091, and the sum of these prime factors is 9414.

Find the sum of the first forty prime factors of R(10^(9)).


To solve this problem we take the number and rearrange it a bit, then check for primes that fit the formula runs in about 12 seconds.

Solution in Python. Requires the RabinMiller.py file.

Problem132.py

from RabinMiller import *
from time import *

t0 = time()
def GeoSum (a,  n,  c):
    if (n == 1):
        return 1
    if (n==2):
        return (1+a)%c
    ret = ((GeoSum(a, n/2, c)%c)*((1+pow(a, n/2, c))%c))%c
    if (n&1):
        ret = (ret+pow(a, n-1, c))%c
    return ret

sum = 0
count = 0
i = 3
while (count <40):
    if (MillerRabin(i) and not GeoSum(10, 10**9,  i)):
        count = count + 1
        sum = sum + i
        print count, " divids by ", i
    i = i + 2
print sum
print time() - t0

RabinMiller.py

import sys
import random

def toBinary(n):
  r = []
  while (n > 0):
    r.append(n % 2)
    n = n / 2
  return r

def test(a, n):
  """
  test(a, n) -> bool Tests whether n is complex.

  Returns:
    - True, if n is complex.
    - False, if n is probably prime.
  """
  b = toBinary(n - 1)
  d = 1
  for i in xrange(len(b) - 1, -1, -1):
    x = d
    d = (d * d) % n
    if d == 1 and x != 1 and x != n - 1:
      return True # Complex
    if b[i] == 1:
      d = (d * a) % n
  if d != 1:
    return True # Complex
  return False # Prime

def MillerRabin(n, s = 3):
  """
    MillerRabin(n, s = 1000) -> bool Checks whether n is prime or not

    Returns:
      - True, if n is probably prime.
      - False, if n is complex.
  """
  for j in xrange(1, s + 1):
    a = random.randint(1, n - 1)
    if (test(a, n)):
      return False # n is complex
  return True # n is prime

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

Project Euler Problem #57

May 6th, 2010

Problem #57 says:

It is possible to show that the square root of two can be expressed as an infinite continued fraction.

√ 2 = 1 + 1/(2 + 1/(2 + 1/(2 + … ))) = 1.414213…

By expanding this for the first four iterations, we get:

1 + 1/2 = 3/2 = 1.5
1 + 1/(2 + 1/2) = 7/5 = 1.4
1 + 1/(2 + 1/(2 + 1/2)) = 17/12 = 1.41666…
1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 = 1.41379…

The next three expansions are 99/70, 239/169, and 577/408, but the eighth expansion, 1393/985, is the first example where the number of digits in the numerator exceeds the number of digits in the denominator.

In the first one-thousand expansions, how many fractions contain a numerator with more digits than denominator?


The fastest way to attack this one is to treat it as a sequence and find the pattern required to expand it. for the numerator to take the previous denominator and multiple it by two then add the previous numerator to it, for the denominator you just add the previous numerator and denominator together. This gets really big really fast so we’ll use the biginteger class in .net 4. then we just check the string length of the numerator and denominator and see which is larger. wash rinse and repeat.

Solution in f# and requires .net 4.0 Runs in (5* (10^-6)) seconds or 65 ticks on my netbook.

open System
open System.Diagnostics

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

let num n d = (2I*d)+n
let den n d = n+d

let rec Answer n d c count =
    if c > 1000 then
        count
    else
        if n.ToString().Length > d.ToString().Length then
            Answer (num n d) (den n d) (c+1) (count+1)
        else
            Answer (num n d) (den n d) (c+1) count

sw.Stop()
Console.WriteLine(Answer 3I 2I 0 0)
Console.WriteLine(sw.Elapsed.Ticks)
Console.ReadLine() |> ignore   

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

Project Euler Problem #75!

May 5th, 2010

Problem #75 says:

It turns out that 12 cm is the smallest length of wire that can be bent to form an integer
sided right angle triangle in exactly one way, but there are many more examples.

12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)

In contrast, some lengths of wire, like 20 cm, cannot be bent to form an integer sided right
angle triangle, and other lengths allow more than one solution to be found; for example,
using 120 cm it is possible to form exactly three different integer sided right angle triangles.

120 cm: (30,40,50), (20,48,52), (24,45,51)

Given that L is the length of the wire, for how many values of L ≤ 1,500,000 can exactly
one integer sided right angle triangle be formed?


Another counting problem but here we have to avoid duplicate entries that can come from a variety
of places. This solution uses the 2 number approach for generating Pythagorean triples, which it then finds
the lengths of the triple and catalogs it in an array to track duplicates, Then counts all the values of the
array that equal 1 for the final answer.

Solution is provided in c#/mono. Runs in under one second.

using System;
using System.Collections;
using System.Diagnostics;

namespace Problem75
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			Stopwatch sw = new Stopwatch();
			sw.Start();
			int limit = 1500000;
			int adjlimit = (int)(Math.Sqrt((float)limit));
			ArrayList tri = new ArrayList(limit);

			for (int st = 0; st < limit+1; st++)
				tri.Add(0);

			for (int i = 1; i <= adjlimit; i=i+2)
			{
				for (int j = 2; j<= adjlimit; j=j+2)
				{
					if (GCD(i,j) == 1)
					{
						int sum = Math.Abs((j*j) - (i*i)) + (2*j*i) + ((i*i) + (j*j));
						for (int s = sum; s < limit; s=s+sum)
						{

							tri[s] = (int)tri[s]+1;
						}
					}
				}
			}
			Console.WriteLine(TallyResults(tri));
			sw.Stop();
			Console.WriteLine(sw.Elapsed);
			Console.ReadLine();
		}
		public static int TallyResults(ArrayList a)
		{
			int RetVal = 0;
			foreach (int i in a)
			{
				if (i == 1)
					RetVal= RetVal+1;
			}
			return RetVal;
		}
		public static int GCD (int a, int b)
		{
			if (b == 0)
				return a;
			else
				return GCD(b,(a%b));
		}
	}
}

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

Project Euler Problem #99!

May 1st, 2010

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

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

Random Snippets.

May 1st, 2010


Here are some random snippets of f# code I had lying around. mostly simple stuff for basic stats calculations. Not the most complex stuff but maybe its useful to somebody trying to learn f#.

// Learn more about F# at http://fsharp.net
open System
open System.Net

//basic arithmetic mean
let Mean (x:List<float>) =
    let sum = x |> Seq.sum
    let items = x |>Seq.length
    sum / (float items)

//basic variance method
let Variance x =
    let m = Mean x
    let deviation =
        x
        |> Seq.map (fun x -> Math.Pow((x-m), 2.))
        |> Seq.toList
    Mean deviation

//Calculates standard deviation using above variance method
let StandardDeviation x =
    Math.Sqrt(Variance x)

//recursively calculates the factorial of x
//will overflow long before running into recursion limits
//at larger values it will wrap back around to 0
let rec Factorial x =
    if x = 0I then
        1I
    else
        x * Factorial(x-1I)

//calculates combinations
let Combinations n r =
    (Factorial n)/ ((Factorial r)*(Factorial (n-r)))

//calculates permutations
let Permutations n r =
    (Factorial n)/ (Factorial (n-r))

//uses Euclids formula for the greatest common divisor
let rec GCD a b =
    if b = 0UL then
        a
    else
        GCD b (a%b)

//returns an array representing the binary version of x
let toBinary (x:uint64) =
    let rec con (x:uint64) r =
        if x = 0UL then
            r
        else
            con (x/2UL) ([x%2UL] @ r)
    con x []

// gets ticker price using yahoo
let GetTickerPrice (ticker:string) =
    let url =
        "http://download.finance.yahoo.com/d/quotes.csv?s="+
        ticker+
        "&f=sl1d1t1c1ohgv&e=.csv"
    let wc = new WebClient()
    let info = wc.DownloadString(url).Split([|','|])
    float info.[1] // info.[1] is the current price will return 0.0 on invalid ticker

Tags: ,
Posted in Software | Comments (0)