Project Euler 29: How many distinct terms are in the sequence generated by a^b for 2 ≤ a ≤ 100 and 2 ≤ b ≤ 100?

Written by on 26 February 2011

Topics: Project Euler

Problem 29 of Project Euler reads

Consider all integer combinations of ab for 2 <= a <= 5 and 2 <= b <= 5:

22=4, 23=8, 24=16, 25=32
32=9, 33=27, 34=81, 35=243
42=16, 43=64, 44=256, 45=1024
52=25, 53=125, 54=625, 55=3125
If they are then placed in numerical order, with any repeats removed, we get the following sequence of 15 distinct terms:

4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125

How many distinct terms are in the sequence generated by ab for 2 <= a <= 100 and 2 <= b <= 100?

It is fairly trivial to realise that the solution is somewhere below 99*99 = 9801, using combinatorics. Based on the amount of possible combinations it should also be fairly easy to implement a brute force solution in C#. That is the first solution we will elaborate on, but I have a pen and paper method as well, where we will use logic to deduce the solution.

Brute Force solution

I intend to use the built in list functionality to hold all the found solutions, and check if newly generated solutions are already found. So the brute force code is fairly simple. I will compare two different parts of the collection framework in C#; the List and the SortedSet.

The code consists of two loops each counting from 2 to 99 for a and b, and if the result does not exist in the list it is added. The code for the list implementation looks like

List set = new List();

for (int a = 2; a <= 100; a++) {
    for (int b = 2; b <= 100; b++) {
        double result = Math.Pow(a, b);
        if (!set.Contains(result)) {
            set.Add(result);
        }
    }
}

The code for the SortedSet implementation is simpler since it cannot hold duplicate entries, so the Contains check is implied when trying to add to the set. The code looks like

SortedSet set = new SortedSet();

for (int a = 2; a <= 100; a++) {
    for (int b = 2; b <= 100; b++) {
        double result = Math.Pow(a, b);
        set.Add(result);
     }
 }

So no significant changes in the code apart from using another list implementation. However, when running the two pieces of code the difference is evident

The number of distinct terms are 9183
Solution took 814 ms
The number of distinct terms are 9183
Solution took 6 ms

with the SortedSet implementation being 100 times faster. I wont go into details about the differences as I have difficulties figuring out exactly how they are implemented. However, one thing I can say is that the SortedSet takes longer to insert items into, but the Contains method is significantly faster. The learning to take away from this, is that the implementation of the data structure matters for execution speed.

The Pen and Paper Way

Although the original meaning of solving the exercises was to implement a solution in C#. As mentioned earlier the maximum number of solutions is 99*99 = 9801, so lets find the duplicated and subtract them from the maximum by hand.

Suppose a is a square of the smaller a, but not a square of a square. In that case we have duplicates when b is between 2-50, a total of 49 times. This happens for the numbers 4, 9, 25, 36, 49, 100 since 4 is a square of 2, and thus 42 = (22)2 = 22*2, and so on.
So 6 a’s which each has 40 duplicates. A total of 6*49 = 294 duplicates.

The reason 16 and 81 is not mentioned, is the condition that we need a number which is not a square of a square.

Same deduction can be made for cubes.  Suppose a is a cube of a smaller a but not cubes of squares. Then we have duplicates when b is between 2 and 33. With the same deduction as above; 43 = (23)2 = 23*2. For b = 34,36,38,…, 66 we have duplicates of the smaller a raised to 3(b/2). This gives us another 17 duplicates. So each cube gives us 32+17=49 duplicates.

This happens for the numbers; 8, 27 which have 32 duplicates each, which is a total of 98 duplicates. The reason 64 is left out, is that it is a cube of a square, and thus will be covered later.

Now we get the pattern, so lets check numbers which are the 4th power of a smaller a or squares of squares if you like.  We have two squares of squares 16 and 81. Whenever b is between 2 and 49, we have a duplicate for the square root of a raised to b*2. And when b is 51, 54, 57, … , 75 we have duplicates of a^(3/4)  raised to b*4/3. An example of that is 1651 = (163/4)51*4/3 = 868.
This gives us a total of 2*(49+16) = 116 duplicates.

