Problem #122 says:

We shall define m(k) to be the minimum number of multiplications to compute \(n^{k}\); for example m(15) = 5.

For 1 ≤ k ≤ 200, find ∑ m(k).

This is a problem about finding Addition Chains that lead to the desired number in the fewest steps. Using the library code I posted yesterday the solution can be written as 5 lines of code (not counting main() and other fluff). If you compile the Addition Chain library so that it doesn’t allow duplicates then the program will give the wrong answer; This is because duplicate paths are needed to ensure that an optimal chain can be found.

Other than that the solution runs in about 45 seconds an decent computer. Solution provided in c#/mono.

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

namespace Problem122
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			AdditionChain ac = new AdditionChain(1);
			int answer = 0;
			for (int i = 1; i<= 200; i++)
				answer += ac.NodeLevel(ac.GetNodeByValue(i));
			Console.WriteLine(answer);
		}
	}
}

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 #98 asks us:

By replacing each of the letters in the word CARE with 1, 2, 9, and 6 respectively, we form a square number: 1296 = \(36^{2}\). What is remarkable is that, by using the same digital substitutions, the anagram, RACE, also forms a square number: 9216 = \(96^{2}\). We shall call CARE (and RACE) a square anagram word pair and specify further that leading zeroes are not permitted, neither may a different letter have the same digital value as another letter.

Using a 16K text file containing nearly two-thousand common English words, find all the square anagram word pairs (a palindromic word is NOT considered to be an anagram of itself).

What is the largest square number formed by any member of such a pair?

NOTE: All anagrams formed must be contained in the given text file.

So we break this down into several parts. 1) Anagrams have to be the same length. 2) Remove words we’ve already paired up or failed to find an anagram for. 3) Each letter must have 1 numerical value assigned to it and used consistently. 4) Likewise each number can only be assigned to one letter.

With those four tips we can structure our solution like this. Build all of the anagram pairs and remove the words we’ve paired up from the main word list. Build a list of all of the square numbers from 2 – 1000000000. For each pair take the first word and find a square that is the same length. Assign each letter to a number without repeats of either letters or numbers (each pairing must be unique in both directions). Map the letter->number combination onto the second word. Parse the mapped second word to a ulong and make sure that it matches the length of the square assigned to word one. If they are the same length check (this will mean that there is not leading zero’s etc) then check the list of square numbers to see if the number created by mapping the second number appears in the list. If the mapped second word is a square then all we have to check is which is higher either the square assigned to the first word or the square we generated with the second. Take the highest one and go to the next anagram word pair.

This approach makes it so we don’t have to check to see if the square number assigned to word one and the generated square are the same because the only way that could happen is if word one and word two are the same words. There are no duplicates in the word list so we can eliminate all checking for that condition. Also the way we move through the list makes it impossible to compare a word to itself meaning we don’t have to waste time checking that either.

Despite the code being a mess it runs far faster than it looks like it should. Solution provided in c#/mono and runs in ~1.5 seconds.

* It’s left as an exercise to the reader to figure out how to get the word list into the words variable. I don’t want to screw up search engines by posting 1786 random words.

using System;
using System.Collections.Generic;

