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

Problem #124 wants us to:

The radical of n, rad(n), is the product of distinct prime factors of n. For example, \(504 = 2^3 × 3^2 × 7\), so rad(504) = 2 × 3 × 7 = 42.

Let E(k) be the kth element in the sorted n column; for example, E(4) = 8 and E(6) = 9.

If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).

This problem is exceedingly easy once we calculate the radicals for the numbers 1-100000. Just create a simple class to remember the n, and rad values and populate those fields and the problem is 9/10ths of the way solved. From here we can use the power of Linq, which I’m liking more and more as time goes on, and do an OrderBy on the rad values. Now this seems like it would be insufficient to solve the problem because if the radical values are the same we need to sort by the n value. Well Linq has us covered; simply add a ThenBy after the OrderBy and you can sort by a second value. This means that the overwhelming majority of the code in the solution is prime number generation, finding distinct factors, and calculation of the radicals.

A second solution without linq would be to sort by the radical value and get the radical value of the 10000th item. Then find the index of the first number with that radical and store it. Then get a list of only the numbers with that radical and sort this new list by the n value. Using the previously stored index you can then calculate what item in the new list should be the answer. But this is alot more complicated than a simple Linq expression.

A third method would be to write a custom IComparer for the class that holds the radical + n values. But this is a lot more work than is needed and would make the solution much less readable in my opinion. Besides an up to date .Net implementation (Microsoft.Net or Mono) should have Linq available.

Solution below written in c#/mono and requires Linq and .Net 4. Runs in ~5 seconds.

using System;
using System.Collections.Generic;
using System.Linq;

namespace Problem124
{
	class MainClass
	{
		class N
		{
			public ulong n = 0;
			public ulong rad = 0;

			public N (ulong num)
			{
				this.n = num;
				this.rad = rad(num);
			}
		}
		public static ulong max = 100000;
		public static List<ulong> primes = GeneratePrimes(max);

		public static void Main (string[] args)
		{
			List<N> NCases = new List<N>();
			for (ulong n = 1; n<=max;n++)
			{
				NCases.Add(new N(n));
			}
			N answer = NCases.OrderBy(x => x.rad)
				.ThenBy(x => x.n).ToList()[10000-1];

			Console.WriteLine(answer.n);
		}
		public static ulong rad(ulong n)
		{
			List<ulong> factors = DistinctFactors(n);
			ulong rad = 1;
			foreach(ulong f in factors)
				rad *= f;
			return rad;
		}
		public static List<ulong> DistinctFactors (ulong n)
		{
			List<ulong> factors = new List<ulong>();

			if(primes.Contains(n))
			{
				factors.Add(n);
				return factors;
			}
			ulong sqrt = (ulong)Math.Sqrt((double)n);
			foreach (ulong p in primes)
			{
				if (p > sqrt)
					break;
				if (n%p == 0)
				{
					factors.Add(p);
					while (n%p == 0)
						n = n/p;
				}
			}
			if (n!= 1ul)
				factors.Add(n);
			return factors;
		}
		static List<ulong> GeneratePrimes(ulong num)
		{
			List<ulong> RetValue = new List<ulong>();
			RetValue.Add((ulong)2);
			RetValue.Add((ulong)3);
			ulong Stepper = 5;
			ulong Check = 1;
			while (Stepper <= num)
			{
				foreach (ulong i in RetValue)
				{
					if (Stepper % i == 0)
					{
						Check = 0;
						break;
					}
					if (Math.Sqrt(Stepper) <= i)
					{
						break;
					}
				}
				if (Check == 1)
				{
					RetValue.Add(Stepper);
				}
				Check = 1;
				Stepper++;
				Stepper++;
			}
			return RetValue;
		}
	}
}

Problem #23 says:

Let p(n) be the nth prime: 2, 3, 5, 7, 11, …, and let r be the remainder when \((p(n)-1)^{n}+(p(n)+1)^{n}\) is divided by p(n)2.

For example, when n = 3, p(3) = 5, and 43 + 63 = 280 ≡ 5 mod 25.

The least value of n for which the remainder first exceeds \(10^{9}\) is 7037.

Find the least value of n for which the remainder first exceeds \(10^{10}\).

This can be solved using a modified version of the solution for Problem #120. We don’t even need to keep track of all of the remainders because the problem only concerns itself with the first remainder to exceed \(10^{10}\). So using our maximum remainder from the previous Problem #120 which was \(2*n*a\) we can modify it to be \(2*p_{n}*n\). From this point we just need a collection of prime numbers large enough to get the Job done; I chose to use a list of primes going up to 300000 because this creates a list of 25998 prime numbers and \(25998 * 2 * p[25998] > 1.5 x 10^{10}\) which gives me all of the data I needed to find a solution exceeding \(10^{10}\).

Solution below in c#/mono. runs in about .113 seconds.

using System;
using System.Collections.Generic;

