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

Addition Chains are paths or little sections of an Addition Tree. To build the tree start with the number 1. To this root add a new node for each ancestor + itself such that the value of the new node is the parent node + ancestor. Thus for the node 1 the new children nodes would be a single node of value 2. The next level then becomes 3(2+1) and 4(2+2). If we then follow the node 3 down to the next level it’s children become 4 (3+1) 5(3+2) & 6(3+3) and we can construct the children of the node 4 with the same technique; we can follow this tree right into infinity if we so choose. To obtain an addition chain we choose a node and get all of its ancestors and that’s the chain. So for the node 6 in our example above we get 6(3+3),3(2+1),2(1+1),1 and this path makes the Addition chain. This algorithm creates duplicate nodes some with longer chains that previously occurring nodes with a give value and the problem of creating optimal addition chains is NP-complete. The below implementation uses this naïve algorithm rather than try and create a single optimal path.

It can be shown that creating a tree containing only optimal paths is impossible with the example triplet 43, 77 and 149; having the first two of these numbers with an optimal path makes it impossible to have an optimal path for the third. Thus a tree can contain many paths and duplicates or not contain optimal paths. In the code below if AllowDuplicates is defined the tree is much smaller but there is no longer a guarantee of an optimal path being found. The library below allows duplicates and when you ask for a node of a given value it gets all nodes of that value and looks for the node with the fewest ancestors which would make that node the optimal path; if there is multiple nodes with the same number of ancestors it defaults to the first path it found.

Now, what good are these things you ask? Well we can use them for Addition Chain exponentiation. Say we wish to calculate \(n^{8}\). Well, we could multiply \(n * n * n * n * n * n * n * n\) and get our answer; this requires 7 multiplications. However, we can create an addition chain for 8 which would look like \(n^{1}, n^{2} (1*1), n^{4} (2*2), n^{8} (4*4)\). This can be used for exponentiation because of the rules for exponents where if we multiply two exponents with a common base we just add the exponent values. This chain method allows us to calculate \(n^{8}\) with only 3 multiplications; which for a small value n isn’t that impressive. But for a larger value this can make quite a difference. However, there is a trade off with the Addition Chain method using more memory and being more complicated to write.

