Volume Number: 16 (2000)
Issue Number: 3
Column Tag: Programmer's Challenge

Programmer's Challenge

by Bob Boonstra, Westford, MA

Rub*k Rotation

Readers sometimes write me to ask whether they can solve the Challenge on a platform other than a Mac. Now, I'm not sure why a reader of this publication would want to do that, but in many cases the Challenge is platform agnostic and quite amenable to solution an another machine. Sometimes I get mail thanking me for creating non-Mac-specific problems. Then again, I also get mail asking why problems that have nothing to do with the Mac appear in a Macintosh publication. Paraphrasing Lincoln, you can't please all of the people all of the time. But it has been a while since we've had a Macintosh-specific problem, so we're going to offer one this month. The last Challenge I competed in required readers to solve the puzzle known as Rub*k's Cube, a 3x3 cube of smaller cubies where each cube face was colored with a different color. Solving that Challenge took me right to the midnight deadline, past many midnights before that one, and actually pushed me over the edge into Challenge retirement (as a contestant, anyway). As I write this column, having watched another midnight come and go, I thought we'd revisit the Cube, but from a completely different perspective. Your Challenge this month is to write a program that will display the cube, animating both rotations of the cube and the moves used to solve the cube. The prototype for the code you should write is:

typedef enum {  /* identify the cube face */
  kFront=0, kBack, kLeft, kRight, kUp, kDown
} CubeFace;

typedef enum {  /* identify the axes of rotation */
  kFrontBack=0, kLeftRight, kUpDown
    /* rotation axis kXY is oriented viewing face X, through the cube, toward face Y */
} CubeAxis;

typedef enum { 
  kClockwise=1, kCounterClockise=-1
} TurnDirection;

void InitCube(
  CWindowPtr cubeWindow,         
    /* window where the rotating cube should be rendered */
  const RGBColor cubeColors[6],   
    /* colors to use in rendering, indexed by CubeFace */
  const short cubieColors[3][3][3], 
    /* cubieColors is the index into cubeColors for individual cubies */
    /* cubeColors[cubieColors[f][r][c]] is the color of the cubie on face f, row r, col c*/
    /* Row and column orientations are as follows: */
    /* cubieColors[kFront][0][0] is adjacent to kUp and kLeft */
    /* cubieColors[kFront][0][2] is adjacent to kUp and kRight */
    /* cubieColors[kBack][0][0] is adjacent to kUp and kRight */
    /* cubieColors[kBack][0][2] is adjacent to kUp and kLeft */
    /* cubieColors[kLeft][0][0] is adjacent to kUp and kBack */
    /* cubieColors[kLeft][0][2] is adjacent to kUp and kFront */
    /* cubieColors[kRight][0][0] is adjacent to kUp and kFront */
    /* cubieColors[kRight][0][2] is adjacent to kUp and kBack */
    /* cubieColors[kUp][0][0] is adjacent to kBack and kLeft */
    /* cubieColors[kUp][0][2] is adjacent to kBack and Right */
    /* cubieColors[kDown][0][0] is adjacent to kFront and kLeft */
    /* cubieColors[kDown][0][2] is adjacent to kFront and kRight */
short cubeWidth,  
    /* size in pixels of a cube in the standard orientation (kFront visible) */
  short stepSize    
    /* granularity to redraw, stepSize steps is one full 360 degree rotation */
);

void QuarterTurn(
  CubeFace face,           /* turn this face one quarter turn */
  TurnDirection direction  /* turn the face in this direction */
    /* turn orientation is looking at the face from outside the cube toward the center 
      of the cube */
);

void RotateCube(
  CubeAxis axis, /* rotate cube about the specified axis */
  TurnDirection direction,  
    /* rotate the cube in this direction about the specified axis */
    /* clockwise about the kFrontBack face is determined viewing from the kFront 
        face, through the cube to the kBack face */
  short stepsToTurn
);

void TermCube(void);

This Challenge will start with a call to your InitCube routine providing a number of problem parameters. The window in which you should render the cube, a color window with a pixelSize of 32 bits, will be provided in cubeWindow. The colors making up the cube faces will be provided in cubeColors, and the individual cubie colors will be specified by cubieColors as indices into cubeColors. You should center your display of the cube in the cubeWindow, and size the cube so that it is cubeWidth pixels on a side, in its initial position, oriented along the u-v axes, and viewed along the cube normal. InitCube should draw the cube in its initial posiiton, viewing the kFront face along the kFrontBack axis, with the kUp face at the top of the cube, perpendicular to the view plane. The cube may be displayed with a perspective projection, from a reasonable distance, or from an orthographic projection from infinity.

