Sam & Polly

I bumped into this over at XKCD today... being bored, I wrote a bit of software to solve it for me (once I had figured out the method). I then had a look at the discussion page (to check my solution) and discovered that there was a conjecture about it working without the x, y <= 97 limit. Since I already had the code, I modded it to work up to an arbitrary values (well... beyond about 46000 it is likely to fall over due to the fact that it works with signed 32-bit numbers). Putting in 1000 for the limit actually gives a second solution: 16, 73. Code below (if you actually use this, please remember to turn on compiler optimisations... it will cut computation time by about 50%). Note that there are two extremely inefficient sections to this code, so up to about N=2000 it works fine but beyond that execution will take more than just a few seconds - I don't know if the algorithm is exponential complexity or just a high order polynomial (and I don't care)... but be warned.

#include <iostream>
using namespace std;


#define N 97
#define NPRIMES (N)
#define NCANDS (N)
#define NFACTS (N)
#define NPC 1000000


struct sPC {
int product, firstsum;
bool unique;
};


int primes[NPRIMES];
int candidates[NCANDS];
int factors[NFACTS];
int factorset[NFACTS];
sPC productcandidates[NPC];


void calculatesolution(sPC *ppc) {
cout << "Solution has a sum of " << ppc->firstsum << " and a product of " << ppc->product << "\n";
cout << "The solution is therefore: " << (ppc->firstsum - ((int) sqrt((double) ((ppc->firstsum * ppc->firstsum)-(4 * ppc->product)))))/2 << ", " << (ppc->firstsum + ((int) sqrt((double) ((ppc->firstsum * ppc->firstsum)-(4 * ppc->product)))))/2 << "\n\n";
}


int main() {
int h, i, j, k; // General iterators
int x, y; // values as per question
int candidatecount, sum, product, factorcount, hfactorcount;
int minx, maxx, factorsetcount, pccount;
bool allindeterminate;


primes[0] = 2;
for(i=3,k=1;i<=N;i++) {
for(j=0;(j<k)&&(i%primes[j]!=0);j++);
if(j==k) primes[k++]=i;
}


// for all possible sums
for(candidatecount=0,sum=6;sum<=(N*2);sum++) {
minx = (sum < (3 + N)) ? 3 : sum - N;
maxx = sum / 2;


// for each possible pair of values for this sum
for(allindeterminate=true,x=minx;x<=maxx;x++) {
y = sum - x;
product = x * y;


// generate a list of all of the prime factors of the product, p
for(factorcount=0, k=0; product!=1; ) {
if((product % primes[k])  == 0) {
factors[factorcount++]=primes[k];
product /= primes[k];
}
else k++;
}
product = x * y;


// work out if this list can generate more than one in-range product
// for every number of primes to have in the product
for(hfactorcount=factorcount/2,factorsetcount=1,i=0;(factorsetcount <= hfactorcount) && (i < 2); factorsetcount++) {


// init a list of indexes to the products
for(j=0;j<factorsetcount;j++) factorset[j]=j;


// check this initial permutation is/isn't a valid product
for(j=0,h=1;j<factorsetcount;j++) h *= factors[factorset[j]];
j = product / h;
if((h>=3)&&(h<=N)&&(j>=3)&&(j<=N)) {
if(i==0) {
i++;
k = h;
}
else if((k!=h)&&(k!=j)) i++;
}


// iterate through all permutations of the list
while((i<2)&&(factorset[0]!=factorcount-factorsetcount)) {
// rearrange index list
for(j=factorsetcount-1;(j>=0)&&(factorset[j]==(factorcount+j-factorsetcount));j--); // find the lattermost "non-stuck" element
for(factorset[j++]++;j<factorsetcount;j++) factorset[j]=factorset[j-1]+1; // move it along one and move all other items to be immediately after it


// check this permutation
for(j=0,h=1;j<factorsetcount;j++) h *= factors[factorset[j]];
j = product / h;
if((h>=3)&&(h<=N)&&(j>=3)&&(j<=N)) {
if(i==0) {
i++;
k = h;
}
else if((k!=h)&&(k!=j)) i++;
}
}
}
if(i!=2) allindeterminate = false;
}


// Check if it is a candidate
if(allindeterminate && (minx != maxx)) candidates[candidatecount++]=sum;
}


// Find the list of all products that can be produced by the candidate sums
// also mark if they are unique or not
for(i=0,pccount=0;i<candidatecount;i++) {
minx = (candidates[i] < (N+3)) ? 3 : candidates[i] - N;
maxx = candidates[i] / 2;
for(x=minx;x<=maxx;x++) {
y = candidates[i] - x;
product = x * y;
for(j=0;(j<pccount)&&(productcandidates[j].product!=product);j++);
if(j==pccount) {
productcandidates[pccount].firstsum = candidates[i];
productcandidates[pccount].product = product;
productcandidates[pccount].unique = true;
pccount++;
}
else productcandidates[j].unique = false;
}
}


// delete all non-unique products
for(i=0,j=0;i<pccount;i++) {
if(productcandidates[i].unique) {
productcandidates[j].product  = productcandidates[i].product;
productcandidates[j].firstsum = productcandidates[i].firstsum;
productcandidates[j].unique   = productcandidates[i].unique;
j++;
}
}
pccount = j;


// find possible solutions by looking for candidate sums with just one unique prouct


if(productcandidates[0].firstsum != productcandidates[1].firstsum) calculatesolution(&productcandidates[0]);
if(productcandidates[pccount-1].firstsum != productcandidates[pccount-2].firstsum) calculatesolution(&productcandidates[pccount-1]);
for(i=1;i<pccount-1;i++) if((productcandidates[i-1].firstsum != productcandidates[i].firstsum)&&(productcandidates[i].firstsum != productcandidates[i+1].firstsum)) calculatesolution(&productcandidates[i]);


return 0;
}