The same analysis can be made when a is the 5th power of a smaller a. This happens for the number 32.  It has duplicates when b is between 2 and 20.  When b is 22, 24, 26, … , 40 we have duplicates for (a2/5)b*5/2, a total of 10 duplicates. When b is 21, 27, 33, 39, 42, 45, 48, 51, 54, 57 , 60 we have duplicates for (a3/5)b*5/3, a total of 11 duplicates. Please note the serie is irregular as some 24, 30 and 36 are covered by previous case. And last but not least we have duplicates for (a4/5)b*5/4 whenever b is 44, 52, 56, 64, 68, 72, 76, 80. A total of 8 duplicates. So the number 32 produces 19+10+11+8 = 48 duplicates.

When a is the 6th power of a smaller a, there are duplicates for the square root of a raised to 2b. This happens for b between 2 and 50.  The we have duplicates for (a4/6)b*6/4 covering 52, 54, 56, … , 66 and we have duplicates for (a5/6)b*6/5 covering 55, 65, 70, 75 , 80. A total of  62 duplicates.

Since 27 = 128 we do not need to check higher powers.

When we sum up the identified duplicates we end up with 294+98+116+48+62=618. Subtract this from 9801 which is the number of possible solutions and we get the result 9801-618 = 9183.

Source code

As usual you can find the source code for the implemented solutions. If you have comments, other solutions or questions, I would be delighted to discuss with you, so post a comment on the blog.