The final parameter provided to InitCube is the stepSize, used by the QuarterTurn and RotateCube calls in animating cube turns and rotations. The stepSize parameter specifies how granular the animations and rotations should be, with stepSize rotation steps constituting one full rotation. The value of stepSize will be a multiple of 4, so that quarter turns will contain an integral number of steps.

The QuarterTurn routine is called to turn the 9 cubies on one face of the cube by 90 degrees relative to the rest of the cube. When QuarterTurn is called, you should rotate the specified face in the specified direction, animated in increments determined by stepSize. Any internal portions of the cube not normally visible, but visible during the turn, should be displayed in a gray or black color of your choice.

The RotateCube routine is called to rotate the entire cube, respectively. When RotateCube is called, you should rotate the entire cube about the specified cube axis in the specified direction. As multiple calls are made to RotateCube, the axes of rotation will become arbitrarily oriented with respect to the view vector, although the cube will remain centered in the cubeWindow.

Finally, the TermCube routine will be called at the end of the test, allowing you to clean up any memory you might have allocated before returning. There may be multiple test cases, each bracketed with an InitCube and a TermCube call.

The commentary in the code for the CubeAxis describe the orientation used to perform clockwise and counterclockwise rotations of the cube. Similarly, the commentary for the CubeFace in the QuarterTurn prototype describes the orientation for performing face quarter turns. If these explanations are not clear, please send me a note on the Challenge mailing list for clarification.

Scoring will be based first on the quality of the accuracy of the display, and the absence of tearing or other display anomalies. Among the accurate entries with acceptable display quality, the winner will be the solution requiring the least execution time.

This will be a native PowerPC Challenge, using the CodeWarrior Pro 5 environment. Solutions may be coded in C, C++, or Pascal. Solutions in Java will also be accepted, but Java entries must be accompanied by a test driver that uses the interface provided in the problem statement.

Congratulations to Ernst Munter (Kanata, ON, Canada) for taking first place in the March Sum Of Powers Challenge, narrowly beating out the second place entry from Miklos Fazekas. You might recall that the March Challenge required you to find a set of terms which, added or subtracted together, formed a given positive integer. The terms were required to be integers raised to a power greater than 1. Scoring was based on the number of terms used to form the result and on the amount of execution time used to calculate the solutions, with a penalty of 1 term per 100 milliseconds of execution time.

After the Challenge was published, several people on the Challenge mailing list pointed out that allowing subtraction of powers of 2 made the problem trivial. Using the fact that (n+1)^2 - n^2 = 2n+1, one can form any odd number by subtracting the squares of two sequential integers, and any even number by adding or subtracting 1 to the squares forming an adjacent odd number. So I amended the problem statement to prohibit subtraction of squared terms.

Ernst's entry starts by calculating the value of all integers raised to a power greater than 2 that fit into a signed long integer. The code offers an option to calculate this array at startup time, but for scoring purposes, this calculation was performed during the first SumOfPowers call. These calculated values are used to to exhaustively look for a two-term solution, and then a three-term solution if no two-term solution is found. Finally, while Ernst conjectures that three terms are the most ever needed, a conjecture that was not disproven by my tests, he provides a recursive solution in the event no three-term answer is found.

Miklos uses some facts from number theory to determine when a number can be written as the sum of two or three squares, and observes that every number can be written as the sum of four squares. He considers only the powers 2, 3, and 5, noting that there are not many numbers raised to the power 7 that fit into a long integer. That assumption was costly, however, as his entry generated an extra term for numbers like 5054. Ernst formed 5054 as 4415^2 - 11^7, while Miklos formed it with an extra term: 71^2 + 3^2 + 2^2.

I tested the entries with 4330 test cases, consisting of 10 sequences of numbers of lengths between 10 and 1000. Ernst and Claes Wihlborg both solved the test cases by generating 11880 total terms, but Ernst's entry took only 10% of the time that Claes' did. Miklos Fazekas' entry actually solved the test cases about 10% faster than Ernst's entry did, but generated enough extra terms to result in a slightly poorer score. The scores were close enough that I considered code size as a tie-breaker, and Ernst's entry was significantly more compact.

