/******************************************************************************/
/* */
/* FIND SOLUTIONS OF (a^3+b^3)/(a+b)=T^3 */
/* 03/21/16 (dkc) */
/* */
/******************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0a.h"
#include "table2.h"
//
// "weak" Furtwangler conditions
//
unsigned int test(unsigned int d) {
unsigned int flag,quotient,f,remd,p,ps;
p=3;
ps=9;
flag=0;
if ((d/2)*2!=d) {
quotient=d/p;
if (quotient*p==d)
f=quotient;
else
f=d;
quotient=f/ps;
remd=f-quotient*ps;
if ((remd!=1)&&(remd!=8))
flag=1;
}
else {
quotient=d/p;
if (quotient*p!=d) {
f=d/2;
quotient=f/ps;
remd=f-quotient*ps;
if ((remd!=1)&&(remd!=8))
flag=1;
}
}
return flag;
}
//
// search for prime factors
//
unsigned int search(unsigned long long s, unsigned int k,
unsigned int *table3) {
unsigned int i,j,l,m,count;
unsigned long long v;
count=0;
for (i=0; i<=k; i++) {
m=0;
l=table3[i];
v=s/(unsigned long long)l;
if (s!=(v*(unsigned long long)l))
continue;
aloop:
m=m+1;
s=v;
v=s/(unsigned long long)l;
if (s==(v*(unsigned long long)l))
goto aloop;
j=m/3;
j=j+j+j;
if (j!=m)
return(0); // not a cube
else
count=count+1;
if (s==1)
return(count);
}
if (s!=1)
return(0); // not completely factored
else
return(count);
}
//
// inner loop
//
unsigned int search(unsigned long long s, unsigned int k,
unsigned int *table3);
unsigned int eloop(unsigned int d, unsigned int k, unsigned int n,
unsigned int outsize, unsigned int *output,
unsigned int *s, unsigned int size,
unsigned int *ptable, unsigned char *flag,
unsigned int *error, unsigned int *table3) {
unsigned int p,f,g,count,eend,fstart,index,rem,temp,et;
int ebeg;
unsigned long long a,a2,c,e,t,u;
ebeg=d-1;
if ((unsigned int)ebeg>size)
eend=ebeg-size+1;
else
eend=1;
fstart=d/size;
if (d!=fstart*size)
fstart=fstart+1;
for (f=fstart; f>0; f--) {
for (g=0; g<size; g++)
flag[g]=0;
for (g=0; g<size; g++) {
p=ptable[g];
if (p>d)
break;
if (d!=(d/p)*p)
continue;
rem=(d-eend)-((d-eend)/p)*p;
while ((eend+rem)<=(unsigned int)ebeg) {
flag[rem]=1;
rem=rem+p;
}
}
index=0;
for (g=0; g<(ebeg-eend+1); g++) {
if (flag[g]==0) {
s[index]=eend+g;
index=index+1;
}
}
for (g=0; g<index; g++) {
et=s[g];
if ((d+et)==((d+et)/3)*3)
flag[g]=1;
else
flag[g]=0;
}
a=(unsigned long long)d;
a2=a*a;
for (g=0; g<index; g++) {
e=s[g];
c=a2-a*e+e*e;
if (flag[g]==1)
c=c/3;
u=(unsigned long long)exp(log((double)c)/3.0);
t=u*u*u;
temp=0;
if (t==c)
temp=1;
u=u+1;
t=u*u*u;
if (t==c)
temp=1;
if (temp!=0) {
count=search(c,k,table3);
if (count==0) {
error[0]=99;
return(n);
}
output[n]=d;
output[n+1]=(unsigned int)e;
output[n+2]=count;
n=n+3;
if (n>outsize)
return(n);
}
}
ebeg=ebeg-size;
if (ebeg<=0)
break;
if ((unsigned int)ebeg>size)
eend=ebeg-size+1;
else
eend=1;
}
return(n);
}
//
// outer loop
//
#include <stdio.h>
unsigned int test(unsigned int d);
unsigned int eloop(unsigned int d, unsigned int k, unsigned int n,
unsigned int outsize, unsigned int *output, unsigned int *s,
unsigned int ssize, unsigned int *ptable, unsigned char *flag,
unsigned int *error, unsigned int *table3);
unsigned int dloop(unsigned int din, unsigned int k, unsigned int n,
unsigned int outsize, unsigned int *output, unsigned int
*s, unsigned int ssize, unsigned int *ptable, unsigned char *flag,
unsigned int dend, unsigned int *error,
unsigned int *table3) {
unsigned int i,oldn,d;
oldn=0;
for (i=din; i>=dend; i--) {
d=i;
if (test(d)!=0)
continue;
n=eloop(i,k,n,outsize,output,s,ssize,ptable,flag,error,table3);
if (error[0]==99)
break;
if (n!=oldn)
printf(" d=%d, n=%d \n",i,n/3);
oldn=n;
if (n>outsize) {
error[0]=6;
break;
}
}
return(n);
}
/*****************************************************************************/
/* */
/* FACTOR (a**p + b**p)/(a + b) */
/* 03/21/16 (dkc) */
/* */
/* This C program finds a and b such that (a**p + b**p)/(a + b) is a cube */
/* or p times a cube. The "weak" Furtwangler conditions are assumed to be */
/* applicable (to reduce execution time). p is set to 3. */
/* */
/*****************************************************************************/
extern char *malloc();
unsigned int dloop(unsigned int a, unsigned int b, unsigned int c,
unsigned int d, unsigned int *e, unsigned int *s,
unsigned int ssize, unsigned int *ptable, unsigned char *flag,
unsigned int dend, unsigned int *error,
unsigned int *table3);
int main ()
{
unsigned int dbeg=100000; // starting "a" value
unsigned int dend=3; // must be greater than 2
unsigned int output[20001*3];
unsigned int error[10];
unsigned int d,e;
unsigned int tsize=17984;
unsigned int t3size=4784;
unsigned int ssize=10000;
unsigned int outsiz=20000*3;
unsigned int n=0;
unsigned int i;
unsigned int *sieve;
unsigned char *flag;
FILE *Outfp;
Outfp = fopen("output.dat","w");
i=(unsigned int)sqrt((double)dbeg);
if (i>table2[tsize-1]) {
printf("ssize too small: %d %d \n",i,table2[tsize-1]);
return(0);
}
sieve=(unsigned int*) malloc((ssize+1)*4);
if (sieve==NULL) {
printf("not enough memory \n");
return(0);
}
flag=(unsigned char*) malloc(ssize+1);
if (flag==NULL) {
printf("not enough memory \n");
return(0);
}
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
d=0;
e=0;
error[0]=0; // clear error array
n=dloop(dbeg,t3size-1,n,outsiz,output,sieve,ssize,table2,flag,dend,error,table3);
output[n]=0xffffffff;
printf("error=%d \n",error[0]);
fprintf(Outfp," count=%d \n",(n+1)/3);
for (i=0; i<(n+1)/3; i++)
fprintf(Outfp," %#10x, %#10x, %d, \n",output[3*i],output[3*i+1],
output[3*i+2]);
fclose(Outfp);
return(0);
}