/*****************************************************************************/
/* */
/* FACTOR (a**p+b**p)/(a+b) */
/* 11/03/06 (dkc) */
/* */
/* This C program finds a and b such that (a**p + b**p)/(a + b) is a cube */
/* or p times a cube. p is set to 3. Only a and b where (a**p+b**p)/(a+b) */
/* has one distinct prime factor are output. */
/* */
/* The output is "a, b". If 2p divides a and p**2 does not divide a, then */
/* an error is indicated ("error[1]" is set to a non-zero value). b and */
/* a-b are treated similarly. If 2p divides a+b and p**3 does not divide */
/* a+b, then an error is indicated. If 2 divides a, p does not divide a, */
/* and a/2 is not a pth power modulus p**2, then an error is indicated */
/* ("error[2]" is set to a non-zero value). If a+b is odd, p divides a+b, */
/* and (a+b)/p is not a pth power modulus p**2, then an error is indicated */
/* ("error[3]" is set to a non-zero value). If a+b is odd, p divides a+b, */
/* and p**2 divides a+b, then an error is indicated. If a is odd and p */
/* divides a, then an error is indicated ("error4" is set to a non-zero */
/* value). b and a-b are treated similarly. */
/* */
/* The rest of Proposition (5) can be verified by examining the output. */
/* "prop5c.c" factors a, b, a-b, and a+b (the output of this program is */
/* input to "prop5c.c" [when "sumdif" is 1]). */
/* */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0a.h"
void dummy(unsigned int a, unsigned int b, unsigned int c);
void sum(unsigned int *addend, unsigned int *augend);
void differ(unsigned int *minuend, unsigned int *subtrahend);
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void quotient(unsigned int *a, unsigned int *b, unsigned int);
unsigned int lmbd(unsigned int mode, unsigned int a);
int main ()
{
unsigned int p=3; // input prime
unsigned int dbeg=5000; // starting "a" value
unsigned int dend=1; // ending "a" value
//unsigned int stop=0x2e831;
unsigned int sumdif=1; // select [(a**p+b**p)/(a+b)] if "sumdif" is non-zero
// or [(a**p-b**p)/(a-b)] otherwise
extern unsigned short table3[];
extern unsigned int output[];
extern unsigned int error[];
unsigned int t3size=2556;
unsigned int outsiz=2000;
unsigned int n=0;
unsigned int d,e,a,b,temp;
unsigned int i,j,k,l,lp,m,ps,pc,flag;
unsigned int S[2],T[2],V[2],X[3];
double croot2,croot4,halfcr4;
FILE *Outfp;
Outfp = fopen("out5a.dat","w");
croot2=1.259921;
croot4=1.587401;
halfcr4=croot4*((double)(0.5));
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
ps=p*p;
pc=ps*p;
for (d=dbeg; d>=dend; d--) {
for (e=d-1; e>0; e--) {
dummy(d,e,2);
// if (e!=stop) continue;
/*******************************/
/* check for common factors */
/*******************************/
if((d==(d/2)*2)&&(e==(e/2)*2)) continue;
if((d==(d/3)*3)&&(e==(e/3)*3)) continue;
if((d==(d/5)*5)&&(e==(e/5)*5)) continue;
if((d==(d/7)*7)&&(e==(e/7)*7)) continue;
/***********************/
/* Euclidean G.C.D. */
/***********************/
a=d;
b=e;
if (b>a) {
temp=a;
a=b;
b=temp;
}
loop: temp = a - (a/b)*b;
a=b;
b=temp;
if (b!=0) goto loop;
if (a!=1) continue;
/************************************/
/* compute (d**p + e**p)/(d + e) */
/************************************/
S[0]=0;
S[1]=d;
for (i=0; i<p-1; i++) {
bigprod(S[0], S[1], d, X);
S[0]=X[1];
S[1]=X[2];
}
T[0]=0;
T[1]=e;
for (i=0; i<p-1; i++) {
bigprod(T[0], T[1], e, X);
T[0]=X[1];
T[1]=X[2];
}
if (sumdif==1) {
sum(S, T);
temp=d+e;
if (((d+e)/p)*p==(d+e))
temp=temp*p;
quotient(T, S, temp);
}
else {
differ(S, T);
temp=d-e;
if (((d-e)/p)*p==(d-e))
temp=temp*p;
quotient(T, S, temp);
}
/************************************************/
/* look for prime factors using look-up table */
/************************************************/
if (S[0]==0)
l = 32 - lmbd(1, S[1]);
else
l = 64 - lmbd(1, S[0]);
j=l-(l/3)*3;
l=l/3;
l = 1 << l;
if (j==0)
lp=(int)(((double)(l))*halfcr4);
if (j==1) {
lp=l;
l=(int)(((double)(l))*croot2);
}
if (j==2){
lp=(int)(((double)(l))*croot2);
l=(int)(((double)(l))*croot4);
}
lp=lp-1;
l=l+1;
if (l>table3[t3size-1]) {
error[0]=5;
goto bskip;
}
else {
j=0;
for (i=0; i<t3size; i++) {
if (table3[i] < lp) j=i;
else break;
}
k=j;
for (i=j; i<t3size; i++) {
if (table3[i] < l) k=i;
else break;
}
}
for (i=j; i<=k; i++) {
m=0;
l = table3[i];
quotient(S, V, l);
bigprod(V[0], V[1], l, X);
if ((S[0]!=X[1]) || (S[1]!=X[2])) continue;
aloop: S[0]=V[0];
S[1]=V[1];
m=m+1;
quotient(S, V, l);
bigprod(V[0], V[1], l, X);
if ((S[0]==X[1]) && (S[1]==X[2])) goto aloop;
if ((m/3)*3!=m)
goto askip;
else
break;
}
if ((m/3)*3!=m) continue;
if ((S[0]!=0) || (S[1]!=1)) continue;
if (n+1>outsiz) {
error[0]=6;
n=n-2;
}
if ((d/2)*2==d) {
if ((d/p)*p==d) {
if ((d/ps)*ps!=d)
error[1]+=1;
}
else {
flag=0;
if ((((d/2)-1)/ps)*ps==((d/2)-1))
flag=1;
if ((((d/2)+1)/ps)*ps==((d/2)+1))
flag=1;
if (flag==0)
error[2]+=1;
}
}
else {
if ((d/p)*p==d)
error[4]+=1;
}
if ((e/2)*2==e) {
if ((e/p)*p==e) {
if ((e/ps)*ps!=e)
error[1]+=1;
}
else {
flag=0;
if ((((e/2)-1)/ps)*ps==((e/2)-1))
flag=1;
if ((((e/2)+1)/ps)*ps==((e/2)+1))
flag=1;
if (flag==0)
error[2]+=1;
}
}
else {
if ((e/p)*p==e)
error[4]+=1;
}
if (sumdif!=0) {
if (((d-e)/2)*2==(d-e)) {
if (((d-e)/p)*p==(d-e)) {
if (((d-e)/ps)*ps!=(d-e))
error[1]+=1;
}
}
else {
if (((d-e)/p)*p==(d-e))
error[4]+=1;
}
if (((d+e)/2)*2==(d+e)) {
if (((d+e)/p)*p==(d+e)) {
if (((d+e)/pc)*pc!=(d+e))
error[1]+=1;
}
}
else {
if (((d+e)/p)*p==(d+e)) {
if (((d+e)/ps)*ps==(d+e))
error[3]+=1;
flag=0;
if (((((d+e)/p)-1)/ps)*ps==(((d+e)/p)-1))
flag=1;
if (((((d+e)/p)+1)/ps)*ps==(((d+e)/p)+1))
flag=1;
if (flag==0)
error[3]+=1;
}
}
}
else {
if (((d+e)/2)*2==(d+e)) {
if (((d+e)/p)*p==(d+e)) {
if (((d+e)/ps)*ps!=(d+e))
error[1]+=1;
}
}
else {
if (((d+e)/p)*p==(d+e))
error[4]+=1;
}
if (((d-e)/2)*2==(d-e)) {
if (((d-e)/p)*p==(d-e)) {
if (((d-e)/pc)*pc!=(d-e))
error[1]+=1;
}
}
else {
if (((d-e)/p)*p==(d-e)) {
if (((d-e)/ps)*ps==(d-e))
error[3]+=1;
flag=0;
if (((((d-e)/p)-1)/ps)*ps==(((d-e)/p)-1))
flag=1;
if (((((d-e)/p)+1)/ps)*ps==(((d-e)/p)+1))
flag=1;
if (flag==0)
error[3]+=1;
}
}
}
output[n]=d;
output[n+1]=e;
n=n+2;
askip:temp=0;
}
}
bskip:
output[n]=-1;
fprintf(Outfp," error0=%d error1=%d error2=%d error3=%d error4=%d \n",
error[0],error[1],error[2],error[3],error[4]);
fprintf(Outfp," count=%d \n",(n+1)/2);
for (i=0; i<(n+1)/2; i++)
fprintf(Outfp," %#10x, %#10x, \n",output[2*i],output[2*i+1]);
fclose(Outfp);
if ((error[1]!=0)||(error[2]!=0)||(error[3]!=0)||(error[4]!=0))
printf(" error \n");
return(0);
}