The table below lists, for each of the solutions submitted, the total score, based on the number of terms generated and total execution time; the execution time itself; the total terms generated for all test cases; and an error indicator. It also provides the code size, data size, and programming language used by each entry. As usual, the number in parentheses after the entrant's name is the total number of Challenge points earned in all Challenges prior to this one.

Top Contestants

Listed here are the Top Contestants for the Programmer's Challenge, including everyone who has accumulated 10 or more points during the past two years. The numbers below include points awarded over the 24 most recent contests, including points earned by this month's entrants.

Here is Ernst's winning Sum Of Powers solution:

SumOfPowers.cp
Copyright © 2000
Ernst Munter

/*

Problem
———-
Given a number (the "result") between 0 and 2^31, find a smallest set of terms which sum to the number.  A term
is sign * x ^ y where y must be 2 or more and sign may only be negative if y >= 3.

Constraint
—————
All intermediate results, as terms are added or subtracted, must not overflow a 32-bit signed integer variable.

Solution
————
Basically by systematic trial and error.  My conjecture is that all results can be made with no more than 3
terms, and a small number of functions are called in succession to try 1-term, 2-term and 3-term sets. A fully
exhaustive test of all 2^31 cases was not done, so my conjecture may be false.  A final recursive function is
provided to find sets larger than 3 terms, just in case.

Assumption
—————
SumOfPowers() will not be called with maxTerms < 3.  This assumption may be overridden by setting
CHECKMAXTERMS_3 = 1.

Arrays of powers of x are precalculated at startup time (takes about 20 millisecs).  If this is considered to be
an unfair precalculation, the macro AUTOINIT should be set to 0, and precalculation will be done as part of the
first call to SumOfPowers().
[ The solution was tested with AUTOINIT set to zero - ChallengeMeister Bob. ]

*/

#include <stdio.h>
#include <stdlib.h>
#define NDEBUG
#include <assert.h>
#include "SumOfPowers.h"

// Turn AUTOINIT on to save about 20msec in the first call.
#define AUTOINIT 0

// Turn CHECKMAXTERMS_3 on to allow tests with maxTerms < 3.
#define CHECKMAXTERMS_3 0

// Count number of powers of N which fit in a 31-bit long.
enum {
	k2 = 46340,		// = sqrt(2^31 - 1) 
	kHigh = 1698,	// = sum of ((2^31 - 1) ^ (1/n)), n = 3 .. 31
		
	kAll=k2+kHigh - 1	
	// = 48037; actually slightly fewer locations will be used 
	// because of entries like 2^6 == 4^3 only one will be stored.
};

struct Entry {
// Similar to struct IntegerPower but tracks v = base ^ power.
	unsigned long	v;
	unsigned short	base;	
	unsigned short	power; 
	Entry(){}
	Entry(unsigned long vx,int b,int p):v(vx),base(b),power(p){}
};

static int Cmp(const void* a,const void* b)
{
// Needed for quicksort
	Entry* pa=(Entry*)a;
	Entry* pb=(Entry*)b;
	return pa->v - pb->v;
// simple delta works here because all values are positive longs
}

Search
inline unsigned long* Search(long key,unsigned long e[],
			int length)
{
// Binary search, returns the array element which is <= key.
	unsigned long l=0,r=length-1,m,v;
	do {
		m=(l+r)>>1;		
		v=e[m];
		if (key==v) return e+m;
		else if (key < v) r=m-1;
		else l=m+1;
	} while(l<=r);
	if (key < v)
		return e+r;
	else
		return e+m;
}

struct PowerSolve
static struct PowerSolve
{
// Collection of data and functions to solve SumOfPowers
	int numAll;			// actual size of array of all entries
	int numHigh;		// actual size of array of high entries
	unsigned long vHigh[kHigh];
	Entry high[kHigh];	// sorted powers x^p with p>=3
	unsigned long vAll[kAll];
	Entry all[kAll];	// sorted powers x^p with p>=2
	
#if AUTOINIT
	PowerSolve(){Create();}
#endif
	void Create();

	bool IsInitialized(){return numAll;}
	
	int Solve(unsigned long result,IntegerPower terms[],
									long maxTerms);
	
