Problem #159 asks:

The function mdrs(n) gives the maximum Digital Root Sum of \(n\). So \(\text{mdrs}(24)=11\).
Find \(\text{mdrs}(n)\) for \(1\lt n\lt 1,000,000\).

First we need to know that the Digital root of a number is equal to

$$
f(x) = \left\{
\begin{array}{lr}
n \pmod{9} & : n \neq 0 \pmod{9}\\
9 & : n = 0 \pmod{9}
\end{array}
\right.
$$

We then need the simple relations that \(\text{mdrs}(n) = n\) for prime \(n\), otherwise we look at the list of divisors such that \( a*b = n\) and apply \(\text{mdrs}(n) = Max(\text{mdrs}(a) + \text{mdrs}(b), \text{drs}[n])\) for each pair a&b. setting this up can be done in two ways: 1) Factor the number and generate each MDRS individually which for the number range we’re looking at isn’t too bad or 2) a dynamic programming approach which avoids factoring and creates all the MDRS’s in one go so we just have to add them up. The solution below uses the second approach, runs in about 3 tenths of a second for the given Maximum.

using System;
using System.Diagnostics;
using System.Linq;


namespace Problem159
{
    class Program
    {
        private const long Max = 1000000;

        static void Main()
        {
            var drs = new long[Max];
            var sw = new Stopwatch();
            sw.Start();
            for (long i = 2;i < Max; i++)
            {
                if (drs[i] == 0)
                    drs[i] = i%9 ==0? 9 : i%9;
                else if (drs[i]<10)
                {
                    long c = i % 9 == 0 ? 9 : i % 9;
                    drs[i] = Math.Max(c, drs[i]);
                }
                long cur = i;
                for (long fac = 2; fac <=i; fac++)
                {
                    cur += i;
                    if (cur >= Max)
                        break;
                    drs[cur] = Math.Max(drs[i] + drs[fac], drs[cur]);

                }
            }
            sw.Stop();
            var sum = drs.Sum();
            Console.WriteLine(sum);
            Console.WriteLine(sw.Elapsed);
            Console.ReadLine();
        }
    }
}

Modular arithmetic also know as congruences… What is it? Well, it shows up a lot in programming. Even in the code on this site I’ve used things like “n % 2 == 0″ which reads “n is 0 modulo 2″. But what does this mean and can it be used for more?

The common example given of modular arithmetic is that of a clock. We count, on a normal civillian clock, 1-12 then instead of going to a thirteenth hour we start back at 1. This however has an obvious problem in that there is no “zeroth” hour when looking at a clock but in modular arithmetic zeros are allowed. Instead think of astronomical time/military time. In this system we go from 00:00 – 23:59 and instead of going to 24:00 start back at 00:00.

So what does this mean in terms of mathematics? Well in math the definition of modular arithmetic is \( n \equiv b \pmod{m} \Leftrightarrow m | (n-b)\). So we consider two numbers modularly equivelent if the modulo (in this case m) divides the difference of the two numbers. If the two numbers are equivilent modulo m then there difference is an integer multiple of m.

Now consider \( a \equiv b \pmod{n}\) and \(c \equiv d \pmod{n}\). From this information and the properties of modular arithmetic we get \( b \equiv a \pmod{n}\), \( a+c \equiv b+d \pmod{n}\) and \( a*c \equiv b*d \pmod{n}\). So we can add, subtract and multiply without any trouble. Now what about division? It turns out that division is rather complicated under modular arithmetic. Consider everyday division, we know that we cannot divide by zero because \( \frac{n}{0}\) is undefined for all \(n\). In modular arithmetic we add to this, consider \(\frac{n}{km}\). Modulo m the demoniator is zero for all \(k\), which is undefined. So now the question becomes can we say when division is allowed?

The answer is yes! Division is allowed when the modular inverse of a number exists modulo \(m\). In other words there exists a number \(a\) such that \(a*n \equiv 1 \pmod{m}\). This inverse exists when \(\text{gcd}(n,m) = 1\). To find this inverse one needs the extended Euclidean algorithm (to be described later). But if this inverse exists then we can effectivly “divide” by multiplying by the modular inverse. One last thing to consider are statements like \( n \equiv -a \pmod{m} \) where \(a> 0\) . This can be rewritten as \( n \equiv m-a \pmod{m} \) and we know that \(m-a \geq 0\) because modulo m we know \(a \in \{1,2,…, m-1\}\).

Now the title mentions Wilson’s theorem. This is a theorem using congruences to make a statement about prime numbers. The theorem states “p is a prime number if and only if \((p-1)! \equiv -1 \pmod{p}\)”. For those just starting in mathematics the statement “if and only if” implies that if one side is true then the other side is also true. In this case if \((p-1)! \equiv -1 \pmod{p}\) then \(p\) is prime. Also if \(p\) is prime then \((p-1)! \equiv -1 \pmod{p}\). This gives us a primality test that only requires multiplication and division among other things.

