﻿ Vandiver's theorem
```/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/23/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.	[(a**p+b**p)/(a+b)] must have two distinct prime     */
/*  factors.  p is set to 3.  The prime factors of [(a**p+b**p)/(a+b)] must  */
/*  be of the form p**2*k+1.						     */
/*									     */
/*  The output is "a, b, split".  If 2p does not divide a, b, a-b, or a+b,   */
/*  then "split" is set to 1, otherwise "split" is set to 0.                 */
/*									     */
/*  If a is even, p does not divide a, and (a/2)**(p-1) is not congruent to  */
/*  1 modulus p**3, then an error is indicated ("error[1]" is set to 9).  b  */
/*  is treated similarly.						     */
/*									     */
/*  If a is odd, p divides a+b, and a**(p-1) is not congruent to 1 modulus   */
/*  p**3, then an error is indicated ("error[1]" is set to 9).  b and a-b    */
/*  are treated similarly.						     */
/*									     */
/*  If p**3 divides a, b, a-b or p**4 divides a+b, 2 does not divide a, and  */
/*  a**(p-1) is not congruent to 1 modulus p**3, then an error is indicated  */
/*  ("error[1]" is set to 8). b and a-b are treated similarly.  If p**3      */
/*  divides a, b, or a-b or p**4 divides a+b, 2 does not divide a+b, p does  */
/*  not divide a+b, and (a+b)**(p-1) is not congruent to 1 modulus p**3,     */
/*  then an error is indicated.  If p**3 divides a, b, or a-b or p**4	     */
/*  divides a+b, 2 does not divide a+b, p divides a+b, and [(a+b)/p]**(p-1)  */
/*  is not congruent to 1 modulus p**3, then an error is indicated.	     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0b.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
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);
int main ()
{
unsigned int p=3;	  // input prime
unsigned int dbeg=10000;  // starting "a" value
unsigned int dend=1;	  // ending "a" value
//unsigned int stop=4160;

extern unsigned int table3[];
extern unsigned int output[];
extern unsigned int error[];
unsigned int t3size=848;
unsigned int outsiz=1999;
unsigned int n=0;
unsigned int d,e,a,b,temp,dsum;
unsigned int i,j,k,l,m;
unsigned int flag,split,qflag,ps,pc,pf;
unsigned int S[2],T[2],V[2],X[3];
double recip7,croot2,croot4;
FILE *Outfp;
Outfp = fopen("out28b.dat","w");
recip7=1.0/7.0;
croot2=1.259921/7.0;
croot4=1.587401/7.0;
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p;
pc=ps*p;
pf=pc*p;
error[0]=0;	// clear error array
error[1]=0;
for (d=dbeg; d>=dend; d--) {
for (e=d-1; e>0; e--) {
dummy(d,e,0);
//    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;
/*******************************/
/* check if 2 and p are split  */
/********************************/
split=1;
if ((d/2)*2==d)
if ((d/p)*p==d)
split=0;
if ((e/2)*2==e)
if ((e/p)*p==e)
split=0;
if (((d-e)/2)*2==(d-e))
if (((d-e)/p)*p==(d-e))
split=0;
if (((d+e)/2)*2==(d+e))
if (((d+e)/p)*p==(d+e))
split=0;
/**********************************************************/
/* check if p**3 divides a, b, or a-b or p**4 divides a+b */
/**********************************************************/
qflag=0;
if ((d/pc)*pc==d)
qflag=1;
if ((e/pc)*pc==e)
qflag=1;
if (((d+e)/pf)*pf==(d+e))
qflag=1;
if (((d-e)/pc)*pc==(d-e))
qflag=1;
/************************************/
/*  compute (d**p + e**p)/(d + e)   */
/************************************/
dsum=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];
}
sum(S, T);
quotient(T, S, dsum);
/************************************************/
/*  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)
l=(int)(((double)(l))*recip7);
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;
}
}
j=0;
for (i=0; 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)
else {
j=i+1;
break;
}
}
if ((m/3)*3!=m) continue;
if ((S[0]==X[1]) && (S[1]==X[2])) continue;
if ((S[0]==0)&&(S[1]==1))
continue;
m=0;
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;
bloop:	 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 bloop;
if ((m/3)*3!=m)
else
break;
}
if ((m/3)*3!=m) continue;
if ((S[0]!=0) || (S[1]!=1)) continue;
//
// p**3
//
if (qflag!=0) {
if ((d/2)*2!=d) {
flag=0;
if (((d-1)/pc)*pc==(d-1))
flag=1;
if (((d+1)/pc)*pc==(d+1))
flag=1;
if (flag==0)
error[1]=8;
}
if ((e/2)*2!=e) {
flag=0;
if (((e-1)/pc)*pc==(e-1))
flag=1;
if (((e+1)/pc)*pc==(e+1))
flag=1;
if (flag==0)
error[1]=8;
}
if (((d-e)/2)*2!=(d-e)) {
flag=0;
if ((((d-e)-1)/pc)*pc==((d-e)-1))
flag=1;
if ((((d-e)+1)/pc)*pc==((d-e)+1))
flag=1;
if (flag==0)
error[1]=8;
}
if (((d+e)/2)*2!=(d+e)) {
if (((d+e)/p)*p==(d+e)) {
flag=0;
if (((((d+e)/p)-1)/pc)*pc==(((d+e)/p)-1))
flag=1;
if (((((d+e)/p)+1)/pc)*pc==(((d+e)/p)+1))
flag=1;
if (flag==0)
error[1]=8;
}
else {
flag=0;
if ((((d+e)-1)/pc)*pc==((d+e)-1))
flag=1;
if ((((d+e)+1)/pc)*pc==((d+e)+1))
flag=1;
if (flag==0)
error[1]=8;
}
}
}
//
// "split" tests
//
if ((d/2)*2==d) {
if ((d/p)*p!=d) {
flag=0;
if ((((d/2)-1)/pc)*pc==((d/2)-1))
flag=1;
if ((((d/2)+1)/pc)*pc==((d/2)+1))
flag=1;
if (flag==0)
error[1]=9;
}
}
else {
if (((d+e)/p)*p==(d+e)) {
flag=0;
if (((d-1)/pc)*pc==(d-1))
flag=1;
if (((d+1)/pc)*pc==(d+1))
flag=1;
if (flag==0)
error[1]=9;
}
}
if ((e/2)*2==e) {
if ((e/p)*p!=e) {
flag=0;
if ((((e/2)-1)/pc)*pc==((e/2)-1))
flag=1;
if ((((e/2)+1)/pc)*pc==((e/2)+1))
flag=1;
if (flag==0)
error[1]=9;
}
}
else {
if (((d+e)/p)*p==(d+e)) {
flag=0;
if (((d-1)/pc)*pc==(d-1))
flag=1;
if (((d+1)/pc)*pc==(d+1))
flag=1;
if (flag==0)
error[1]=9;
}
}
if (((d-e)/2)*2!=(d-e)) {
if (((d+e)/p)*p==(d+e)) {
flag=0;
if ((((d-e)-1)/pc)*pc==((d-e)-1))
flag=1;
if ((((d-e)+1)/pc)*pc==((d-e)+1))
flag=1;
if (flag==0)
error[1]=9;
}
}
if (n+2>outsiz) {
error[0]=6;
goto bskip;
}
output[n]=d;
output[n+1]=e;
output[n+2]=split;
n=n+3;
}
}
bskip:
output[n]=-1;
fprintf(Outfp," error0=%d error1=%d \n",error[0],error[1]);
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);
if (error[1]!=0)
printf(" error \n");
return(0);
}
```