	void CopyPositiveTerm(unsigned long* v,IntegerPower & ip)
	{
// Copies an Entry into an IntegerPower
		Entry* p=&all[v-vAll];
		ip.value=p->base;
		ip.power=p->power;
		ip.sign=1;
	}
	
	void CopyNegativeTerm(unsigned long* v,IntegerPower & ip)
	{
		Entry* p=&high[v-vHigh];
		ip.value=p->base;
		ip.power=p->power;
		ip.sign=-1;
	}
	
Find2
	bool Find2(unsigned long V,unsigned long* v0,
							IntegerPower terms[])
	{
// Exhaustive search for a pair of terms to match V.
// Trying X+Y:
		unsigned long* saveV=v0;
		unsigned long* v1=v0=Search(V/2,vAll,1+v0-vAll);		
		if (V == *v0+*v1) 
		{
			CopyPositiveTerm(v0,terms[0]);
			CopyPositiveTerm(v1,terms[1]);				
			return true;
		}
		do {		
			long z,delta=V-*v1;
			while ((z=*++v0) <= delta) {
				if (z == delta) 
				{
					CopyPositiveTerm(v0,terms[0]);
					CopyPositiveTerm(v1,terms[1]);					
					return true;
				}
			}					
						
			if (v1<=vAll) break; 
			delta=V-z;
			while ((v1>vAll) && ((z=*—v1) >= delta)) {	
				if (z == delta) 
				{
					CopyPositiveTerm(v0,terms[0]);
					CopyPositiveTerm(v1,terms[1]);		
					return true;
				}
			}
			
		} while (v0 < saveV);
		
// Trying X-Y:		
		v1=vHigh;
		unsigned long sum=V+*v1;
		v0=Search(sum,vAll,numAll);	
		unsigned long* endAll=vAll+numAll-1;	
		unsigned long* endHigh=vHigh+numHigh-1;
		for (;;) 
		{						
			if (sum==*v0) 
			{
				CopyPositiveTerm(v0,terms[0]);
				CopyNegativeTerm(v1,terms[1]);			
				return true;
			}
			if (v1>=endHigh) break;
				
			sum=V+*++v1;
			while ((v0<endAll) && (v0[1] <= sum)) 
				v0++;
		} 
		return false;
	}
	
FindSum3
	bool FindSum3(long V,unsigned long* v0,IntegerPower terms[])
	{
// Trying X + Y + Z	
		unsigned long* v2=vAll;
		unsigned long* v10=Search(V-*v0,vAll,numAll);
		for (;v0>=vAll;v0—)
		{
			long diff=V - *v0;
			unsigned long* v1=Search(diff,vAll,numAll);
			v2=vAll;
			do {
				diff=V - *v0 - *v1;
				if (diff<0) break;
				
				long z;
				while ((v2<=v1) && ((z=*v2) < diff))
					v2++;
			
				if (z == diff)
				{
					CopyPositiveTerm(v0,terms[0]);
					CopyPositiveTerm(v1,terms[1]);
					CopyPositiveTerm(v2,terms[2]);		
					return true;
				}
			} while (—v1 >= v2);
		}
		return false;
	}
		
FindDelta3
	bool FindDelta3(long V,IntegerPower terms[])
	{
// Trying X - Y + Z	
		for (unsigned long* v1=vHigh;v1<vHigh+numHigh;v1++) {	
			unsigned long sum=V+*v1;
			unsigned long* v0=Search(sum/2,vAll,numAll);
			long diff=sum-*v0;
			unsigned long* v2=Search(diff,vAll,numAll);
			
			if (*v2 == diff)
			{
				CopyPositiveTerm(v0,terms[0]);
				CopyNegativeTerm(v1,terms[1]);
				CopyPositiveTerm(v2,terms[2]);			
				return true;
			}
			for (;;)
			{
				v0++;
				if (v0>=vAll+numAll) break;
				diff=sum-*v0;
				long z;
				while ((v2>=vAll) && ((z=*v2) > diff)) 
					v2—;
				
				if (z == diff)
				{
					CopyPositiveTerm(v0,terms[0]);
					CopyNegativeTerm(v1,terms[1]);
					CopyPositiveTerm(v2,terms[2]);		
					return true;
				}
			}
		} 
		return false;
	}
	
FindFinal
	int FindFinal(unsigned long V,unsigned long* v0,
		IntegerPower terms[],long maxTerms)
	{
// Subtracts the largest available term and solves for the remainder.
// Always finds a (perhaps suboptimal) set of terms,
//		unless maxTerms is too small.
		CopyPositiveTerm(v0,terms[0]);
		int subTerms=Solve(V-*v0,terms+1,maxTerms-1);
		
		if (subTerms) return subTerms+1;
		return 0;
	}

} PS;