namespace Problem98
{
	class MainClass
	{
		public static List<string> words = new List<string>();
		public static List<List<string>> AnagramWordPairings = new List<List<string>>();
		public static List<ulong> squares = new List<ulong>();
		public static void Main (string[] args)
		{

			for (ulong i = 2; i <= 31700; i++)
			{
				squares.Add(i*i);
			}
			Max("rat","tar");
			while (words.Count > 0)
			{
				string primary = words[0];
				List<string> anagramset = new List<string>();
				anagramset.Add(primary);

				for (int i = 1; i < words.Count;i++)
				{
					if (Anagrams(primary,words[i]))
						anagramset.Add(words[i]);
				}

				if (anagramset.Count ==1)
					words.Remove(primary);
				else
				{
					foreach (string s in anagramset)
						words.Remove(s);
					AnagramWordPairings.Add(anagramset);
				}

			}
			ulong max = 0;
			foreach(List<string> s in AnagramWordPairings)
			{
				ulong tmax = Max(s[0],s[1]);
				if (tmax > max)
					max = tmax;
			}
			Console.WriteLine(max);
		}
		public static ulong Max (string s1,string s2)
		{
			ulong max = 0;
			foreach(ulong sq in squares)
			{
				if (sq.ToString().Length == s1.Length)
				{
					char[] premap = s1.ToCharArray();
					int digits = sq.ToString().Length;
					Dictionary<char,char> map = new Dictionary<char, char>();
					foreach (char c in premap)
					{
						if (!map.ContainsKey(c) &&
                                                    !map.ContainsValue(sq.ToString()[s1.IndexOf(c.ToString())]))
							map.Add(c,sq.ToString()[s1.IndexOf(c.ToString())]);
					}
					string ts2 = (string)s2.Clone();
					foreach (char c in map.Keys)
					{
						ts2 = ts2.Replace(c,map[c]);
					}
					try
					{
					ulong tp = ulong.Parse(ts2);
					if (tp.ToString().Length != sq.ToString().Length)
						continue;
					if (squares.Contains(tp))
						max = Math.Max(sq,tp);
					}
					catch{}
				}
			}
			return max;
		}

		public static bool IsSquare (ulong i)
		{
			double sqrt = Math.Sqrt((double)i);
			return Math.Ceiling(sqrt) == Math.Floor(sqrt);
		}
		public static bool Anagrams (string s1, string s2)
		{
			if (s1.Length != s2.Length)
				return false;
			char[] c1 = s1.ToUpper().ToCharArray();
			char[] c2 = s2.ToUpper().ToCharArray();

			Array.Sort(c1);
			Array.Sort(c2);
			string ns1 = new string(c1);
			string ns2 = new string(c2);
			return (ns1==ns2);
		}
	}
}

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.

Problem #342 says:

Consider the number 50.
\(50^{2} = 2500 = 2^{2} × 5^{4}\), so \(φ(2500) = 2 × 4 × 5^{3} = 8 × 5^{3} = 2^{3} × 5^{3}.\)
So 2500 is a square and φ(2500) is a cube.

Find the sum of all numbers n, 1 < n < \(10^{10}\) such that \(φ(n^{2})\) is a cube.

This solution makes much use of the following property of the totient function \(φ(p^{a}) = p^{a-1} (p-1)\). With a little effort one can see that if you have the prime factorization of a number N you can calculate the prime factorization of the totient of \(n^{2}\) and if one has that factorization the to make sure that it is a cube is as simple as making sure that the power of each exponent in the factorization is 0 mod 3.

Thus for the example problem above we have the number n = 50. Which factors to \(2 * 5^{2}\). Using the above identity on the 5 factor we get \(φ(5^{2}) = 5^{(2*2)-1} (5-1)\) which when we factor any nonfactored terms comes out to be \(2^{2} * 5^{3}\) which is then multiplied by the result of the same property applied to the factor 2 gets us the result \( 2^{3} * 5^{3}\). We can see that this is a cube because all factors are raised to a multiple of three.

With that basic set my solution is relying on factoring numbers and as most people know that’s rather difficult to do fast. Especially when memorizing previously factored numbers becomes impossible due to the size of the problem set. However this is a massively parallel operation for each value of n. So I took this as an opportunity to learn some OpenMPI and see what I could do with at bunch of old computers. However this means that solution only runs in a reasonable time if you have a lot of hardware backing it and a network for the various computers to talk over. It’s actually quite easy to setup a cluster from old computers but maybe not for the faint of heart or those unused to Linux.

Solution written in c++ using g++ 4.5.1 and OpenMPI 1.4.1. Runtime depends on how many processors you give the program; but I will say that it’s far from the fastest solution unless you have an account on a TOP500 supercomputer.