Why are we interested in it? Consider \((p-5)! \equiv -1 \pmod{p}\) for a very large prime number \(p\). Doing all of the multiplications to calculate the factorial could take a long time, we could have problems storing the number in a common integer type and might have to move to a slower BigInt class, etc… However if we know that \(p \) is prime then we know by Wilson’s theorem that \((p-1)! \equiv -1 \pmod{p}\) we can then apply the deffinition of the factorial function to this to get \((p-1)(p-2)(p-3)(p-4)(p-5)! \equiv -1 \pmod{p}\). But now we require division! How can we be sure that these numbers have an inverse modulo \(p\)? Well it’s simple to prove that for a prime number \(p\) every integer less than \(p\) is coprime to \(p\). And \(p > (p-1)\) etc… so we know that \(p-1,p-2,p-3,p-4\) all have inverses and this if we can find these inverses we can find \((p-1)(p-2)(p-3)(p-4)(p-5)! \equiv -1 \pmod{p}\) with only a handful of multiplications compared to the enourmous number required to calculate the actual factorial.

I’ll return to this subject in another post where I’ll show how to find these modular inverses.

Problem #91 says

The points \(P (x_1, y_1)\) and \(Q (x_2, y_2)\) are plotted at integer co-ordinates and are joined to the origin, \(O(0,0)\), to form \(ΔOPQ\).
There are exactly fourteen triangles containing a right angle that can be formed when each co-ordinate lies between 0 and 2 inclusive; that is,
\( 0 \leq x_1, y_1, x_2, y_2 \leq 2 \).
Given that \( 0 \leq x_1, y_1, x_2, y_2 \leq 50 \), how many right triangles can be formed?

A good way to tackle this problem is with good old Pythagoras’s equation. \(a^2 + b^2 = c^2\) but we need to avoid floating point operations which might screw up our results! The best way to do this is to avoid taking square roots because we know all of our points are located on integral spots on the grid. So we create little triangle objects that hold the origin and the two other points and we check the Pythagoras equation and if it passes we add it to our collection after checking for duplicates. Simply iterate through all of the points and count the triangles in out collection to get our answer.

So somethings that could be changed to make it faster. This solution creates a lot of objects and doesn’t really need them; this whole program could be rewritten to run without the OriginTri and Point classes but it’s easier to think about with them in my opinion. Removing the object creation however could shave quite a bit off of the run time. Solution in c# and runs in 13 seconds or so on my machine.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;
namespace Problem91
{
    class Program
    {
        static void Main(string[] args)
        {
            Stopwatch sw = new Stopwatch();
            sw.Start();
            OriginTri testcase;
            HashSet<string> triangles = new HashSet<string>();
            int max = 50;
            int count = 0;
            for (int x1 = 0; x1 <= max; x1++)
            {
                for (int y1 = 0; y1 <= max; y1++)
                {
                    for (int x2 = 0; x2 <= max; x2++)
                    {
                        for (int y2 = 0; y2 <= max; y2++)
                        {
                            testcase = new OriginTri(new Point(x1, y1), new Point(x2, y2));
                            if (testcase.IsRightTri() && !triangles.Contains(testcase.ToString()))
                            {
                                count++;
                                triangles.Add(testcase.ToString());
                                triangles.Add((new OriginTri(new Point(x2, y2), new Point(x1, y1))).ToString());
                            }
                        }
                    }
                }
            }
            sw.Stop();
            Console.WriteLine(count);
            Console.WriteLine(sw.Elapsed);
            Console.ReadLine();
        }
    }
    class OriginTri
    {
        Point origin = new Point(0,0);
        Point p1;
        Point p2;
        public OriginTri(Point p, Point q)
        {
            p1 = p;
            p2 = q;
        }
        public bool IsRightTri()
        {
            if (p1.ToString() == p2.ToString())
                return false;
            if (p1.ToString() == origin.ToString() || p2.ToString() == origin.ToString())
                return false;
        
            int S1 = p1.SquaredDistToPoint(origin);
            int S2 = p2.SquaredDistToPoint(origin);
            int S3 = p1.SquaredDistToPoint(p2);

            if (S1 + S2 == S3)
                return true;
            if (S2 + S3 == S1)
                return true;
            if (S3 + S1 == S2)
                return true;
            return false;
        }
        public override string ToString()
        {
            if (p1.MagFromZero() < p2.MagFromZero())   
                return "(0,0)" + p1.ToString() + p2.ToString();
            if (p1.MagFromZero() >= p2.MagFromZero())
                return "(0,0)" + p2.ToString() + p1.ToString();
            return "error";
        }
    }
    class Point
    {
        public int X;
        public int Y;
        public Point(int x, int y)
        {
            this.X = x;
            this.Y = y;
        }
        public int SquaredDistToPoint(Point p)
        {
            int tmpy = p.Y - this.Y;
            int tmpx = p.X - this.X;

            return (tmpx * tmpx) + (tmpy * tmpy);
        }
        public double MagFromZero()
        {
            return Math.Sqrt((double)(X * X) + (double)(Y * Y));
        }
        public override string ToString()
        {
            return "(" + X.ToString() + "," + Y.ToString() + ")";
        }
    }
}