Compact
inline int Compact(Entry e[],const int numE)
{
// Compresses Entry array e[] by removing duplicate entries (same v)
	int numOut=1;
	Entry* input=e;
	Entry* out=e;
	for (int i=1;i<numE;i++)
	{
		unsigned long v=(++input)->v;
		if (v != out->v)
			*++out=*input;
	}
	return 1+out-e;
} 

Create
// The following two functions are not inlined.  This is to control the
// compiler's (CodeWarrior Pro 5.3) inlining behaviour of the rest.			
void PowerSolve::Create()
{
// Constructs all[] and high[] arrays by exhaustive enumeration
	int nAll=0;
	int nHigh=0;

// high[] array starts with 1^3 
// The high[] array is constructed by multiplying with base 
	high[nHigh++]=Entry(1,1,3);		
	long base;
	for (base=2;base<=kHigh;base++)
	{
		unsigned long v=base*base;
		unsigned long maxV=0x7FFFFFFF/base;
		int power=2;
		for (;v<=maxV;) {						
			v *= base;
			assert(nHigh<kHigh);
			high[nHigh++]=Entry(v,base,++power);							
		}
	}
		
// Sort and compress the high[] array:
	qsort(high,nHigh,sizeof(Entry),Cmp);
	nHigh=Compact(high,nHigh);
	
		
// The all[] array is constructed by additive accumulation of the
// square terms and merging with the terms from the sorted high[]
	unsigned long value=4;
	unsigned long delta=5;
	unsigned long lastV=0;
	Entry* hp=high;
	for (base=2;base<=k2;base++)
	{
		while ((hp<(high+nHigh)) && (value >= hp->v))
		{
			if (hp->v != lastV)
			{
				assert(nAll<kAll);
				vAll[nAll]=lastV=hp->v;
				all[nAll++]=*hp;	
			}
			hp++;
		}
		if (value != lastV)
		{
			assert(nAll<kAll);
			vAll[nAll]=lastV=value;
			all[nAll++]=Entry(value,base,2); 
		}
		value+=delta;
		delta+=2;			
	}
	for (int i=0;i<nHigh;i++)
		vHigh[i]=high[i].v;
		
	numHigh=nHigh;
	numAll=nAll;
}

Solve
int PowerSolve::
Solve(unsigned long result,IntegerPower terms[],long maxTerms)
{
// Solve() function solves for result by finding the fewest possible
// terms, up to maxTerms.  

	unsigned long* v0=Search(result,vAll,numAll);
// A 1-term result is found directly in the all[] array
	if (*v0 == result)
	{
		CopyPositiveTerm(v0,terms[0]);
		return 1;
	}
	
#if CHECKMAXTERMS_3
	if (maxTerms < 2) return 0;
#endif
		
// A 2-term array may be a sum or difference of terms		
	if (Find2(result,v0,terms)) return 2;		
				
#if CHECKMAXTERMS_3
	if (maxTerms < 3) return 0;
#endif
		
// 3-term results may be:
//		all[x] + all[y] + all[z]		
	if (FindSum3(result,v0,terms)) return 3;
			
// 	or	all[x] - high[y] + all[z]
//  (very few results require this form)			
	if (FindDelta3(result,terms)) return 3;
			
//  or	all[x] - high[y] - high[z]
//  but this variation does not seem to be needed				
				
//  Lengthy but not exhaustive tests have not turned up a single
//	instance of result that could not be solved with three or fewer terms.  
//  But just in case there is a term that needs more terms,				
//  FindFinal() will provide a general solution.
	return FindFinal(result,v0,terms,maxTerms);
}

SumOfPowers
long /* number of factors */ SumOfPowers (
	long result,			/* terms need to sum to this result */
	IntegerPower terms[],	/* return terms sign*integer^power here */
	long maxTerms			/* maximum number of terms allowed */
) {	
		
#if CHECKMAXTERMS_3
	if (maxTerms<=0) return 0;
#endif
		
	if (!PS.IsInitialized()) PS.Create();
	
	return PS.Solve(result,terms,maxTerms);
}