14 Comments For This Post I'd Love to Hear Yours!

  1. Anonymous says:

    Very nice.

  2. Leo says:

    i have this code ! but i really don’t know WHY IT DOESN’t Work !!
    can someone really help ?
    /////////////////////////////////
    #include
    #include
    #include
    using namespace std;
    int A[100];
    vector prime;
    vector<vector<vector<pair > > > mother;
    vector<pair > magh(int,int),temp;
    int vecmogh(vector<pair > , vector<pair > ),gcd(int ,int ),ans,test;
    using namespace std;

    int main()
    {
    cin>>test;
    for (int i=0;i<100;i++)
    {
    A[i]=i;
    }

    for (int i=2;i<50;i++)
    {
    if (A[i])
    {
    for (int j=2;j99)
    break;

    A[A[i]*j]=0;
    }
    }
    }
    mother.resize(101);
    for (int i=0;i<=100;i++)
    {
    mother[i].resize(101);
    }
    for (int i=2;i<100;i++)
    {
    if (A[i])
    prime.push_back(A[i]);
    }

    for (int i=2;i<=100;i++)
    {
    for (int j=2;j<=100;j++)
    {
    mother[i][j]=magh(i,j);
    }
    }

    for (int i=2;i<=test;i++)
    {
    for (int j=2;j<=test;j++)
    {
    for (int h=2;h<=test;h++)
    {
    if (h==i)
    continue;

    if (!(gcd(i,h)==h || gcd(i,h)==i))
    continue;

    for (int q=2;q<=test;q++)
    {
    if (vecmogh(mother[i][j],mother[h][q])==1)
    {
    ans++;
    cout<<" i : "<<i<<"\t j : "<<j<<"\t h : "<<h<<"\t q : "<<q<<endl;

    }
    }
    }
    }
    }
    ans/=2;
    ans=(test-1)*(test-1)-ans;
    cout<<endl<>ans;
    }

    vector<pair > magh(int key,int tavan)
    {
    int counter;
    vector<pair > ans;
    for (int i=0;i<prime.size();i++)
    {
    if (key==1)
    break;

    if (key%prime[i]==0)
    {
    counter=0;

    while(key % prime[i]==0)
    {
    key/=prime[i];
    counter++;
    }
    counter*=tavan;
    ans.push_back(make_pair(prime[i],counter));

    }
    }
    return ans;

    }

    int vecmogh(vector<pair > one, vector<pair > two)
    {
    if (one.size()!=two.size())
    return 0 ;

    for (int i=0;i<one.size();i++)
    {
    if (one[i].first != two[i].first || one[i].second!=two[i].second)
    return 0 ;
    }

    return 1 ;

    }

    int gcd(int a,int b)
    {
    if (a%b==0)return b;
    a%=b;
    a=gcd(b,a);
    }

  3. Scott Dennison says:

    I’ve tried to convert your pen and paper solution into an algorithm, so that it can work with any maximum bounds of A or B (initially caused by my misreading <=100 and <=1000).

    Please could you look over the code and see if you believe I am on the right track. It SEEMS to work.

    int intPossibilites = (Problem029.MAX_A – 1) * (Problem029.MAX_B – 1);

    int intMaximumPower = 0;
    int intPowerValue = 1;
    while (true) {
    intPowerValue < Problem029.MAX_B) {
    break;
    }
    else {
    intMaximumPower++;
    }
    }

    int[] arrPossibilitiesToRemovePerPower = new int[intMaximumPower+1];
    for (int intNthPower=2; intNthPower<=intMaximumPower; intNthPower++) {
    Set setPossibilitesToRemoveForThisPower = new HashSet();
    for (int intSubPower=1; intSubPower<intNthPower; intSubPower++) {
    int intGoUpTo = (Problem029.MAX_B * intSubPower) / intNthPower;
    int intIncrease = intSubPower / BasicMathUtilities.gcd(intSubPower, intNthPower);
    int intValue = intSubPower;
    if (intSubPower == 1) {
    intValue = 2;
    }
    while (intValue <= intGoUpTo) {
    setPossibilitesToRemoveForThisPower.add(intValue);
    intValue += intIncrease;
    }
    }
    arrPossibilitiesToRemovePerPower[intNthPower] = setPossibilitesToRemoveForThisPower.size();
    }

    for (int intA=2; intA=2; intPower–) {
    double dblApproximateNRoot = Math.pow(intA,(1d / intPower)); // Don’t trust the doubles!.
    int intExactNRoot = -1;
    if (Math.pow(Math.floor(dblApproximateNRoot),intPower) == intA) {
    intExactNRoot = (int)Math.floor(dblApproximateNRoot);
    }
    else if (Math.pow(Math.ceil(dblApproximateNRoot),intPower) == intA) {
    intExactNRoot = (int)Math.ceil(dblApproximateNRoot);
    }
    if (intExactNRoot != -1) {
    int intExactNRootSquareRoot = (int)Math.sqrt(intExactNRoot);
    if ((intExactNRootSquareRoot*intExactNRootSquareRoot) != intExactNRoot) {
    intPossibilites -= arrPossibilitiesToRemovePerPower[intPower];
    break;
    }
    }
    }
    }
    return intPossibilites;

  4. Jean-Marie Hachey says:

    Project Euler – Problem 29

    Summary of tests on the algorithms presented here.

    1) Kristian’s algo
    (26 February 2011)
    http://www.mathblog.dk/files/euler/Problem29.cs
    Number of errors detected : 0

    2) Leo’s algo
    (January 7, 2014 at 16:14)
    Number of errors detected : 77

    3) Scott Dennison’s algo
    (June 24, 2014 at 13:58)
    Number of errors detected : 62

    ___

    Source :
    1) Microsoft Visual C# 2010 Express

  5. Scott Dennison says:

    Just a note, the code I posted was Java not C#, and also relies on a class named BasicMathUtilities with a method named, which can be substituted with your preferred gcd algorithm. It also uses two variables Problem029.MAX_A and Problem029.MAX_B which were not provided.

    For helpfulness, I have included these in the below code block:

    package uk.co.scottdennison.java.soft.tasks.projecteuler.solutions;
    
    import java.util.HashSet;
    import java.util.Set;
    
    public class Problem029 {
    	private static final int MAX_A = 100;
    	private static final int MAX_B = 100;
    
    	public static int gcd(int intLeft, int intRight) {
    		int intTemp;
    		intLeft = Math.abs(intLeft);
    		intRight = Math.abs(intRight);
    
    		while (true) {
    			if (intLeft == 0) {
    				return intRight;
    			}
    			else if (intRight == 0) {
    				return intLeft;
    			}
    			else {
    				intTemp = intLeft % intRight;
    				intLeft = intRight;
    				intRight = intTemp;
    			}
    		}
    	}
    
        // Losely based on the pen-and-paper theory at 'http://www.mathblog.dk/project-euler-29-distinct-terms-sequence-ab/', combined with a lot of my own research into the relationships between nth powers which I had failed to realise was such.
    	public int solve() {
    		int intPossibilites = (Problem029.MAX_A - 1) * (Problem029.MAX_B - 1);
    
    		int intMaximumPower = 0;
    		int intPowerValue = 1;
    		while (true) {
    			intPowerValue &lt;&lt;= 1;
    			if (intPowerValue &gt; Problem029.MAX_B) {
    				break;
    			}
    			else {
    				intMaximumPower++;
    			}
    		}
    
            int[] arrPossibilitiesToRemovePerPower = new int[intMaximumPower+1];
    		for (int intNthPower=2; intNthPower&lt;=intMaximumPower; intNthPower++) {
    			Set&lt;Integer&gt; setPossibilitesToRemoveForThisPower = new HashSet&lt;&gt;();
    			for (int intSubPower=1; intSubPower&lt;intNthPower; intSubPower++) {
    				int intGoUpTo = (Problem029.MAX_B * intSubPower) / intNthPower;
                    int intIncrease = intSubPower / Problem029.gcd(intSubPower, intNthPower);
                    int intValue = intSubPower;
    				if (intSubPower == 1) {
    					intValue = 2;
    				}
    				while (intValue &lt;= intGoUpTo) {
    					setPossibilitesToRemoveForThisPower.add(intValue);
    					intValue += intIncrease;
    				}
    			}
    			arrPossibilitiesToRemovePerPower[intNthPower] = setPossibilitesToRemoveForThisPower.size();
    		}
    
    
    		for (int intA=2; intA&lt;=Problem029.MAX_A; intA++) {
    			for (int intPower=intMaximumPower; intPower&gt;=2; intPower--) {
    				double dblApproximateNRoot = Math.pow(intA,(1d / intPower)); // Don't trust the doubles!.
    				int intExactNRoot = -1;
    				if (Math.pow(Math.floor(dblApproximateNRoot),intPower) == intA) {
    					intExactNRoot = (int)Math.floor(dblApproximateNRoot);
    				}
    				else if (Math.pow(Math.ceil(dblApproximateNRoot),intPower) == intA) {
    					intExactNRoot = (int)Math.ceil(dblApproximateNRoot);
    				}
    				if (intExactNRoot != -1) {
    					int intExactNRootSquareRoot = (int)Math.sqrt(intExactNRoot);
    					if ((intExactNRootSquareRoot*intExactNRootSquareRoot) != intExactNRoot) {
                            intPossibilites -= arrPossibilitiesToRemovePerPower[intPower];
                            break;
    					}
    				}
    			}
    		}
    		return intPossibilites;
    	}
    }
    

    @Jean-Marie Hachey
    By errors, what exactly do you mean. Compilation errors? Incorrect results (and if so, with what dataset?)? Something else.

  6. Jean-Marie Hachey says:

    @ Scott Dennison
    Thank you for your code block in Java.
    Sorry but as indicated the test was done with Microsoft Visual C# 2010 Express.
    (So it doesnt apply to Java)

  7. Jean-Marie Hachey says:

    Table 1
    Distinct terms in the sequence generated by a^b
    for 2 <= a <= 200 and 2 <= b <= 200.

    http://img4.hostingpics.net/pics/514596pe29tab1dpwr.jpg
    ___

    Sources:
    1) Project Euler – Problem 29
    Kristian's algorithm
    http://www.mathblog.dk/files/euler/Problem29.cs
    2) Microsoft Visual C# 2010 Express

  8. Scott Dennison says:

    @Jean-Marie Hachey
    Here is an extended table using my algorithm.

    http://oi57.tinypic.com/o7n9dk.jpg

    Note, I get a different result for 200, which was worrying until I investigated further.
    Using double, which Kristian’s algorithm does (Kristian has rightly admitted hating using doubles), the imprecision causes incorrect results when max_a=max_b=200

    This java code:

            SortedSet&lt;BigInteger&gt; setResultsBigInteger = new TreeSet&lt;&gt;();
            SortedSet&lt;Double&gt; setResultsDouble = new TreeSet&lt;&gt;();
            for (int intA=2; intA&lt;=200; intA++) {
                for (int intB=2; intB&lt;=200; intB++) {
                    BigInteger biResult = BigInteger.valueOf(intA).pow(intB);
                    Double dblResult = Math.pow(intA,intB);
                    setResultsBigInteger.add(biResult);
                    setResultsDouble.add(dblResult);
                }
            }
            System.out.println(setResultsBigInteger.size());
            System.out.println(setResultsDouble.size());
    

    Gives exactly the same result as your table for double, but when using biginteger, which is guaranteed to be precise for this problem, you get over 7000 more results.
    It can be shown to be even worse when max_a=max_b=1000, where the double solution has 118168 results, but the biginteger solution has over 8 times more – 977358.

  9. Jean-Marie Hachey says:

    @Scott Dennison
    Thank you for your Java code and the extended table.
    I have tried an algo in Python written by Lucas Willems:
    http://blog.lucaswillems.com/1216/project-euler-probleme-29

    All the results you got with Java are obtained with this Python program with the limit at 2001.
    (No memory for new parser at 2001).
    The biggest result from this Python program at 1001: 977358.

  10. sinan says:
                Stopwatch sw = new Stopwatch();
                sw.Restart();
                int[] powNum = new int[bound+1]; /*if a number is not a power of a smaller number, it will add 99 to the result
                                                  * otherwise, further calculation will be necessary*/
                int result = 0;
    
                for (int i = 2; i &lt;= bound / 2; i++) {
                    if(powNum[i] != 0)/*i want to represent numbers in the smallest base possible
                                       * i do not want 4^2 to overwrite 2^4 for the number 16 for example*/
                        continue;
                    int current = i * i;
                    int pow = 2;
                    while (current &lt;= bound) {
                        powNum[current] = pow;
                        current *= i;
                        pow++;
                    }
                }
    
                for (int i = 2; i &lt;= bound; i++) {
                    int[] lookupValues = new int[10];/*what is the maximum exponent a number can have in this problem? Math.log(2,bound) could be an upperbound 
                                                      *for this bound value, this is diminishing return if there were many repeating exponent values,
                                                      *it would have been smarter to record them*/
                    if (powNum[i] == 0) { //it is a good idea to check for this since most of the numbers are not a power of anything
                        result += bound - 1;
                    }
                    else if (powNum[i] == 2) { //this also speeds up the code and it is very easy to see that the result will be incremented by 50 (for bound = 100)
                        result += (int)(bound / 2 + 0.99);
                    }
                    else {
                        if (lookupValues[powNum[i]] == 0) {/*for number 8=2^3, the value that needs to be incremented for exponent 3 will not be determined by now
                                                            * so I will count it*/
                            bool[] alreadyOccurred = new bool[bound + 1];
                            for (int j = 1; j &lt; powNum[i]; j++) {
                                for (int k = 1; k &lt;= bound; k++) {
                                    if ((j * k) % powNum[i] == 0) {
                                        alreadyOccurred[j * k / powNum[i]] = true;
                                    }
                                }
                            }
                            for (int j = 1; j &lt;= bound; j++) {
                                if (!alreadyOccurred[j])
                                    lookupValues[powNum[i]]++;
                            }
    
                            result += lookupValues[powNum[i]];
                        }
                        else {/*for 27=3^3, the value for exponent 3 will aready have been recorded from 2^3=8
                               * but as I have mentioned, for such a low frequency of repeated values, this actually slows down the code*/
                            result += lookupValues[powNum[i]];
                        }
                    }
                }
    
                sw.Stop();
                Console.Write(&quot;time:\t{0}\n\n&quot;, ((double)(sw.Elapsed.TotalMilliseconds)).ToString(&quot;0.000 000 ms&quot;));//0.017ms=17us on my machine
  11. Jonathan de La Marche says:

    If you factorize the base number a, we can easily see why 16^51 = 8^68 or (2^4)^51 = 2^204 = (2^3)^68. I wrote a little program in Python which uses this method and keeps the factorized numbers as a string in a set.

    from collections import Counter
    
    def prime_factors(n):
    	factors = []
    	d = 2
    	while n &gt; 1:
    		while n % d == 0:
    			factors.append(d)
    			n /= d
    		d = d + 1
    		if d*d &gt; n:
    			if n &gt; 1: factors.append(n)
    			break
    	return factors
    
    def CountPower(c,b):
    	s = []
    	for k in c:
    		s += [str(k) + '^' + str(b*c[k])]
    		S = '*'.join(s)
    	return S
    
    N = 1000
    A = set()
    for a in range(2,N+1):
    	f = prime_factors(a)
    	c = Counter(f)
    	for b in range(2,N+1):
    		A.add(CountPower(c,b))
    
    print len(A)
  12. Chas Busenburg says:

    Why does the set have to be a SortedSet?

    If you used a HashSet ( if c# has those, running on java knowledge), wouldn’t the Big Oh of the insertion be better, and retrieval is about the same, but you don’t need to retrieve, just the size of the set afterwards. Order doesn’t really matter for this problem, or at least from what i’m seeing.

  13. Wei says:

    2*(49+16) = 116 is wrong for 16 and 81

  14. Krasney says:

    Wei is correct, 2* (49 + 9) is the result for dups (2^3)^(4n/3)

    51, 54, 57, 60, 63, 66, 69, 72, 75

Trackbacks For This Post

  1. Project Euler 29 – Distinct Powers | DL-UAT

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

This site uses cookies.