﻿ proposition 16
```/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/29/06 (dkc)							     */
/*									     */
/*  This C program finds a and b such that [(a**p+b**p)/(a+b)] and	     */
/*  [(a**p-b**p)/(a-b)] are primes of the form p**2*k+1 and p**2 divides     */
/*  a, b, a-b, or a+b if p>3, or p**2 divides a or b, or p**3 divides a-b or */
/*  a+b if p=3.  Whether 2 is a pth power modulo [(a**p+b**p)/(a+b)]	     */
/*  or [(a**p-b**p)/(a-b)] is determined.  Whether [(a**p+b**p)/(a+b)] is    */
/*  a pth power modulo [(a**p-b**p)/(a-b)] and whether [(a**p-b**p)/(a-b)]   */
/*  is a pth power modulo [(a**p+b**p)/(a+b)] is determined.		     */
/*									     */
/*  The output is "(a<<16)|b, (flag0<<16)|flag1".  "flag0" consists of two   */
/*  bits; if the most significant bit is set, 2 is a pth power modulo	     */
/*  [(a**p-b**p)/(a-b)] and if the least significant bit is set, 2 is a pth  */
/*  power modulo [(a**p+b**p)/(a+b)].  "flag1" consists of two bits; if      */
/*  the most significant bit is set, [(a**p+b**p)/(a+b)] is a pth power      */
/*  modulo [(a**p-b**p)/(a-b)] and if the least significant bit is set,      */
/*  [(a**p-b**p)/(a-b)] is a pth power modulo [[(a**p+b**p)/(a+b)].  If      */
/*  "flag0" is not equal to "flag1", an error is indicated ("error[1]" is    */
/*  is set to a non-zero value).					     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table3b.h"
void hugeprod(unsigned int a0, unsigned int a2, unsigned int a4,
unsigned int a6, unsigned int *product, unsigned int m0);
void dummy(unsigned int a, unsigned int b, unsigned int c);
unsigned int lmbd(unsigned int mode, unsigned int a);
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);
void hugeres(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
unsigned int *e, unsigned int f, unsigned int g);
void bigbigs(unsigned int *a, unsigned int *b);
void bigbigd(unsigned int *a, unsigned int *b);
void bigbigq(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
unsigned int *e, unsigned int f, unsigned int g);

int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 65000.
//	 The maximum "dbeg" value for p=5 is about 65000.
//	 The maximum "dbeg" value for p=7 is about 1000.
//	 The maximum "dbeg" value for p=11 is about 50.
//
unsigned int p=3;	    // prime value
unsigned int dbeg=1000;	    // starting "a" value
unsigned int dend=1;	    // ending "a" value
unsigned int rflag=1;	    // if set, p**2 must divide d, e, d+e, or d-e
//unsigned int stop=0xf4;

extern unsigned short table[];	 // prime look-up table
extern unsigned int tmptab[];	 // extended prime look-up table
extern unsigned int output[];	 // output array
extern unsigned int error[];	 // error array
unsigned int tmpsav;		 // size of extended prime look-up table
extern unsigned int count;	 // output count
unsigned int tsize=1228;	 // size of prime look-up table
unsigned int tmpsiz;		 // size of extended prime look-up table
unsigned int outsiz=1000*2;	 // size of output array
unsigned int d,e;
unsigned int a,b,temp;
unsigned int i,j,k,l;
unsigned int flag,sumdif,ps,pc;
unsigned int tflag,sflag;
unsigned int P[4],Q[2],R[2],S[2],T[2],U[2],X[3],Y[4],Z[4];
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out16.dat","w");
error[0]=0;	// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
/*********************************/
/*  extend prime look-up table	 */
/*********************************/
//
// copy prime look-up table
//
for (i=0; i<tsize; i++) tmptab[i] = (unsigned int)(table[i]);
tmpsiz=tsize;
//
// check for small prime factors
//
for (d=10001; d<160000; d++) {
if(d==(d/2)*2) continue;
if(d==(d/3)*3) continue;
if(d==(d/5)*5) continue;
if(d==(d/7)*7) continue;
if(d==(d/11)*11) continue;
if(d==(d/13)*13) continue;
if(d==(d/17)*17) continue;
if(d==(d/19)*19) continue;
//
// find size of largest possible prime factor
// find maximum index into prime look-up table
//
l = (int)(100.0 + sqrt((double)d));
k=0;
if (l>table[tsize-1]) {
error[0]=1;
goto bskip;
}
else {
for (i=0; i<tsize; i++) {
if (table[i]<l) k=i;
else break;
}
}
flag=1;
l=k;      // maximum index needed
//
//  check for prime factors using look-up table
//
for (i=0; i<=l; i++) {
k=table[i];
if ((d/k)*k==d) {
flag=0;
break;
}
}
if (flag==1) {
tmptab[tmpsiz]=d;
tmpsiz=tmpsiz+1;
}
}
tmpsav=tmpsiz;
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p; 	// p**2
pc=ps*p;	// p**3
count=0;
for (d=dbeg; d>=dend; d--) {	    // loop through d values
for (e=d-1; e>0; e--) {	    // loop through e values
//    if (e!=stop) continue;
/*********************************************/
/*  check if p**2 divides d, e, d+e or d-e   */
/*********************************************/
if (rflag!=0) {
if ((d/ps)*ps==d)	     // check if p**2 divides d
goto zskip;
if ((e/ps)*ps==e)	     // check if p**2 divides e
goto zskip;
if (p!=3) {
if (((d+e)/ps)*ps==(d+e))	// check if p**2 divides d+e
goto zskip;
if (((d-e)/ps)*ps!=(d-e))	// check if p**2 divides d-e
continue;
}
else {
if (((d+e)/pc)*pc==(d+e))	// check if p**3 divides d+e
goto zskip;
if (((d-e)/pc)*pc!=(d-e))	// check if p**e divides d-e
continue;
}
}
/****************************************/
/* check for common factors of d and e	*/
/****************************************/
zskip: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;
if((d==(d/11)*11)&&(e==(e/11)*11)) continue;
if((d==(d/13)*13)&&(e==(e/13)*13)) continue;
if((d==(d/17)*17)&&(e==(e/17)*17)) continue;
if((d==(d/19)*19)&&(e==(e/19)*19)) 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)  */
/*******************************/
tflag=0;
sflag=0;
Y[0]=0;
Y[1]=0;
Y[2]=0;
Y[3]=d;
for (i=0; i<p-1; i++)
hugeprod(Y[0], Y[1], Y[2], Y[3], Y, d);  // multiple-word product
Z[0]=0;
Z[1]=0;
Z[2]=0;
Z[3]=e;
for (i=0; i<p-1; i++)
hugeprod(Z[0], Z[1], Z[2], Z[3], Z, e);  // multiple-word product
P[0]=Z[0];
P[1]=Z[1];
P[2]=Z[2];
P[3]=Z[3];
bigbigs(Y, Z);				  // multiple-word sum
bigbigd(Y, P);				  // multiple-word difference
temp=d+e;
if (((d+e)/p)*p==(d+e))
temp=temp*p;
bigbigq(Z[0], Z[1], Z[2], Z[3], Y, 0, temp);  // multiple-word quotient
R[0]=Y[2];
R[1]=Y[3];
temp=d-e;
if (((d-e)/p)*p==(d-e))
temp=temp*p;
bigbigq(P[0], P[1], P[2], P[3], Y, 0, temp);  // multiple-word quotient
Q[0]=Y[2];
Q[1]=Y[3];
/***********************/
/*  begin inner loop   */
/***********************/
for (sumdif=0; sumdif<2; sumdif++) {
if (sumdif==0) {		// select between q and r
S[0]=Q[0];
S[1]=Q[1];
}
else {
S[0]=R[0];
S[1]=R[1];
}
/*********************************************************/
/* find size of largest possible prime factor		 */
/* find maximum index into extended prime look-up table  */
/*********************************************************/
if (S[0]==0) {
l = (33 - lmbd(1, S[1]))/2;     // left-most bit detection
l = 1 << l;
}
else {
l = (65 - lmbd(1, S[0]))/2;     // left-most bit detection
l = 1 << l;
}
k=0;
if (l>tmptab[tmpsiz-1]) {
k=tmpsiz-1;
}
else {
for (i=0; i<tmpsiz; i++) {
if (tmptab[i] < l) k=i;
else break;
}
}
l=k;	    // maximum index needed
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
for (i=0; i<=l; i++) {
k = tmptab[i];			     // load prime
quotient(S, T, k);			     // multiple-word quotient
bigprod(T[0], T[1], k, X);		     // multiple-word product
if ((S[0]==X[1]) && (S[1]==X[2])) goto askip;  // check if factor
}
/***********************************************************/
/*  continue looking for factors (not necessarily prime)   */
/***********************************************************/
//
// find size of largest possible factor
//
if (S[0]==0) {
j = (33 - lmbd(1, S[1]))/2;      // left-most bit detection
j = 1 << j;
}
else {
j = (65 - lmbd(1, S[0]))/2;      // left-most bit detection
j = 1 << j;
}
//
// loop through potential factors
//
for (i=tmptab[tmpsiz-1]; i<j; i+=2) {
quotient(S, T, i);		      // multiple-word quotient
bigprod(T[0], T[1], i, X);	      // multiple-word product
if ((X[1]==S[0]) && (X[2]==S[1])) // check if factor
}
/****************************************************************************/
/* check if [(a**p+b**p)/(a+b)] or [(a**p-b**p)/(a-b)] is of form p**2*k+1  */
/****************************************************************************/
T[0]=0;
T[1]=1;
differ(S, T);			      // multiple-word difference
quotient(T, T, ps);		      // multiple-word quotient
bigprod(T[0], T[1], ps, X);	      // multiple-word product
T[0]=0;
T[1]=1;
differ(S, T);			      // multiple-word difference
if ((X[1]!=T[0]) || (X[2]!=T[1]))    // check if p**2 is factor
/***********************************/
/* check if proposition satisfied  */
/***********************************/
quotient(T, T, p);			   // (q-1)/p
hugeres(T[0], T[1], S[0], S[1], U, 0, 2); // 2**((q-1)/p) modulus q
if ((U[0]!=0)||(U[1]!=1))	   // check if 2**((q-1)/p)==1(mod q)
tflag=tflag|1;		   // set flag
if (sumdif==0)
hugeres(T[0], T[1], S[0], S[1], U, R[0], R[1]); // least residue
else
hugeres(T[0], T[1], S[0], S[1], U, Q[0], Q[1]); // least residue
if ((U[0]!=0)||(U[1]!=1))	    // check if r**((q-1)/p)=1(mod q)
sflag=sflag|1;		    // set flag
if (sumdif==0) {
tflag=tflag<<1;		    // left-shift flag for packing
sflag=sflag<<1;		    // left-shift flag for packing
}
else { 			    // output results
if (n+1>outsiz) {
error[0]=6;
goto bskip;
}
output[n]=(d<<16) | e;	    // output d and e
output[n+1]=(tflag<<16)|sflag;  // output packed flags
n=n+2;
count=count+1;		    // increment count
if (tflag!=sflag)
error[1]=8;
}
}
}
}
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d \n",error[0],error[1]);
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)
printf(" error \n");
return(0);
}
```