/*****************************************************************************/
/*									     */
/*  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.	The "weak" Furtwangler conditions must be fulfilled. */
/*  p is set to 3.							     */
/*									     */
/*  Note:  The maximum value of d**2-d*e+e**2 is (3/4)*d**2 if d is even or  */
/*  (3/4)*(d**2-1)+1 if d is odd.					     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0a.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
void sum(unsigned int *addend, unsigned int *augend);
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void quotient(unsigned int *a, unsigned int *b, unsigned int);
void midprod(unsigned int a, unsigned int b, unsigned int *output);
unsigned int dloop(unsigned int a, unsigned int b, unsigned int c,
		   unsigned int d, unsigned int *e);

int main ()
{
unsigned int dbeg=1000000;    // starting "a" value
unsigned int tflag=0;	   // leave set to 0

extern unsigned short table3[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int d;
extern unsigned int e;
unsigned int t3size=2556;
unsigned int outsiz=1999*3;
unsigned int n=0;
unsigned int i,j,k,l;
unsigned int T[2],V[2],X[3];
double croot2,croot4;
FILE *Outfp;
Outfp = fopen("output.dat","w");
croot2=1.259921;
croot4=1.587401;
/*****************************************/
/* find minimum index into look-up table */
/*****************************************/
midprod(dbeg, dbeg, V);
bigprod(V[0], V[1], 3, X);
V[0]=X[1];
V[1]=X[2];
quotient(V, V, 4);
if (tflag!=0) {
   T[0]=0;
   T[1]=1;
   quotient(V, V, 7);
   sum(T, V);
   quotient(V, V, 7);
   sum(T, V);
   quotient(V, V, 7);
   sum(T, V);
   if (tflag==2) {
      quotient(V, V, 7);
      sum(T, V);
      quotient(V, V, 7);
      sum(T, V);
      quotient(V, V, 7);
      sum(T, V);
      }
   }
if (V[0]==0)
   l = 32 - lmbd(1, V[1]);
else
   l = 64 - lmbd(1, V[0]);
j=l-(l/3)*3;
l=l/3;
l = 1 << l;
if (j==1)
   l=(int)(((double)(l))*croot2);
if (j==2)
   l=(int)(((double)(l))*croot4);
l=l+1;
if (l>table3[t3size-1]) {
   error[0]=5;
   goto bskip;
   }
else {
   k=0;
   for (i=0; i<t3size; i++) {
      if (table3[i] < l) k=i;
      else break;
      }
   }
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
d=0;
e=0;
error[0]=0;	// clear error array
n=dloop(dbeg, k, n, outsiz, output);
bskip:
output[n]=0xffffffff;
fprintf(Outfp," count=%d \n",(n+1)/3);
for (i=0; i<(n+1)/3; i++)
   fprintf(Outfp," %#10x %#10x %#10x \n",output[3*i],output[3*i+1],
	   output[3*i+2]);
fclose(Outfp);
return(0);
}