Problem #125 says

The palindromic number 595 is interesting because it can be written as the sum of consecutive squares: \(6^2 + 7^2 + 8^2 + 9^2 + 10^2 + 11^2 + 12^2\).

There are exactly eleven palindromes below one-thousand that can be written as consecutive square sums, and the sum of these palindromes is 4164. Note that \(1 = 0^2 + 1^2\) has not been included as this problem is concerned with the squares of positive integers.

Find the sum of all the numbers less than \(10^8\) that are both palindromic and can be written as the sum of consecutive squares.

So the first thing we need to notice is that numbers like \(121 = 11^2\) don’t count for this problem because the number must be created from two or more squares. Next we need a formula for generating sums of squares and we can do some math and end up with \(\frac{(n(1+n)((2n)+1))}{6}\) and with mathematical induction we can show that this works for all natural numbers. But that formula starts at 1 and goes to n squaring each term; what about numbers like \(6^2 + 7^2 + 8^2 + 9^2 + 10^2 + 11^2 + 12^2\)? Well, if we generate all of the sums from 1 to n, we can subtract smaller sums to get the terms that don’t start at 1. Doing this for summation with all of the summations smaller than the first gives us the missing terms. From then we just need to check for palindromes and add the filtered terms.

To accomplish this we create a helper class that keeps track of the sum and where the sum starts and stops. For example let \(P(x,y) := x^2 + (x+1)^2 + … + y^2\) with \(x < y\) and \(x \neq y\). Then \(P(1,4) = 1^2 + 2^2 + 3^2 + 4^2 = 30\) and \(P(1,8) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 + 6^2 + 7^2 + 8^2 = 204\) and from this we can make \(P(5,8) = P(1,8) – P(1,4) =174\) just wash, rinse and repeat to find all of the valid sums. To make sure a sum is valid we need its starting point minus its ending point to be at least 2. We take these and check for the palindrome property and from there it’s an addition problem.

Solution provided in C#/mono and runs in 5 seconds or so with most of the time being spent checking to make sure we’re not adding duplicates. This program could potentially be made much faster if we were to make the math slightly more complicated to avoid creating all of these summation helper objects. But its fast enough for this problem.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace Problem125
{
    class Program
    {
        static void Main(string[] args)
        {
            Stopwatch sw = new Stopwatch();
            sw.Start();
            List<SumHolder> Possibles = new List<SumHolder>();
            long max = 100000000000L;//100000000L;
			
			long t = 1L;
			while (SumOfCSquares(t-1L) <=max)
			{
				SumHolder s = new SumHolder(SumOfCSquares(t),t,1L);
				Possibles.Add(s);
				t++ ;
			}
			List<SumHolder> Generated = new List<SumHolder>();
			foreach (SumHolder s in Possibles)
			{
				foreach (SumHolder s2 in Possibles)
				{
					if (s.Sum > s2.Sum)
					{
						if (s.Sum != s2.Sum)
						{
							SumHolder temp = new SumHolder(s.Sum -s2.Sum, s.UpperLimit,s2.UpperLimit);
							Generated.Add(temp);
						}
					}
				}
			}
			long sum = 0L;
			int count = 0;
            
                        HashSet<long> Duplicates = new HashSet<long>();
			foreach (SumHolder s in Possibles)
			{
				if (Duplicates.Contains(s.Sum))
				    continue;
				if (s.UpperLimit - s.LowerLimit >= 1)
				{
					if (s.Sum <= 100000000L)
					{
						if (Palindrome(s.Sum))
						{
							sum += s.Sum;
							Console.WriteLine("Found :: "+s.Sum.ToString());
							Duplicates.Add(s.Sum);
							count++;
						}
					}
				}
			}
			foreach (SumHolder s in Generated)
			{
				if (Duplicates.Contains(s.Sum))
				    continue;
				if (s.UpperLimit - s.LowerLimit >= 2)
				{
					if (s.Sum <= 100000000L)
					{
						if (Palindrome(s.Sum))
						{
							sum += s.Sum;
							Console.WriteLine("Found :: "+s.Sum.ToString());
							Duplicates.Add(s.Sum);
							count++;
						}
					}
				}
			}
            sw.Stop();
            Console.WriteLine(sum);
			Console.WriteLine(count);
            Console.WriteLine(sw.Elapsed);
            Console.ReadLine();

        }
        static bool Palindrome(long n)
        {
            string s = n.ToString();
            string rs = Reverse(n.ToString());
            return rs == s;
        }
        static string Reverse(string str)
        {
            char[] array = str.ToCharArray();
            Array.Reverse(array);
            return new string(array);
        }
        static long SumOfCSquares(long n)
        {

            return (n*(1+n)*(1+(2*n)))/6L;

        }

    }
	class SumHolder
	{
		public long Sum { get; set; }
		public long UpperLimit { get; set; }
		public long LowerLimit { get; set; }
		public SumHolder(long s, long u, long l)
		{
			this.Sum = s;
			this.UpperLimit = u;
			this.LowerLimit = l;
		}
	}
	
}

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