namespace Problem123
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			ulong upperbound = 10000000000ul;
			List<ulong> p = GeneratePrimes(300000);
			p.Insert(0,0UL);
			Console.WriteLine("primes done");
			for (uint n = 1; n < 300000ul; n++)
			{
				if (n%2 ==0ul)
					continue;
				ulong r = 2 * p[(int)n] * n;
				if (r > upperbound)
				{
					Console.WriteLine(n);
					Console.WriteLine(p[(int)n]);
					Console.WriteLine(r);
					break;
				}
			}

		}
		static List<ulong> GeneratePrimes(ulong num)
		{
			List<ulong> RetValue = new List<ulong>();
			RetValue.Add((ulong)2);
			RetValue.Add((ulong)3);
			ulong Stepper = 5;
			ulong Check = 1;
			while (Stepper <= num)
			{
				foreach (ulong i in RetValue)
				{
					if (Stepper % i == 0)
					{
						Check = 0;
						break;
					}
					if (Math.Sqrt(Stepper) <= i)
					{
						break;
					}
				}
				if (Check == 1)
				{
					//Console.WriteLine(((float)Stepper / (float)num) * 100);
					RetValue.Add(Stepper);
				}
				Check = 1;
				Stepper++;
				Stepper++;
			}
			return RetValue;
		}
	}
}

Also just as a note I’ve moved from jsMath to MathJax to see if it solves any of the issues I’ve been having. If you see something that might be an equation but it’s all screwy please point it out in the comments.

Problem #80 is stated as:

It is well known that if the square root of a natural number is not an integer, then it is irrational. The decimal expansion of such square roots is infinite without any repeating pattern at all.

The square root of two is 1.41421356237309504880…, and the digital sum of the first one hundred decimal digits is 475.

For the first one hundred natural numbers, find the total of the digital sums of the first one hundred decimal digits for all the irrational square roots.

This problem actually lead me to a new and interesting way of finding square roots described in a paper titled “Square roots by subtraction” by Frazer Jarvis. Frazer claims that it is an old math algorithm from Japan but it’s not the algorithms history that makes it useful. Instead it’s the fact that the algorithm allows us to get arbitrarily precise square roots using basic subtraction and multiplication (I felt the title was a little odd but oh well).

The algorithm works like this take your number \(n\) and make a new number \(a=5n\) and let \(b=5\). The you apply one of two rules depending on whether \(a\) is larger than \(b\). If \(a\) is larger then you apply Rule 1: which is to set \(a=a-b\) and set \(b=b+10\). If \(b\) is larger you apply Rule 2: set \(a= a*100\) and set \(b=((b-5)*10)+5)\). Just repeat those rules using the \((a,b)\) from the previous round until your answer (which is \(b\)) has the desired precision.

So for the square root of 2 we get (10,5) \(Rule 1\atop \longrightarrow\) (5,15) \(Rule 2\atop \longrightarrow\) (500,105)\(Rule 1\atop \longrightarrow\) (395,115) \(Rule 1\atop \longrightarrow\) (280,125) \(Rule 1\atop \longrightarrow\) (155,135) \(Rule 1\atop \longrightarrow\) (20,145) \(Rule 2\atop \longrightarrow\) (2000,1405) \(Rule 1\atop \longrightarrow\) (595,1415) and as we continue b gets closer and closer to the square root of 2. With this algorithm we are only limited by the size of the integer class we are using. As an added bonus this algorithm uses an integer class so we don’t have to worry about rounding errors. This means that with c# and .Net 4 we can using the BigInteger class from System.Numerics and generate the digits of the square root of 2 forever*.

So once we have this method all figured out what’s left to do? Well as it turns out not a whole lot; we just need a digital sum and a list of perfect squares to exclude and we’re golden. Solution provided in c#/mono and runs in 1 second or so. Could easily be made faster (at the very least one could add some parallelism).

using System;
using System.Numerics;
using System.Collections.Generic;

namespace Problem80
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			List<int> squares = new List<int>();
			for (int i = 1; i <= 10; i++)
			{
				squares.Add(i*i);
			}
			uint sum = 0;
			for (int i = 1; i <=100; i++)
			{
				if (!squares.Contains(i))
					sum += DigitalSum(i);
			}
			Console.WriteLine(sum);

		}
		public static uint DigitalSum (int i)
		{
			BigInteger digits = BigInteger.Parse(SquareRoot(i).ToString().Substring(0,100));
			uint digitsum = 0;
			while (digits != BigInteger.Zero)
			{
				uint temp = (uint)(digits%10);
				digitsum += temp;
				digits = (digits -temp) / 10;
			}

			return digitsum;
		}
		public static BigInteger SquareRoot (int i)
		{
			BigInteger a = 5 * i;
			BigInteger b = 5;
			while (b.ToString().Length < 102)
			{
				if (a >= b)
				{
					a = a-b;
					b += 10;
				}
				else
				{
					a *= 100;
					b = ((b-5)*10)+5;
				}
			}
			return b;
		}
	}
}

* Memory limits and the physics of our universe prevent a person or computer from actually doing this.