P.S. This program divides work with a fairly simple algorithm and works best with a multiple of 10 for the number of processes created by mpiexec. If you don’t use a multiple of 10 then it may skip numbers or come up short of the actual answer.

#include <iostream>
#include <vector>
#include <map>
#include <string.h>
#include <time.h>
#include <algorithm>
#include <mpi.h>
#include <math.h>

using namespace std;

vector<unsigned long> get_primes(unsigned long max){
    vector<unsigned long> primes;
    char *sieve;
    sieve = new char[max/8+1];
    // Fill sieve with 1
    memset(sieve, 0xFF, (max/8+1) * sizeof(char));
    for(unsigned long x = 2; x <= max; x++)
        if(sieve[x/8] & (0x01 << (x % 8))){
            primes.push_back(x);
            // Is prime. Mark multiplicates.
            for(unsigned long j = 2*x; j <= max; j += x)
                sieve[j/8] &= ~(0x01 << (j % 8));
        }
    delete[] sieve;
    return primes;
}

vector<unsigned long> primes = get_primes(100000);

map<unsigned long, unsigned long> Factor (unsigned long n)
{
    map<unsigned long, unsigned long> Factors;
    unsigned long sqrtn = sqrt(n);
    for (unsigned int i = 0; i < primes.size(); i++)
    {
        if (primes[i] > sqrtn)
        {
            break;
        }
        else
        {
            if (n % primes[i] == 0)
            {
                Factors[primes[i]] = 0ul;
                while(n % primes[i] == 0)
                {
                    Factors[primes[i]] = Factors[primes[i]] + 1ul;
                    n = n/primes[i];
                }
            }
        }
    }
    if (n != 1)
        Factors[n]=1;
    if (Factors.size() == 0)
    {
        Factors[n] = (unsigned long)1;
    }
    return Factors;
}

bool IsTotientCube(unsigned long n)
{
    map<unsigned long, unsigned long> factors = Factor(n);
    map<unsigned long, unsigned long>::iterator pos;
    map<unsigned long, unsigned long> FactorsOfSquaredTotient;
    for (pos = factors.begin(); pos != factors.end(); ++pos)
    {

        map<unsigned long, unsigned long> pm1factors = Factor((pos->first)-1);
        map<unsigned long, unsigned long>::iterator pm1pos;
        for (pm1pos = pm1factors.begin(); pm1pos != pm1factors.end(); ++pm1pos)
        {

            if (pm1pos->first != 1ul)
                FactorsOfSquaredTotient[pm1pos->first] += pm1pos->second;
        }
        FactorsOfSquaredTotient[pos->first] += (2* (pos->second)) -1;
    }
    map<unsigned long, unsigned long>::iterator possquaretotient;
    for (possquaretotient = FactorsOfSquaredTotient.begin();
         possquaretotient != FactorsOfSquaredTotient.end();
         ++possquaretotient)
    {

        if (possquaretotient->second % 3 != 0)
        {
            return false;
        }

    }
    return true;
}

int main(int argc, char *argv[])
{
    int numprocs, rank, namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Get_processor_name(processor_name, &namelen);

    unsigned long count = 0ul;
    unsigned long max = 10000000000;

    unsigned long workload = max / numprocs;

    unsigned long responses[numprocs];

    for (unsigned long i = (workload * rank); i < ((workload * rank) + workload);++i)
    {
        if (IsTotientCube(i) == true)
        {
            count += i;
            cout << i << endl;
        }
    }

    MPI_Barrier(MPI_COMM_WORLD);

    MPI_Gather(&count,1,MPI_UNSIGNED_LONG,&responses,1,MPI_UNSIGNED_LONG,0,MPI_COMM_WORLD);
    unsigned long answer = 0ul;
    for (int i = 0; i < numprocs;++i)
    {
        answer += responses[i];
    }

    cout << "Rank "<< rank << " Found #"<< count << endl;

    if (rank == 0)
        cout << "Answer "<< answer << endl;
    MPI_Finalize();
    return 0;
}