Project Euler Problem #86!

July 20th, 2011
by Serinox

Problem #86 says:

A spider, S, sits in one corner of a cuboid room, measuring 6 by 5 by 3, and a fly, F, sits in the opposite corner. By travelling on the surfaces of the room the shortest “straight line” distance from S to F is 10 and the path is shown on the diagram.

However, there are up to three “shortest” path candidates for any given cuboid and the shortest route is not always integer.

By considering all cuboid rooms with integer dimensions, up to a maximum size of M by M by M, there are exactly 2060 cuboids for which the shortest distance is integer when M=100, and this is the least value of M for which the number of solutions first exceeds two thousand; the number of solutions is 1975 when M=99.

Find the least value of M such that the number of solutions first exceeds one million.

I thought this problem was going to be more interesting than it turned out being. Basically we’re looking for triplets (a,b,M) where 1<=a<=b<=M that make the equation \((a+b)^{2}+n^{2}\) a perfect square. We just do this for every value M=1-?? and count the number of solutions for that M value until the sum of all previous solution counts equals one million. Basically we take M = 3 which has a solution count of 3 and is the first none zero solution count and M = 4 with a solution count of 1; We add the solutions counts for M = 1-4 which are {0,0,2,1} and check to see if the sum of that series is greater than one million. If that sum isn’t one million we calculate the next value of M and add that to the series and check the series sum again; wash, rinse, repeat.

This problem became a lot less interesting when I found the formula for the shortest point and googled it and came back with two OEIS sequences (A143714,A143715) which made the solution trivial to write because it showed the method of building on the previous M values. On the other hand trivial to program solutions mean a lot less time spent frustrated trying to get the program to work.

Solution provided in c#/mono and runs in ~30ish seconds. Not the fastest but it works.

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

namespace Problem86
{
	class MainClass
	{
		public static List<int> A143714 = new List<int>();

		public static void Main (string[] args)
		{
			A143714.Add(0);
			A143714.Add(0);

			while(A143714.Sum() < 1000000)
			{
				int n = A143714.Count+1;

				int count = 0;
				for (int b = 1;b<=n;b++)
				{
					for (int a = 1;a<=b;a++)
					{
						if (IsPerfectSquare(((a+b)*(a+b))+(n*n)))
							count++;
					}
				}

				A143714.Add(count);
			}
			Console.WriteLine(A143714.Count);
		}

		public static bool IsPerfectSquare(int input)
		{
		    var sqrt = Math.Sqrt((double)input);
		    return Math.Ceiling(sqrt) == Math.Floor(sqrt);
		}

	}
}

With this solution I have only 11 problems left between 1-100 and only 39 left until I hit level 4!

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

Project Euler Problem #80!

July 18th, 2011
by Serinox

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.

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

Project Euler Problem #68!

July 15th, 2011
by Serinox

Problem #68 boils down to:

… By concatenating each group it is possible to form 9-digit strings; the maximum string for a 3-gon ring is 432621513.

Using the numbers 1 to 10, and depending on arrangements, it is possible to form 16- and 17-digit strings. What is the maximum 16-digit string for a “magic” 5-gon ring?

Some rules apply. First all of the spars must add up to the same number and you need to pay close attention to how you build the final answer to put into the answer box or you’ll get told your correct polygon is wrong. (See the Project Euler page for details and an example)

This is a question that is easier solved on paper with a pencil than any program. Just a few simple assumptions will make solving this problem easier: 1) The outer ring has a higher place value compared to it’s inner ring counterparts. 2) You want a large number so it’s easy see you want the numbers 6-10 on the outer rings. 3) We want a 16 digit number which means that the 10 must be on the outer ring else it gets counted twice and you’ll end up with a 17 digit polygon.

There is no solution presented for this problem; however with the above hints it’s very easy to build the answer.

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

Project Euler Problem #342!

July 5th, 2011
by Serinox

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;
}

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

Project Euler Problem #95!

June 28th, 2011
by Serinox

Problem #95 says:

The proper divisors of a number are all the divisors excluding the number itself. For example, the proper divisors of 28 are 1, 2, 4, 7, and 14. As the sum of these divisors is equal to 28, we call it a perfect number.

Interestingly the sum of the proper divisors of 220 is 284 and the sum of the proper divisors of 284 is 220, forming a chain of two numbers. For this reason, 220 and 284 are called an amicable pair.

Perhaps less well known are longer chains. For example, starting with 12496, we form a chain of five numbers:

12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → …)

Since this chain returns to its starting point, it is called an amicable chain.

Find the smallest member of the longest amicable chain with no element exceeding one million.

An extremely easy problem to solve despite the fact that it’s quite far down the list when sorted by difficulty. We can instantly stop generating the chain if it exceeds one million at any point in the chain. Also we need to check to make sure the chain ends up where we started or it isn’t valid for the problem; this means that we need to make sure we don’t get stuck in an infinite loop where the chain starts over and we miss it. Solution is written using a simple sum of divisors function that simply brute forces the sum instead of a more elaborate prime factor function; I don’t think the prime factor version would be too much faster because the magnitude of the numbers (under one million) just isn’t enough to warrant such an effort.

Solution written in mono/c# with a helping hand from the linq name-space (mainly for list.min and list.max functions without having to write my own) runs in about 6 seconds on a modern machine. Interestingly enough this could be made parallel if one so desired with a few simple modifications but it runs fast enough already.

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

namespace Problem95
{
	class MainClass
	{
		public static void Main (string[] args)
		{

			List<long> chain = new List<long>();
			List<long> tmp = new List<long>();

			for (long i = 1; i <= 100000; i++)
			{
				tmp = GenerateChain(i);
				if (tmp.Count() > chain.Count())
				{
					List<long> ch = tmp.Where(x=> x > 1000000).ToList();
					if (ch.Count() == 0)
						chain = GenerateChain(i);
				}
			}
			Console.WriteLine(chain[0].ToString() + " :: " + chain.Min().ToString());
		}
		public static List<long> GenerateChain (long num)
		{
			List<long> chain = new List<long>();
			long previous = num;
			while (!chain.Contains(previous) && previous < 1000000L)
			{
				chain.Add(previous);

				previous = NextInChain(previous);
			}

			if (previous != chain[0])
				chain.Clear();

			if (previous > 1000000L)
				chain.Add(9999999L);

			return chain;
		}
		public static long NextInChain (long num)
		{
			return ProperDivisors(num).Sum();
		}
		public static List<long> ProperDivisors (long num)
		{
			long max = (long)Math.Ceiling(Math.Sqrt((double)num));

			List<long> RetValue = new List<long>();
			RetValue.Add(1L);
			if (num % max == 0d)
				RetValue.Add(max);

			for(long i = 2; i < max; i++)
			{
				if (num % i == 0L)
				{
					RetValue.Add(i);
					RetValue.Add(num/i);
				}
			}

			return RetValue.Distinct().ToList();
		}
	}
}

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