/* lp04.c -- lattice point triangle with integer sides by Michael Somos */
/* use "gcc -o lp04 lp04.c" to compile program */

/* Sun, Sep  8, 1997 Cleveland State University <somos@cis.csuohio.edu> */
/* Wed, Aug 22, 2007 fix subtle bug in fsqrt() for(;;) sequencing code */

/* This program steps through all triangles with integer sides and
   selects those with integer areas and no common side factor. Call these
   "reduced Heronian triangles". Prints out an analysis of the integers
   constructed from the side lengths.
*/
#include <stdio.h>
#define AMAX 300

/*------ simple subroutines --------*/
/* sq2(x,y) = x*x + y*y */
int sq2(int x,int y) { return x*x+y*y ; }

/* max(x,y) = maximum of x and y integers */
int max(int x,int y) { return (x>y)?x:y ; }

/* isqrt(x) = floor(sqrt(x)) where x>0 is an integer */
int isqrt(int x) { int y,t; for(y=(1+x)/2;t=(y+x/y)/2,t<y;y=t); return y; }

/* fsqrt(x) = x^.5 where x>0 is a double */
double fsqrt(double x)
 {double y,t;for(y=(1+x)/2;t=(y+x/y)/2,t<y;y=t);return y;}

/* gcd(x,y) = greatest common divisor of two integers */
int gcd(int x,int y) { int t; while (y!=0) { t= x%y ;x=y;y=t;} return abs(x); }

/* ggcd(ar,ai,br,bi,&cr,&ci) sets GCD of two gaussian integers */
void ggcd(int ar,int ai,int br,int bi,int *cr,int *ci) { int t,tr,ti,ur,ui;
while (br!=0 || bi!=0) { t=br*br+bi*bi; tr=ar*br+ai*bi; ti=ai*br-ar*bi;
tr=( tr+ t*( 1-2*(tr<0) )/2 )/t; ti=( ti+ t*( 1-2*(ti<0) )/2 )/t;
ur=ar-tr*br+ti*bi; ui=ai-tr*bi-ti*br; ar=br;ai=bi; br=ur;bi=ui; }
while (ar<=0 || ai<0) { t=-ai; ai=ar; ar=t; } *cr=ar; *ci=ai; }

int main() {
double dd;
int n,at,d,D;
int a1,a2,a3; /* side lengths */
int s0,s1,s2,s3; /* semiperimeter lengths */
int p1,p2,p3,q1,q2,q3; /* angle numerator and denominators */
int W1,W2,W3;
int S0,S1,S2,S3;
int U1,U2,U3,V1,V2,V3;
int m1,m2,m3,n1,n2,n3;
int T; double R;
int tx,ty,Mx,My;
int E1,E2,E3;
int x1,x2,x3,y1,y2,y3; /* side vectors */
int X1,X2,X3,Y1,Y2,Y3; /* vertex coordinates */
n = 0;
for ( a3=1; a3 <= AMAX; a3++)
 for ( a2=1+(a3/2); a2 <= a3; a2++)
  for ( a1=a3-a2+1; a1 <= a2; a1++)
   if ( gcd ( gcd ( a3 , a2 ) , a1 ) == 1 ) {
   dd=((a1+a2+a3)*(double)(a1+a2-a3))*((a1+a3-a2)*(double)(a2+a3-a1));
   d=(.5+fsqrt(dd)); if (d*(double)d==dd) {
     /* (a1,a2,a3) is a reduced Heronian triangle */
     n = n+1;
     s0=(a1+a2+a3)/2; s1=(-a1+a2+a3)/2; s2=(a1-a2+a3)/2; s3=(a1+a2-a3)/2;
     D=d/4;
     V1=gcd(s0,s1); U1=gcd(s2,s3); E1=U1*V1; W1=a1/E1;
     V2=gcd(s0,s2); U2=gcd(s1,s3); E2=U2*V2; W2=a2/E2;
     V3=gcd(s0,s3); U3=gcd(s1,s2); E3=U3*V3; W3=a3/E3;
     S0=isqrt(s0/(V1*V2*V3));
     S1=isqrt(s1/(V1*U2*U3));
     S2=isqrt(s2/(U1*V2*U3));
     S3=isqrt(s3/(U1*U2*V3));
     T=S0*S1*S2*S3; R=W1*W2*W3/(4.0*T);
     q1=V1*S0*S1; q2=V2*S0*S2; q3=V3*S0*S3;
     p1=U1*S2*S3; p2=U2*S1*S3; p3=U3*S1*S2;
     ggcd(q3,-p3,q2,p2,&m1,&n1);
     ggcd(q1,-p1,q3,p3,&m2,&n2);
     ggcd(q2,-p2,q1,p1,&m3,&n3);
     if (n2*m3 > m2*n3) at=1; else if (n3*m1 > m3*n1) at=2; else at=3;
     x1 = E1*(m1*m1-n1*n1) ; y1 = E1*(2*m1*n1) ; if (at==1) {x1=-x1;y1=-y1;}
     x2 = E2*(m2*m2-n2*n2) ; y2 = E2*(2*m2*n2) ; if (at==2) {x2=-x2;y2=-y2;}
     x3 = E3*(m3*m3-n3*n3) ; y3 = E3*(2*m3*n3) ; if (at==3) {x3=-x3;y3=-y3;}
     if (x3>=0 && x2<0) { X1=0 ; X2=x3 ; X3=-x2 ; }
     if (y3>=0 && y2<0) { Y1=0 ; Y2=y3 ; Y3=-y2 ; }
     if (x1>=0 && x3<0) { X1=-x3 ; X2=0 ; X3=x1 ; }
     if (y1>=0 && y3<0) { Y1=-y3 ; Y2=0 ; Y3=y1 ; }
     if (x2>=0 && x1<0) { X1=x2 ; X2=-x1 ; X3=0 ; }
     if (y2>=0 && y1<0) { Y1=y2 ; Y2=-y1 ; Y3=0 ; }
     Mx = max ( X1 , max ( X2 , X3 ) ) ; My = max ( Y1 , max ( Y2 , Y3 ) ) ;

     printf("%3d %3d %3d %3d",n,a1,a2,a3);
     printf(" : %3d %5d",2*s0,D);

     printf(" : %3d %3d %3d %3d %3d %3d",X1,Y1,X2,Y2,X3,Y3);
/*   printf(" : %3d %3d %3d %3d %3d %3d",p1,q1,p2,q2,p3,q3); */
     printf(" : %2d %2d %2d %2d %2d %2d",n1,m1,n2,m2,n3,m3);

     printf(" : %d",at);
     printf("\n");
   } /* end if d*d==dd */
  } /* end gcd = 1 */
return 0;
} /* end main */