So where is this all going? Well I found them interesting and wanted to share some code to build them. But more practically there is a Project Euler problem (Problem #122) that asks us to find a group of optimal chains.

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

namespace AdditionChain
{
    class AdditionChain
    {
        public AdditionChainNode Root
        {
            get;
            set;
        }
        private List<int> _knownValues = new List<int>();
        private List<AdditionChainNode> _nodeListing = new List<AdditionChainNode>();
        public AdditionChain(int rootVal)
        {
            this.Root = new AdditionChainNode(rootVal);
            _nodeListing.Add(this.Root);
            _knownValues.Add(this.Root.Value);
        }
        public void AddChild(AdditionChainNode parent, int value)
        {
			#if AllowDuplicates
            if (_knownValues.Contains(value))
                return;
            else
            {
			#endif
                AdditionChainNode n = new AdditionChainNode(value);
                n.Parent = parent;
                _knownValues.Add(value);
                _nodeListing.Add(n);
                parent.Children.Add(n);
			#if AllowDuplicates
            //}
			#endif
        }
        public int NodeLevel(AdditionChainNode start)
        {
            if (start == null)
                return int.MaxValue;
            int Val = 0;
            while (start.Parent != null)
            {
                start = start.Parent;
                Val++;
            }
            return Val;
        }
        public AdditionChainNode GetNodeByValue(int val)
        {
            while (!_knownValues.Contains(val))
                CreateNextLevel();
            AdditionChainNode m = null;
            foreach(AdditionChainNode n in _nodeListing)
            {
                //do the cheap check first
                if (n.Value == val)
                {
                    if (NodeLevel(n) < NodeLevel(m))
                        m = n;
                }
            }

            return m;
        }
        public void CreateNextLevel()
        {
            List<AdditionChainNode> lowest = GetLowestLevelNodes();
            foreach (AdditionChainNode n in lowest)
            {
                List<int> ancestors = GetAncestorValues(n);
                foreach (int i in ancestors)
                {
					#if AllowDuplicates
                    //if (!_knownValues.Contains(n.Value + i))
                    //{
					#endif
                        AddChild(n, n.Value + i);
					#if AllowDuplicates
                    //}
					#endif
                }
            }
        }

        private List<AdditionChainNode> GetLowestLevelNodes()
        {
            if (Root.Children.Count == 0)
                return new List<AdditionChainNode>() { Root };
            else
            {
                List<AdditionChainNode> Vals = new List<AdditionChainNode>();

                foreach (AdditionChainNode n in _nodeListing)
                {
                    if (n.Children.Count == 0)
                        Vals.Add(n);
                }

                return Vals.OrderBy(i=> i.Parent.Value).ToList();
            }
        }
        public List<int> GetAncestorValues(AdditionChainNode start)
        {
            if (start.Parent == null)
                return new List<int>(){start.Value};
            List<int> Vals = new List<int>();

            AdditionChainNode p = start;
            Vals.Insert(0, p.Value);

            while (p.Parent != null)
            {
                p = p.Parent;
                Vals.Insert(0, p.Value);
            }

            return Vals;
        }
    }
    class AdditionChainNode
    {

        public List<AdditionChainNode> Children
        {
            get;
            set;
        }
        public AdditionChainNode Parent
        {
            get;
            set;
        }
        public int Value
        {
            get;
            set;
        }
        public override string ToString()
        {
            return this.Value.ToString();
        }
        public AdditionChainNode(int v)
        {
            this.Value = v;
            this.Children = new List<AdditionChainNode>();
        }
    }
}

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

Problem #131 says:

There are some prime values, p, for which there exists a positive integer, n, such that the expression \(n^{3} + n^{2}p\) is a perfect cube.

For example, when \(p = 19, 8^{3} + 8^{2}×19 = 12^{3}\).

What is perhaps most surprising is that for each prime with this property the value of n is unique, and there are only four such primes below one-hundred.

How many primes below one million have this remarkable property?

This is a problem that would drive a person insane if they tried to brute force a solution; even though the search area is (relatively for a Project Euler problem) small. This problem has a very elegant solution as long as you know how to break down the problem to it’s base components. Writing the problem out gives us \(n^{3} + n^{2}p = y^{3}\) which we can rearrange to look like \(n^{2}(n+p) = y^{3}\). Playing with a few values leads us to the conclusion that both \(n^{2}\) and \((n+p)\) have to be cubes. This makes p a difference of cubes such that \( (a+1)^{3} – a^{3} = p \)

So we take this new knowledge and figure out when \( (a+1)^{3} – a^{3} \) becomes greater than 1,000,000 and find that this happens around a = 577. This means that we just put the numbers from 1 to 577 into our difference of squares equation and count which ones return a prime.

Solution provided in c#/mono runs in about .038 seconds.

using System;
using System.Collections.Generic;
namespace Problem131
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			List<int> primes = GeneratePrimes(1000);
			int count = 0;

			for (int i = 1; i < 577; i++)
			{
				int check = ((i+1)*(i+1)*(i+1)) - (i*i*i);

				if (IsPrime(check,primes))
				{
					count++;
				}
			}
			Console.WriteLine(count);
			Console.ReadLine();
		}
		static bool IsPrime(int i,List<int> primes)
		{
			foreach(int p in primes)
			{
				if (i%p == 0 && i!= p)
					return false;
			}

			return true;
		}
		static List<int> GeneratePrimes(int num)
        {
            List<int> RetValue = new List<int>();
            RetValue.Add((int)2);
            RetValue.Add((int)3);
            int Stepper = 5;
            int Check = 1;
            while (Stepper <= num)
            {
                foreach (int 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;
        }
	}

}

So according to my website stats people are ending up on this site after looking for the keywords “check if triangle contain the origin” (yeah I know that’s not good English but it’s what people enter into the search engine), no doubt looking for help with Project Euler Problem #102! (warning that link contains a solution to the problem; don’t click if you want to solve it yourself). However maybe they just want something to check a point and a triangle and see if the triangle encompasses the point. So I went a head and extracted my function for testing whether a triangle encompasses a point from my Project Euler solution and translated it to f# and cleaned it up a bit.

This function uses the half-plane test to check if a triangle contains a point (said point could be the origin if one wants.) this method could be extended for n-polygon shapes but as provided it’s limited to triangles.

Code provided in f#.

#light
open System

type Point(x:float, y) =
  member p.X = x
  member p.Y = y

let Origin = new Point(0.0,0.0)

let Sign (p1:Point) (p2:Point) (p3:Point) =
    (p1.X - p3.X) * (p2.Y - p3.Y) - (p2.X - p3.X) * (p1.Y - p3.Y)

let CheckPoint pt p1 p2 p3 =
    let test1 = (Sign pt p1 p2) < 0.0
    let test2 = (Sign pt p2 p3) < 0.0
    let test3 = (Sign pt p3 p1) < 0.0
    ((test1 = test2) && (test2 = test3))

let A = new Point(-340.0,495.0)
let B = new Point(-153.0,-910.0)
let C = new Point(835.0,-947.0)

let X = new Point(-175.0,41.0)
let Y = new Point(-421.0,-714.0)
let Z = new Point(574.0,-645.0)

Console.WriteLine (CheckPoint Origin A B C)
Console.WriteLine (CheckPoint Origin X Y Z)