DiophantineQuadraticEquation.java

// http://www.alpertron.com.ar/QUAD.HTM
// <XMP>
// Diophantine Quadratic Equation Solver
// ax^2 + bxy + cy^2 + dx + ey + f = 0 (unknowns x,y integer numbers)
//
// Written by Dario Alejandro Alpern (Buenos Aires - Argentina)
// Last updated December 15th, 2003
//
// No part of this code can be used for commercial purposes without
// the written consent from the author. Otherwise it can be used freely.

import java.applet.*;
import java.util.*;
import java.awt.*;
import java.math.*;

public final class quad extends Applet implements Runnable {

private BigInteger Primes[];
private int Exponents[];
private BigInteger PrimesBak[];
private int ExponentsBak[];
private int digitsInGroup;
private Vector sortedSolsX = new Vector(50, 50);
private Vector sortedSolsY = new Vector(50, 50);
private boolean allSolsFound;
private long A, B, C, D, E, F;
private long Xi;
private long Xl;
private long Yi;
private long Yl;
private long CY1, CY0;
private boolean also, ExchXY, teach;
private long SQD;
private long NUM[] = new long[6];
private long DEN[] = new long[6];
private long DET;
private long Mi=1000000000;
private long Bi=Mi*Mi;
private long DosALa32 = (long)1 << 32;
private long DosALa31 = (long)1 << 31;
private double dDosALa32 = (double)DosALa32;
private double dDosALa64 = dDosALa32 * dDosALa32;
private double dDosALa96 = dDosALa64 * dDosALa32;
private double dDosALa128 = dDosALa96 * dDosALa32;
private double dDosALa160 = dDosALa128 * dDosALa32;
private long DosALa32_1 = DosALa32 - 1;
private long A1;
private long A2;
private long B1;
private long B2;
private String UU="";
private String VU="";
private String UL="";
private String VL="";
private String UL1="";
private String VL1="";
private String FP="";
private String msg="There are no solutions !!!",sq;
private String txt="";
private String divgcd="Dividing the equation by the greatest common divisor we obtain:<BR>";
private long CX2;
private long CXY;
private long CY2;
private long CX;
private long CY;
private long C1;
private long H1[];
private long H2[];
private long K1[];
private long K2[];
private long L1[];
private long L2[];
private int NbrSols, NbrCo, NbrEqs, EqNbr;
private StringBuffer info;
private BigInteger ValA, ValB, ValC, ValD, ValE, ValF;
private volatile Thread calcThread;
private String textError;

private void SolveEquation() {
  byte b;
  long Fact1, Fact2, Tmp;
  long biA[] = new long[6];
  long biB[] = new long[6];
  long biC[] = new long[6];
  boolean teachaux;
  long D1, E1, F1, G, H, K, P, P1, P2, Q, Q1, Q2, R, R1, R2, X1, Y1;
  long S, S1, S2, T, g, r, s, t, u, w1, w2;
  String t1,x,y,x1,y1;
  CX2=ValA.longValue();
  CXY=ValB.longValue();
  CY2=ValC.longValue();
  CX=ValD.longValue();
  CY=ValE.longValue();
  C1=ValF.longValue();
  also = ExchXY = false;
  sortedSolsX.removeAllElements();   /* Start fresh for new equation */
  sortedSolsY.removeAllElements();
  allSolsFound = false;
  NbrCo = -1;
  w("<TITLE>"+txt+"</TITLE><CENTER><B>");
  A=CX2;B=CXY;C=CY2;D=CX;E=CY;F=C1;
  ShowEq(A,B,C,D,E,F,"x","y");
  w(" = 0</B></CENTER><P><I>by Dario Alejandro Alpern</I><P>");
  t=gcd(A,gcd(B,gcd(C,gcd(D,E))));
  if (teach) {
    w("First of all we must determine the <B>gcd</B> of all coefficients but the constant term, that is: <B>gcd</B>("+A+", "+B+", "+C+", "+D+", "+E+") = "+t+".<P>");
    }
  if (t!=0) {
    if (F%t!=0) {
      NoGcd(F);return;
      }
    else {
      A/=t;B/=t;C/=t;D/=t;E/=t;F/=t;
      if (teach) {
        w(divgcd);
        ShowEq(A,B,C,D,E,F,"x","y");
        w(" = 0</B><P>");
        }
      }
    }
  if (D==0 && A!=0 && C!=0) {
    if (CheckMod(A,B,C,E,F)) {
      return;
      }
    }
  if (E==0 && A!=0 && C!=0) {
    if (CheckMod(B,C,A,D,F)) {
      return;
      }
    }
  S=B*B-4*A*C;
  if (S>0 && sqrt(S)*sqrt(S)!=S && D==0 && E==0 && F!=0)
    {
    teachaux = teach;
    if (abs(F)!=1) {teach = false;}
    LongToDoublePrecLong(A, biA);
    LongToDoublePrecLong(B, biB);
    LongToDoublePrecLong(C, biC);
    GetRoot(biA, biB, biC);
    teach = teachaux;
    G=H=F;K=1;T=3;
    while(G%4==0) {
      G/=4;K*=2;
      }
    while(abs(G)>=T*T) {
      while(G%(T*T)==0) {
        G/=T*T;K*=T;
        }
      T+=2;
      }
    for (T=1;T*T<=K;T++) {
      if (K%T==0) {
        SolContFrac(H,T,A,B,C,"");
        }
      }
    for (T=T-1;T>0;T--) {
      if (K%T==0 && T*T<K) {
        SolContFrac(H,K/T,A,B,C,"");
        }
      }
    w("<P>");
    if (also==false) {
      return;
      }
    if (teach) {
      w("<TABLE BORDER=1><TR><TH>");
      }
    else {
      ShowAllLargeSolutions();
      }
    if (teach==false) {
      w("If <B>(x,y)</B> is a solution, <B>(-x,-y)</B> is also a solution.<P>");
      }
    ShowRecursion((byte)0);
    also=true;
    return;
    }
  if (A==0 && C==0) {
    if (B==0) {
      b=Linear(D,E,F);
      if (teach) {
        w("<TABLE BORDER=1><TR><TH>");
        }
      PrintLinear(b,"t");
      if (teach) {
        w("</TABLE>");
        }
      return;
      }
    else {
      R=D*E-B*F;
      if (teach) {
        w("Multiplying by "+B+" we obtain:<BR>");
        ShowEq(0,B*B,0,D*B,E*B,0,"x","y");
        w(" = "+(-F*B)+"<P>");
        w("Adding "+D*E+" to both sides of the equal sign:<BR>");
        ShowEq(0,B*B,0,D*B,E*B,D*E,"x","y");
        w(" = "+R+"<P>");
        w("Now the left side can be factored as follows:<BR>");
        w("(");
        Show1(E,Show(B," x",(byte)0));
        w(") (");
        Show1(D,Show(B," y",(byte)0));
        w(") = "+R+"<P>");
        if (R!=0) {
          w("Then ");
          Show1(E,Show(B," x",(byte)0));
          w(" must be a factor of "+R+", so we must find the factors of "+R+":<P><UL>");
          }
        else {
          w("One of the parentheses must be zero, so:<P>");
          Show1(E,Show(B," x",(byte)0));
          w(" = 0");
          if (E%B==0) {
            w(" means that <B>x = "+(-E/B)+"</B> and <B>y</B> could be any integer.<P>");
            also=true;
            }
          else {
            w("This equation cannot be solved in integers.<P>");
            }
          Show1(D,Show(B," y",(byte)0));
          w(" = 0");
          if (D%B==0) {
            w(" means that <B>y = "+(-D/B)+"</B> and <B>x</B> could be any integer.<P>");
            also=true;
            }
          else {
            w("This equation cannot be solved in integers.<P>");
            }
          return;
          }
        }
      if (R!=0) {
        S=sqrt(abs(R));
        for (T=1; T<=S; T++) {
          if (R%T==0) {
            SolByFact(R,T,B,D,E);
            SolByFact(R,-T,B,D,E);
            if (T*T!=abs(R)) {
              SolByFact(R,R/T,B,D,E);
              SolByFact(R,-R/T,B,D,E);
              }
            }
          } /* end for */
        if (teach) {
          w("</UL>");
          }
        return;
        }
      if (E%B==0) {
        Xi=-E/B;Xl=0;Yi=0;Yl=1;PrintLinear((byte)0,"t");
        }
      if (D%B==0) {
        Xi=0;Xl=1;Yi=-D/B;Yl=0;PrintLinear((byte)0,"t");
        }
      return;
      }
    }
  if (teach) {
    w("We try now to solve this equation module 9, 16 and 25.<P>");
    }
  if (Mod(9)) {return;}
  if (Mod(16)) {return;}
  if (Mod(25)) {return;}
  if (teach) {
    w("There are solutions, so we must continue.<P>");
    }
  if (A==0) {
    T=A;A=C;C=T;T=D;D=E;E=T;ExchXY=true;
    }
  D1=4*A*C-B*B;
  E1=4*A*E-2*B*D;
  F1=4*A*F-D*D;
  x=(ExchXY?"y":"x");
  y=(ExchXY?"x":"y");
  x1=x+"´";
  y1=y+"´";
  if (D1==0) {              /* Parabolic case */
    r=gcd(2*A,B);
    s=2*A/r;
    P=r/2;
    Q=D;
    R=(2*A*E-B*D)/r;
    S=2*A*F/r;
    if (teach) {
      if (s!=1) {
        w("Multiplying the equation by "+par(s)+":<BR>");
        ShowEq(A*s,B*s,C*s,D*s,E*s,F*s,x,y);
        w(" = 0<P>");
        }
      if (r!=2) {
        w("Extracting the factor "+r/2+" in the quadratic terms:<BR>");
        w(par1(r/2)+" (");
        ShowEq(s*s,2*s*B/r,2*s*C/r,0,0,0,x,y);
        w(")");
        b=Show(D*s," "+x,(byte)1);
        b=Show(E*s," "+y,b);
        Show1(F*s,b);
        w(" = 0<P>");
        }
      if (B!=0) {
        w(par1(r/2)+"(");ShowLin(s,B/r,0,x,y);w(")"+sq);
        if (D!=0 || E!=0 || F!=0) {
          w(" + (");
          ShowLin(D*s,E*s,F*s,x,y);
          w(")");
          }
        w(" = 0<P>");
        if (D!=0) {
          w("Adding and subtracting "+par1(B*D/r)+""+y+":<BR>");
          w(par1(r/2)+"(");ShowLin(s,B/r,0,x,y);w(")"+sq);
          w(" + "+par1(D)+" (");
          ShowLin(s,B/r,0,x,y);
          w(")");
          }
        if (E1!=0 || F!=0) {
          w(" + (");
          ShowLin(0,R,F*s,x,y);
          w(")");
          }
        }
      else {
        w(par1(r/2)+"(");ShowLin(s,B/r,0,x,y);w(")"+sq);
        if (D!=0) {
          w(" + "+par1(D)+" (");
          ShowLin(s,B/r,0,x,y);
          w(")");
          }
        if (E1!=0 || F!=0) {
          w(" + (");
          ShowLin(0,R,F*s,x,y);
          w(")");
          }
        }
      w(" = 0<P>Now we perform the substitution:<BR>");
      w(x1+" = ");
      ShowLin(s,B/r,0,x,y);
      w("<P>This gives:<BR>");
      ShowEq(r/2,0,0,D,R,F*s,x1,y);
      w(" = 0<P>");
      }
    if (E1==0) {
      if (teach) {
        w("This can be solved by the standard quadratic equation formula:<P>");
        w("The roots are: "+x1+" = ");
        if (D!=0) {
          w("("+par(-D)+" ± sqrt("+(-F1)+"))");
          }
        else {
          w(" ± sqrt("+(-F1)+")");
          }
        w(" / "+par(r)+"<P>");
        }
      if (F1>0) {
        if (teach) {
          w("This quadratic equation has no solution in reals, so it has no solution in integers.<P>");
          also=true;
          }
        return;
        }
      T=(long)Math.floor((-D+Math.sqrt(-F1))/r+0.5);
      if (r*T*T/2+D*T+2*A*F/r==0) {
        if (teach) {
          w("The first root is: "+x1+" = "+T+"<P>");
          ShowLin(s,B/r,0,x,y);
          w(" = "+T+"<P>");
          }
        b=Linear(s,B/r,-T);
        if (teach) {
          w("<TABLE BORDER=1><TR><TH>");
          }
        PrintLinear(b,"t");
        if (teach) {
          w("</TABLE>");
          }
        }
      else {
        w("The first root is not an integer.<P>");
        }
      if (F1==0) {
        return;
        }
      T=(long)Math.floor((-D-Math.sqrt(-F1))/r+0.5);
      if (r*T*T/2+D*T+2*A*F/r==0) {
        if (teach) {
          w("The second root is: "+x1+" = "+T+"<P>");
          ShowLin(s,B/r,0,x,y);
          w(" = "+T+"<P>");
          }
        b=Linear(s,B/r,-T);
        if (teach) {
          w("<TABLE BORDER=1><TR><TH>");
          }
        PrintLinear(b,"t");
        if (teach) {
          w("</TABLE>");
          }
        }
      else {
        w("The second root is not an integer.<P>");
        }
      return;
      }
    long N=DivideGcd(P,Q,R,S,x1,y);
    if (N==0) {
      return;
      }
    P/=N;Q/=N;R/=N;S/=N;
    int numEq=1;
    for (u=0;u<abs(R);u++) {
      T=P*u*u+Q*u+S;
      if (T%R==0) {
        int numEq2=3;
        P1=B*P*R/r;Q1=abs(R)+B*(2*P*u+Q)*abs(R)/R/r;
        R1=-s;S1=B*T/R/r+u;
        P2=-P*R;Q2=(2*P*u+Q)*abs(R)/-R;S2=-T/R;
        if (teach) {
          t1=((u==0)?"":u+" + ")+abs(R)+"t";
          w("<LI>"+x1+" = "+t1+"   <B>("+numEq+")</B>");
          w("<P>Replacing this in the equation shown above:<BR>");
          w(par1(-R)+y+" = ");
          ShowEq(P,0,0,Q,0,S,"("+t1+")",y);
          w("<P>"+par1(-R)+y+" = ");
          ShowEq(P*R*R,0,0,(2*P*u+Q)*abs(R),0,T,"t",y);
          w("<P>"+y+" = ");
          ShowEq(P2,0,0,Q2,0,S2,"t",y);
          w("   <B>("+(numEq+1)+")</B><P>From <B>("+numEq+")</B>: ");
          ShowLin(s,B/r,0,x,y);
          w(" = "+t1+"<P>Replacing <B>("+(numEq+1)+")</B> here:<BR>");
          w(par1(s)+x+" = ");
          ShowEq(P1,0,0,Q1,0,S1,"t",y);
          w("   <B>("+(numEq+2)+")</B><P>");
          }
        long N1=DivideGcd(P1,Q1,R1,S1,"t",x);
        if (N1==0) {
          continue;
          }
        P1/=N1;Q1/=N1;R1/=N1;S1/=N1;
        for (long u1=0;u1<abs(R1);u1++) {
          long T1=P1*u1*u1+Q1*u1+S1;
          if (T1%R1==0) {
            if (teach && abs(R1)!=1) {
              t1=((u1==0)?"":u1+" + ")+abs(R1)+"u";
              w("<LI>t = "+t1+"   <B>("+(numEq+numEq2)+")</B>");
              w("<P>Replacing this in <B>("+(numEq+2)+")</B>:<BR>");
              w(par1(-R1)+x+" = ");
              ShowEq(P1,0,0,Q1,0,S1,"("+t1+")",y);
              w("<P>"+par1(-R1)+x+" = ");
              ShowEq(P1*R1*R1,0,0,(2*P1*u1+Q1)*abs(R1),0,T1,"u","");
              w("<P>"+x+" = ");
              ShowEq(-P1*R1,0,0,(2*P1*u1+Q1)*abs(R1)/-R1,0,-T1/R1,"u",y);
              w("<P>From <B>("+(numEq+1)+")</B> and <B>("+(numEq+numEq2)+")</B>:<BR>"+y+" = ");
              ShowEq(P2,0,0,Q2,0,S2,"("+t1+")",y);
              w("<P>"+y+" = ");
              ShowEq(P2*R1*R1,0,0,(2*P2*u1+Q2)*abs(R1),0,P2*u1*u1+Q2*u1+S2,"u","");
              w("<P>");
              numEq2++;
              }
            if (teach) {
              w("<TABLE BORDER=1><TR><TH>");
              }
            showAlso();
            w("x = ");
            if (ExchXY==false) {
              ShowEq(-P1*R1,0,0,(2*P1*u1+Q1)*abs(R1)/-R1,0,-T1/R1,"u","");
              w("<BR>y = ");
              ShowEq(P2*R1*R1,0,0,(2*P2*u1+Q2)*abs(R1),0,P2*u1*u1+Q2*u1+S2,"u","");
              }
            else {
              ShowEq(P2*R1*R1,0,0,(2*P2*u1+Q2)*abs(R1),0,P2*u1*u1+Q2*u1+S2,"u","");
              w("<BR>y = ");
              ShowEq(-P1*R1,0,0,(2*P1*u1+Q1)*abs(R1)/-R1,0,-T1/R1,"u","");
              }
            if (teach) {
              w("</TABLE>");
              }
            w("<P>");
            }
          }
        if (teach && abs(R1)!=1) {
          w("</UL>");
          }
        numEq+=numEq2;
        }
      }
    if (teach) {
      w("</UL>");
      }
    return;
    }
  g=gcd(D1,E1/2);
  CY1=D1/g;
  CY0=E1/2/g;
  long D0=D1;
  long N0=CY0*CY0*g-CY1*F1;
  long h=gcd(CY1,gcd(g,N0));
  double sqrt = Math.sqrt((double)g*(double)N0);
  double R3=(-E1/2-sqrt)/D1;
  R1=(long)Math.ceil(R3);
  double R4=(-E1/2+sqrt)/D1;
  R2=(long)Math.floor(R4);
  if (teach) {
    w("We want to convert this equation to one of the form:<BR>");
    w(x1+sq+" + B "+y+sq+" + C "+y+" + D = 0<P>");
    w("Multiplying the equation by "+par(4*A)+":<BR>");
    ShowEq(4*A*A,4*A*B,4*A*C,4*A*D,4*A*E,4*A*F,x,y);
    w(" = 0<P>");
    ShowLin(4*A*A,0,0,x+sq,y);
    if (B!=0 || D!=0) {
      w(" + (");
      ShowLin(0,4*A*B,4*A*D,x,y);
      w(")"+x);
      }
    if (C!=0 || E!=0 || F!=0) {
      w(" + (");
      ShowEq(0,0,4*A*C,0,4*A*E,4*A*F,x,y);
      w(") = 0<P>");
      }
    if (B!=0 || D!=0) {
      w("To complete the square we should add and subtract:<BR>(");
      ShowLin(0,B,D,x,y);
      w(")"+sq+"<P>Then the equation converts to:<BR>(");
      ShowLin(2*A,B,D,x,y);
      w(")"+sq+" + (");
      ShowEq(0,0,4*A*C,0,4*A*E,4*A*F,x,y);
      w(") - (");
      ShowEq(0,0,B*B,0,2*B*D,D*D,x,y);
      w(") = 0<P>");
      }
    w("(");
    ShowLin(2*A,B,D,x,y);
    w(")"+sq+" + (");
    ShowEq(0,0,D1,0,E1,F1,x,y);
    w(") = 0<P>Now we perform the substitution:<BR>");
    w(x1+" = ");
    ShowLin(2*A,B,D,x,y);
    w("<P>This gives:<BR>");
    ShowEq(1,0,D1,0,E1,F1,x1,y);
    w(" = 0<P>");
    if (D0>0) {         /* elliptical case */
      w("Since "+x1+sq+" is always greater than, or equal to zero,<BR>");
      ShowEq(0,0,D1,0,E1,F1,x1,y);
      w(" must be less than, or equal to zero. This is verified in the segment limited by the roots.<P>");
      if (N0<0) {
        w("The polynomial in "+y+" is always positive,");
        NoSol();
        return;
        }
      w("The roots are: (-"+par(E1)+" - sqrt("+E1+sq+" - 4 * "+par(D1)+" * "+par(F1)+")) / (2 * "+par(D1)+") = "+R3+"<BR>");
      w("and: (-"+par(E1)+" + sqrt("+E1+sq+" - 4 * "+par(D1)+" * "+par(F1)+")) / (2 * "+par(D1)+") = "+R4+"<P>");
      if (R2<R1) {
        w("There are no integers in this range,");
        NoSol();
        return;
        }
      w("All values of "+y+" from "+R1+" to "+R2+" should be replaced in <BR>");
      ShowEq(0,0,D1,0,E1,F1,x1,y);
      w(". The result should be the negative of a perfect square.<P>");
      b=0;
      for (u=R1;u<=R2;u++) {
        w1=-D1*u*u-E1*u-F1;w2=sqrt(w1);
        if (w2*w2==w1) {
          if (b!=0) {
            w(", "+u);
            }
          else {
            w("The values of "+y+" are: "+u);
            b=1;
            }
          }
        }
      if (b==0) {
        w("This is not satisfied by any value of "+y);
        NoSol();
        return;
        }
      w("<P><UL>");
      }
    }
  if (D0>0 && N0<0) {
    return;
    }
  if (D0>0) {
    for (u=R1;u<=R2;u++) {
      w1=-D1*u*u-E1*u-F1;w2=sqrt(w1);
      if (w2*w2==w1) {
        if (teach) {
          w("<LI>"+y+" = "+u+"<BR>");
          w(x1+" = ");
          ShowLin(2*A,B,D,x,y);
          w(" = ±sqrt("+(w2*w2)+") = ±"+w2+"<P><UL>");
          }
        ShowElipSol(A,B,D,u,x,x1,y,w2);
        if (w2!=0) {
          ShowElipSol(A,B,D,u,x,x1,y,-w2);
          }
        if (teach) {
          w("</UL>");
          }
        }
      }
    if (teach) {
      w("</UL>");
      }
    return;
    }
  if (teach) {
    if (D1!=g*h) {
      w("Multiplying the equation by "+CY1/h+":<BR>");
      ShowEq(CY1/h,0,D1*CY1/h,0,D1*E1/g/h,D1*F1/g/h,x1,y);
      }
    w(" = 0<P>");
    if (E1!=0) {
      Show(g/h,"(",Show(CY1/h,x1+sq,(byte)0));
      ShowEq(CY1*CY1,0,0,2*CY0*CY1,0,0,y,"");
      w(")");
      Show1(D1*F1/g/h,(byte)1);
      w(" = 0<P>");
      Show(g/h,"(",Show(CY1/h,x1+sq,(byte)0));
      if (CY1 != 1) {
        w(par(CY1)+sq+" ");
        }
      w(y+sq+" + 2*");
      if (CY1 != 1) {
        w(par(CY1)+"*");
        }
      w(par(CY0)+" "+y+")");
      Show1(D1*F1/g/h,(byte)1);
      w(" = 0<P>");
      w("Adding and subtracting "+(g==h?"":g/h+" * ")+par(E1/2/g)+sq+":<BR>");
      }
    Show(g/h,"(",Show(CY1/h,x1+sq,(byte)0));
    if (CY1 != 1) {
      w(par(CY1)+sq+" ");
      }
    w(y+sq+" + 2*");
    if (CY1 != 1) {
      w(par(CY1)+"*");
      }
    w(par(CY0)+" "+y+" + "+par(CY0)+sq+")");
    Show1(D1*F1/g/h,(byte)1);
    w(" - "+(g==h?"":g/h+" * ")+par(E1/2/g)+sq+" = 0<P>");
    Show(g/h,"(",Show(CY1/h,x1+sq,(byte)0));
    ShowLin(0,CY1,CY0,x,y);
    w(")"+sq);
    Show1(-N0/h,(byte)1);
    w(" = 0<P>Making the substitution "+y1+" = ");
    ShowLin(0,CY1,CY0,x,y);
    w(":<BR>");
    ShowLin(CY1/h,g/h,-N0/h,x1+sq,y1+sq);
    w(" = 0<P>");
    }
  long Sqd=sqrt(-D0);
  long N1=abs(N0/h);
  long Xc=sqrt(abs(CY1/h));
  long Yc=sqrt(g/h);
  if (Sqd*Sqd==-D0) {
    if (teach) {
      w("(");
      ShowLin(Yc,Xc,0,y1,x1);
      w(") (");
      ShowLin(Yc,-Xc,0,y1,x1);
      w(") = "+(N0/h)+"<P>");
      }
    if (N0==0) {
      if (teach) {
        w("<P>One of the parentheses must be zero, so:<P><UL>");
        }
      for (T=(Sqd==0?1:0); T<2; T++) {
        if (teach) {
          w("<LI>");
          ShowLin(Yc,Xc,0,y1,x1);
          w(" = 0<P>"+par(Yc)+" (");
          ShowLin(0,CY1,CY0,x,y);
          w(") + "+par1(Xc)+" (");
          ShowLin(2*A,B,D,x,y);
          w(") = 0<P>");
          ShowLin(2*A*Xc,Xc*B+CY1*Yc,D*Xc+CY0*Yc,x,y);
          w(" = 0<P>");
          }
        b=Linear(2*A*Xc,Xc*B+CY1*Yc,D*Xc+CY0*Yc);
        if (teach) {
          w("<TABLE BORDER=1><TR><TH>");
          }
        PrintLinear(b,"t");
        if (teach) {
          w("</TABLE><P>");
          }
        Xc=-Xc;
        }
      Sqd=abs(Sqd);
      if (teach) {
        w("</UL>");
        }
      return;
      }
    S=sqrt(N1);
    if (teach) {
      w("Now we have to find all factors of "+N1+".<P><OL>");
      }
    for (T=1;T<=S;T++) {
      if (N1%T==0) {
        Fact1=T;
        Fact2=N0/h/T;
        if (teach) {
          w("<LI>Since "+(Fact1*Fact2)+" is equal to "+Fact1+" times "+Fact2+", we can set:<BR>");
          ShowLin(Yc,Xc,0,y1,x1);
          w(" = "+Fact1+"<BR>");
          ShowLin(Yc,-Xc,0,y1,x1);
          w(" = "+Fact2+"<BR>");
          if ((Fact1-Fact2)%(2*Xc)==0 && (Fact1+Fact2)%(2*Yc)==0) {
            X1=(Fact1-Fact2)/(2*Xc);
            Y1=(Fact1+Fact2)/(2*Yc);
            w(x1+" = "+X1+"<BR>"+y1+" = "+Y1+"<P><B>");
            ShowX1Y1(X1,Y1,A,B,D,CY1,CY0);
            w("</B>");
            }
          else {
            w("Solving this system we do not obtain integer values for "+x1+" and "+y1+".<P>");
            }
          }
        else {
          if ((Fact1-Fact2)%(2*Xc)==0 && (Fact1+Fact2)%(2*Yc)==0) {
            X1=(Fact1-Fact2)/(2*Xc);
            Y1=(Fact1+Fact2)/(2*Yc);
            ShowX1Y1(X1,Y1,A,B,D,CY1,CY0);
            }
          }
        }
      }
    if (teach) {
      w("</OL><P>");
      }
    return;
    }
  if (N0==0) {
    ShowX1Y1(0,0,A,B,D,D1,E1/2);
    return;
    }
  /* Test if we need two cycles or four cycles */
  LongToDoublePrecLong(A, biA);
  LongToDoublePrecLong(C, biB);
  MultDoublePrecLong(biA, biB, biC);
  LongToDoublePrecLong(1, biA);
  LongToDoublePrecLong(B, biB);
  GetRoot(biA, biB, biC);
  LongToDoublePrecLong(A, biA);
  ContFrac(biA,(byte)5,(byte)1,0,B*B-4*A*C,1,A); /* A2, B2 solutions */
  DET=B*B-4*A*C;
  G=(2*A2+B*B2)%DET;
  H=(B*A2+2*A*C*B2)%DET;
  if (((C*D*(G-2)+E*(B-H))%DET!=0 || (D*(B-H)+A*E*(G-2))%DET!=0) &&
      ((C*D*(-G-2)+E*(B+H))%DET!=0 || (D*(B+H)+A*E*(-G-2))%DET!=0)) { 
      NbrCo*=2;
    }
  LongToDoublePrecLong(D1/g/h, biA);
  LongToDoublePrecLong(0, biB);
  LongToDoublePrecLong(g/h, biC);
  GetRoot(biA, biB, biC);
  G=H=-N0/h;K=1;T=3;
  while(G%4==0) {
    G/=4;K*=2;
    }
  while(abs(G)>=T*T) {
    while(G%(T*T)==0) {
      G/=T*T;K*=T;
      }
    T+=2;
    }
  for (T=1;T*T<=K;T++) {
    if (K%T==0) {
      SolContFrac(H,T,D1/g/h,0,g/h,"'");
      }
    }
  for (T=T-1;T>0;T--) {
    if (K%T==0 && T*T<K) {
      SolContFrac(H,K/T,D1/g/h,0,g/h,"'");
      }
    }
  w("<P>");
  if (also==false) {
    return;
    }
  else {
    if (teach==false) {
      ShowAllLargeSolutions();
      }
    }
  ShowRecursion((byte)1);
  also=true;
  return;
  }
private void NoSol() {
   w(" so there are no integer solutions.<P>");
   also=true;
   }

private void NoGcd(long F) {
  if (teach) {
    w("This <B>gcd</B> is not a divisor of the constant term ("+F+"),");
    NoSol();
    }
  }

private void w(String texto) {
  info=info.append(texto);
  }

private long gcd(long M,long N) {
  long P=M;
  long Q=N;
  while (P!=0) {
    long R=Q%P;
    Q=P;
    P=R;
    }
  return abs(Q);
  }

private long abs(long num) {
  return (num<0)?-num:num;
  }

private long floordiv(long num,long den) {
  if ((num<0 && den>0 || num>0 && den<0) && num%den!=0) {
    return num/den-1;
    }
  return num/den;
  }

private long ceildiv(long num,long den) {
  if ((num>0 && den>0 || num<0 && den<0) && num%den!=0) {
    return num/den+1;
    }
  return num/den;
  }

private long sqrt(long num) {
  long num1=0;
  long num2=(long)65536*(long)32768;  /* 2^31 */
  while (num2!=0) {
    if ((num1+num2)*(num1+num2)<=num) {
      num1+=num2;
      }
    num2/=2;
    }
  return num1;
  }

/* Calculate factor1*factor2 mod Mod */
private long MultMod(long factor1, long factor2, long Mod) {
  long aux;
  aux=factor1*factor2-Mod*(long)(((double)factor1*(double)factor2)/(double)Mod);
  if (aux>=Mod) {return aux-Mod;}
  if (aux<0) {return aux+Mod;}
  return aux;
  }

private long ModPow(long Base, long Exp, long Mod) {
  long Pot, Pwr, mask, value;
  if (Exp==0) {return 1L;}
  mask = 1L;
  Pot = 1L;
  Pwr = Base;
  value = 0;
  while (true) {
    if ((Exp & mask) != 0) {
      Pot = MultMod(Pot,Pwr,Mod);
      value += mask;
      if (value == Exp) {return Pot;}
      }
    mask *= 2L;
    Pwr = MultMod(Pwr,Pwr,Mod);
    }
  }

private long ModInv(long Val, long Mod) {
  long U1,U3,V1,V3,Aux,Q;
  U1=1;U3=Val;V1=0;V3=Mod;
  while (V3!=0) {
    Q=U3/V3;
    Aux=U1-V1*Q;U1=V1;V1=Aux;
    Aux=U3-V3*Q;U3=V3;V3=Aux;
    }
  return (U1+Mod)%Mod;
  }

private boolean Mod(long mod) {
  for (long x=0;x<mod;x++) {
    long z=A*x*x+D*x+F;
    long t=B*x+E;
    for (long y=0;y<mod;y++) {
      if ((z+y*(t+C*y))%mod == 0) {
        return false;
        }
      }
    }
  if (teach) {
    w("No solutions found using mod "+mod+",");
    NoSol();
    }
  return true;
  }

private long DivideGcd(long P,long Q,long R,long S,String x1,String y) {
  byte t;
  long N=gcd(P,gcd(Q,R));
  if (N!=1) {
    if (teach) {
      w("We must get the gcd of all terms except the constant:<BR>");
      w("gcd("+P);
      if (Q!=0) {
        w(", "+Q);
        }
      if (R!=0) {
        w(", "+R);
        }
      w(") = "+N+"<P>");
      }
    if (S%N!=0) {
      NoGcd(S);
      return 0;
      }
    if (teach) {
      w(divgcd);
      P/=N;Q/=N;R/=N;S/=N;
      ShowEq(P,0,0,Q,R,S,x1,y);
      w(" = 0<P>");
      }
    }
  if (teach) {
    if (abs(R)!=1) {
      w("This means that ");
      ShowEq(P,0,0,Q,0,S,x1,y);
      w(" should be a multiple of "+abs(R)+"<P>");
      w("To determine this, we should try all values of "+x1+" from 0 to "+(abs(R)-1)+" to check if the condition holds.<BR>");
      t=0;
      for (long u=0;u<abs(R);u++) {
        if ((P*u*u+Q*u+S)%R == 0) {
          if (t!=0) {
            w(", "+u);
            }
          else {
            w("The values of "+x1+" (mod "+abs(R)+") are: "+u);
            t=1;
            }
          }
        }
      if (t==0) {
        w("The modular equation is not satisfied by any "+x1);
        NoSol();
        return 0;
        }
      w("<UL>");
      }
    else {
      t=1;
      }
    }
  return N;
  }

private String par(long num) {
  if (num<0) {
    return "("+num+")";
    }
  return ""+num;
  }

private String par1(long num) {
  return (num==1)?"":par(num);
  }

private byte Linear(long D,long E,long F) {
  long Tx;
  int t;
  if (teach) {
    w("This is a linear equation ");
    ShowLin(D,E,F,"x","y");
    w(" = 0<P>");
    }
  if (D==0) {
    if (E==0) {
      if (F!=0) {
        return 1;             // No solutions
        } 
      else {
        return 2;             // Infinite number of solutions
        }
      }
    if (F%E!=0) {
      return 1;               // No solutions
      } 
    else {
      Xi=0;Xl=1;Yi=-F/E;Yl=0;
      return 0;               // Solution found
      }
    }
  if (E==0) {
    if (F%D!=0) {
      return 1;               // No solutions
      } 
    else {
      Xi=-F/D;Xl=0;Yi=0;Yl=1;
      return 0;               // Solution found
      }
    }
  long Q=gcd(D,E);
  if (Q!=1 && Q!=-1) {
    if (teach) {
      w("To solve it, we first find the <B>gcd</B> of the linear coefficients, that is: <B>gcd</B>("+D+", "+E+") = "+Q+".<P>");
      }
    if (F%Q!=0) {
      NoGcd(F);
      return 1;               // No solutions
      }
    D=D/Q;E=E/Q;F=F/Q;
    }
  if (teach) {
    w(divgcd);
    ShowLin(D,E,F,"x","y");
    w(" = 0<P>");
    w("Now we must apply the Generalized Euclidean algorithm:<P>");
    }
  long U1=1;long U2=0;long U3=D;
  long V1=0;long V2=1;long V3=E;
  t=1;
  while (V3!=0) {
    if (teach) {
      w("Step "+t+": "+par(U1)+" * "+par(D)+" + "+par(U2)+" * "+par(E)+" = "+U3+"<BR>");
      }
    long q=floordiv(U3,V3);
    long T1=U1-q*V1;long T2=U2-q*V2;long T3=U3-q*V3;
    U1=V1;U2=V2;U3=V3;
    V1=T1;V2=T2;V3=T3;
    t++;
    }
  Xi=-U1*F/U3;Xl=E;Yi=-U2*F/U3;Yl=-D;
  if (teach) {
    w("Step "+t+": "+par(U1)+" * "+par(D)+" + "+par(U2)+" * "+par(E)+" = "+U3+"<BR>");
    w("<BR>Multiplying the last equation by "+par(-F/U3)+" we obtain:<BR>");
    w(par(Xi)+" * "+par(D)+" + "+par(Yi)+" * "+par(E)+" = "+(-F)+"<P>");
    w("Adding and subtracting "+par(D)+" * "+par(E)+" t' we obtain:<BR>");
    w("("+Xi+" + "+par(E)+" t') * "+par(D)+" + "+Yi+" - "+par(D)+" t') * "+par(E)+" = "+(-F)+"<P>");
    w("So, the solution is given by the set:<BR>");
    PrintLinear((byte)0,"t'");
    }
  V1 = D*D + E*E;
  Tx = floordiv((D*Yi-E*Xi)+V1/2,V1);
  if (teach) {
    w("By making the substitution t = "+Tx+" + t' we finally obtain:<P>");
    }
  Xi+=E*Tx;Yi+=-D*Tx;
  if (Xl<0 && Yl<0) {
    Xl=-Xl;Yl=-Yl;
    }
  return 0;
  }

private void showAlso() {
  if (also) {
    w("and also:<BR>");
    } 
  else {
    also=true;
    }
  }

private byte PrintLinear(byte Ret, String va) {
  if (Ret==1) {
    return 0;
    }
  if (va.equals("t")) {
    showAlso();
    }
  if (Ret==2) {
    w("x, y: any integer");
    return 1;
    }
  if (ExchXY) {
    long T=Xi;Xi=Yi;Yi=T;T=Xl;Xl=Yl;Yl=T;
    }
  w("x = ");
  if (Xi==0 && Xl==0) {
    w("0");
    }
  if (Xi!=0) {
    w(""+Xi);
    }
  if (Xl<0) {
    w(" - ");
    }
  if (Xl>0 && Xi!=0) {
    w(" + ");
    }
  if (Xl!=0) {
    if (abs(Xl)==1) {
      w(" "+va);
      }
    else {
      w(abs(Xl)+" "+va);
      }
    }
  w("<BR>y = ");
  if (Yi==0 && Yl==0) {
    w("0");
    }
  if (Yi!=0) {
    w(""+Yi);
    }
  if (Yl<0) {
    w(" - ");
    }
  if (Yl>0 && Yi!=0) {
    w(" + ");
    }
  if (Yl!=0) {
    if (abs(Yl)==1) {
      w(" "+va);
      }
    else {
      w(abs(Yl)+" "+va);
      }
    }
  w("<P>");
  return 0;
  }

private void PrintQuad(long A,long B,long C) {
  if (A==1) {
    w(" t"+sq);
    }
  else {
    if (A==-1) {
      w("-t"+sq);
      }
    else {
      if (A!=0) {
        w(A+" t"+sq);
        }
      }
    }
  if (B<0) {
    w(" - ");
    }
  if (B>0 && A!=0) {
    w(" + ");
    }
  if (abs(B)==1) {
    w("t ");
    }
  else {
    if (B!=0) {
      w(abs(B)+" t");
      }
    }
  if (C<0) {
    w(" - ");
    }
  if (C>0 && (B!=0 || A!=0)) {
    w(" + ");
    }
  if (C!=0) {
    w(""+abs(C));
    }
  w("<BR>");
  }

private void ShowXY(long X,long Y) {
  showAlso();
  if (ExchXY==false) {
    w("x = "+X);
    w("<BR>y = "+Y+"<P>");
    }
  else {
    w("x = "+Y);
    w("<BR>y = "+X+"<P>");
    }
  }

private void ShowX1Y1(long X1,long Y1,long A,long B,long D,long D1,long E1) {
  if ((Y1-E1)%D1==0) {
    long Y=(Y1-E1)/D1;
    if ((X1-B*Y-D)%(2*A)==0) {
      long X=(X1-B*Y-D)/(2*A);
      ShowXY(X,Y);
      }
    if (X1!=0 && (-X1-B*Y-D)%(2*A)==0) {
      long X=(-X1-B*Y-D)/(2*A);
      ShowXY(X,Y);
      }
    }
  if (Y1!=0 && (-Y1-E1)%D1==0) {
    long Y=(-Y1-E1)/D1;
    if ((X1-B*Y-D)%(2*A)==0) {
      long X=(X1-B*Y-D)/(2*A);
      ShowXY(X,Y);
      }
    if (X1!=0 && (-X1-B*Y-D)%(2*A)==0) {
      long X=(-X1-B*Y-D)/(2*A);
      ShowXY(X,Y);
      }
    }
  }

private byte Show(long num, String str, byte t) {
  if (t==2) {
    txt="";
    }
  if (num!=0) {
    String str1="";
    if ((t&1)!=0 && num>0) {
      str1+=" +";
      }
    if (num<0) {
      str1+=" -";
      }
    if (abs(num)!=1) {
      str1+=" "+abs(num);
      }
    str1+=str;
    if ((t&2)!=0) {
      txt+=str1;
      } 
    else {
      w(str1);
      }
    return (byte)(t|1);
    }
  return t;
  }

private void Show1(long num, byte t) {
  byte u=Show(num,"",t);
  if ((u&1)==0 || abs(num)==1) {
    if ((u&2)!=0) {
      txt+=abs(num);
      } 
    else {
      w(""+abs(num));
      }
    }
  }

private void ShowElipSol(long A,long B,long D,long u,String x,String x1,String y,long w2) {
  if (teach) {
    w("<LI>");
    ShowLin(2*A,B,D,x,y);
    w(" = "+w2+"<BR>");
    ShowLin(2*A,B,D,x,par(u));
    w(" = "+w2+"<BR>");
    w(par1(2*A)+x+" = "+(w2-D-B*u)+"<BR>");
    if ((w2-D-B*u)%(2*A)==0) {
      w(x+" = "+((w2-D-B*u)/(2*A)));
      }
    else {
      w("There is no integer solution for "+x+"<P>");
      }
    }
  if ((w2-D-B*u)%(2*A)==0) {
    if (teach) {
      w("<TABLE BORDER=1><TR><TH>");
      }
    ShowXY((w2-D-B*u)/(2*A),u);
    if (teach) {
      w("</TABLE><P>");
      } 
    else {
      w("<P>");
      }
    also=true;
    }
  }

private void SolByFact(long R,long T,long B,long D,long E) {
  if (teach) {
    w("<LI>"+T+" is a factor of "+R+", so we can set:<BR>");
    Show1(E,Show(B," x",(byte)0));
    w(" = "+T);
    if (E!=0) {
      w("<BR>This means that ");
      Show(B," x = ",(byte)0);
      w(""+(T-E));
      }
    w("<P>");
    Show1(D,Show(B," y",(byte)0));
    w(" = "+R+"/"+par(T)+" = "+R/T);
    if (D!=0) {
      w("<BR>This means that ");
      Show(B," y = ",(byte)0);
      w(""+(R/T-D));
      }
    w("<P>");
    }
  if ((T-E)%B==0 && (R/T-D)%B==0) {
    if (teach) {
      w("<TABLE BORDER=1><TR><TH>");
      }
    ShowXY((T-E)/B,(R/T-D)/B);
    if (teach) {
      w("</TABLE>");
      }
    }
  else {
    if (teach) {
      w("These equations are not valid in integers.<P>");
      }
    }
  }

private void ShowLin(long D,long E,long F,String x,String y) {
  byte t=Show(D," "+x,(byte)0);
  t=Show(E," "+y,t);
  Show1(F,t);
  }

private void ShowEq(long A,long B,long C,long D,long E,long F,String x,String y) {
  byte t=Show(A," "+x+sq,(byte)0);
  t=Show(B," "+x+y,t);
  t=Show(C," "+y+sq,t);
  t=Show(D," "+x,t);
  t=Show(E," "+y,t);
  Show1(F,t);
}

public void init() {
  H1=new long[20000];
  H2=new long[20000];
  K1=new long[20000];
  K2=new long[20000];
  L1=new long[20000];
  L2=new long[20000];
  teach=false;
  }

private BigInteger value(String nbr) {
  String InputField;
  int ExpressionRC;
  BigInteger [] ExpressionResult = new BigInteger[1];
  InputField = nbr.trim();
  if (InputField.equals("")) {
    return BigInteger.valueOf(0L);
    }
  try {
    ExpressionRC = expression.ComputeExpression(InputField, 1, ExpressionResult);
    } catch (OutOfMemoryError e) {
      textError = "Out of memory.";
      return BigInteger.valueOf(0L);
    }
    switch (ExpressionRC) {
      case -2:
        textError = "Number too high (more than 1000 digits).";
        break;
      case -3:
        textError = "Intermediate expression too high (more than 2000 digits).";
        break;
      case -4:
        textError = "Non-integer division.";
        break;
      case -5:
        textError = "Parentheses mismatch.";
        break;
      case -6:
        textError = "Syntax error.";
        break;
      case -7:
        textError = "Too many parentheses.";
        break;
      case -8:
        textError = "Invalid parameter.";
        break;
      default:
        textError = "";
    }
    return ExpressionResult[0];
  }

public String startCalc(String valA, String valB, String valC, String valD, String valE, String valF, String valDig, String teachMode, String superMode) {
  BigInteger ValDig;

  if (calcThread != null) {
//    TerminateThread = true;
    try {
      calcThread.join();        /* Wait until the solving thread dies */
      } catch (InterruptedException ie) {};
    }
  calcThread = new Thread(this);  /* Start solving thread */
  try {
    ValA = value(valA);
    if (textError.length() > 0) {return "Coefficient x<SUP>2</SUP>: "+textError;}
    ValB = value(valB);
    if (textError.length() > 0) {return "Coefficient xy: "+textError;}
    ValC = value(valC);
    if (textError.length() > 0) {return "Coefficient y<SUP>2</SUP>: "+textError;}
    ValD = value(valD);
    if (textError.length() > 0) {return "Coefficient x: "+textError;}
    ValE = value(valE);
    if (textError.length() > 0) {return "Coefficient y: "+textError;}
    ValF = value(valF);
    if (textError.length() > 0) {return "Constant term: "+textError;}
    ValDig = value(valDig);
    if (textError.length() > 0) {return "Number of digits: "+textError;}
    } catch (Exception e) {return "Invalid data entered";}
  teach = teachMode.equals("step");
  sq = (superMode.equals("yes")?"<SUP>2</SUP>":"^2");
  digitsInGroup = ValDig.intValue();
  if (digitsInGroup < 0) {
    return "Number of digits: It should be positive.";
    }
  calcThread.start();
  return "";
  }

public String resultCalc() {
  if (calcThread == null) {
    return info.toString();
    }
  return "";
  }

public void run() {
  info=new StringBuffer();
  Date OldDate=new Date();
  long Old=OldDate.getTime();
  NbrSols=0;
  try {
    SolveEquation();
    }
  catch (ArithmeticException e) {
    w("<BR>exception: " + e.getMessage());
    e.printStackTrace();
    }
  if (also==false) {
    w(msg);
    }
  Date NewDate=new Date();
  long New=NewDate.getTime();
  w("<P>Calculation time: ");
  int t=(int)(((New-Old)/1000+86400)%86400);
  w(t/3600+"h "+((t%3600)/60)+"m "+(t%60)+"s");
  w("<P>If you found any mistake or any solution is missing, please send me an <A HREF=\"mailto:alpertron@hotmail.com?subject=Error solving: "+txt.replace(' ',' ')+"\">e-mail</A><P>");
  w("<A HREF=\"GBOOK.HTM\">Click here</A> to view or sign my Guestbook.<P>");
  calcThread = null;
  }

private void ChangeSign(long Nbr[],byte type) {
  long Cy=Bi;
  if (LargeNumberIsZero(Nbr)) {return;}
  for (int s=1;s<=abs(Nbr[0]);s++) {
    Nbr[s]=Cy-Nbr[s];
    if (Nbr[s]==Bi) {
      Nbr[s]=0;
      }
    else {
      Cy=Bi-1;
      }
    }
  if (type==1) {Nbr[0]=-Nbr[0];}
  }

private void ShowLargeNumber(long Nbr[]) {
  String value;
  int i;

  if (Nbr[0]<0) {
    w("-");
    ChangeSign(Nbr,(byte)0);
    }
  int s=(int)abs(Nbr[0]);
  while(s>1 && Nbr[s]==0) {
    s--;
    }
  int Dig=0;
  for (long R=Nbr[s];R>=1;R/=10) {
    Dig++;
    }
  value = Nbr[s]+"";
  while (s>1) {
    Dig += 18;
    value += (Bi+Nbr[--s]+"").substring(1);
    }
  i = (value.length()-1)%digitsInGroup + 1;
  w(value.substring(0,i));
  while (i<value.length()) {
    w(" "+value.substring(i,i+digitsInGroup));
    i += digitsInGroup;
    }
  if (Dig>6) {
    w(" <FONT SIZE=\"-1\">("+Dig+" digits)</FONT>");
    }
  if (Nbr[0]<0) {
    ChangeSign(Nbr,(byte)0);
    }
  }

private void ExtendSign(long Nbr[]){
  if(Nbr[0]<0) {
    Nbr[1-(int)Nbr[0]]=Nbr[2-(int)Nbr[0]]=Bi-1;
    }
  else {
    Nbr[1+(int)Nbr[0]]=Nbr[2+(int)Nbr[0]]=0;
    }
  }

private void AdjustSign(long Nbr[]) {
  int S=(int)abs(Nbr[0])+2;
  if(Nbr[S]*2>=Bi) {
    while(Nbr[S]==Bi-1 && --S>abs(Nbr[0])){
      }
    Nbr[0]=-S;
    }
  else {
    while(Nbr[S]==0 && --S>abs(Nbr[0])){
      }
    Nbr[0]=S;
    }
  }

private boolean LargeNumberIsZero(long Nbr[]) {
  int i;
  if (Nbr[0]>0) {
    for (i = (int)abs(Nbr[0]); i>0; i--) {
      if (Nbr[i]!=0) {
        return false;
        }
      }
    return true;
    }
  return false;
  }

private long DivLargeNumber(long Nbr[],long Coef,long Dest[]) {
  long Q, R, Rem, C;
  int i;
  if (LargeNumberIsZero(Nbr)) {return 0;};
  C=(Coef < 0? -Coef: Coef);
  Rem=(Nbr[0]<0? C - 1: 0);
  for (i = (int)(Nbr[0]<0? -Nbr[0]: Nbr[0]); i>0; i--) {
    R=Rem*Mi+Nbr[i]/Mi;
    Q=(long)((double)R/(double)C);
    R=R-C*Q;
    if (R>=C && R<C+C) {Q++;R-=C;}
    else {
      if (R>=-C && R<0) {Q--;R+=C;}
      }
    Rem=Q;
    R=R*Mi+Nbr[i]%Mi;
    Q=(long)((double)R/(double)C);
    R=R-C*Q;
    if (R>=C && R<C+C) {Q++;R-=C;}
    else {
      if (R>=-C && R<0) {Q--;R+=C;}
      }
    Dest[i]=Rem*Mi+Q;
    Rem=R;
    }
  Dest[0]=Nbr[0];
  AdjustSign(Dest);
  if (Coef<0) {
    ChangeSign(Dest,(byte)1);
    }
  return Rem;
  }


private void Mult2LargeNumbers(long Nbr1[],long Nbr2[],long Dest[]) {
  long LSM1, MSM1, LSM2, MSM2;
  long Re2, Re1, Re4, Tmp, MaxJ;
  int i,j;
  if (Nbr1[0]<0) {
    ChangeSign(Nbr1,(byte)0);
    }
  if (Nbr1!=Nbr2 && Nbr2[0]<0) {
    ChangeSign(Nbr2,(byte)0);
    }
  i=(int)abs(Nbr1[0]);
  j=(int)abs(Nbr2[0]);
  Nbr1[i+1]=0;
  Nbr1[i+2]=0;
  Nbr2[j+1]=0;
  Nbr2[j+2]=0;
  MaxJ = j+2;
  for (i = i+j+3; i > 0; i--) {
    Dest[i]=0;
    }
  for (i = 1; i <= abs(Nbr1[0])+2; i++) {
    LSM1=Nbr1[i]%Mi;MSM1=Nbr1[i]/Mi;
    long CyL=0;long CyH=0;
    for (j = 1; j <= MaxJ; j++) {
      LSM2=Nbr2[j]%Mi;MSM2=Nbr2[j]/Mi;
      Tmp=LSM1*LSM2+CyL+Dest[i+j-1];Re2=Tmp/Mi;Re1=Tmp%Mi;
      Tmp=MSM1*LSM2+LSM1*MSM2+CyH+Re2;Re4=Tmp/Mi;Re1+=(Tmp%Mi)*Mi;
      Tmp=MSM1*MSM2+Re4+Re1/Bi;CyH=Tmp/Mi;CyL=Tmp%Mi;
      Dest[i+j-1]=Re1%Bi;
      }
    }
  Dest[0] = abs(Nbr1[0]) + abs(Nbr2[0])+2;
  if (Nbr1[0]*Nbr2[0]<0) {
    ChangeSign(Dest,(byte)0);
    }
  Dest[0]-=2;
  if (Nbr1!=Nbr2 && Nbr2[0]<0) {
    ChangeSign(Nbr2,(byte)0);
    }
  if (Nbr1[0]<0) {
    ChangeSign(Nbr1,(byte)0);
    }
  AdjustSign(Dest);
  }

private void MultLargeNumber(long Coef,long Nbr[],long Dest[]) {
  long LSM, MSM, Re2, Re1, Re4, Tmp;
  ExtendSign(Nbr);
  long CfH=Coef/Mi;long CfL=Coef%Mi;
  long CyL=0;long CyH=0;
  for(int S=1;S<=abs(Nbr[0])+2;S++){
    LSM=Nbr[S]%Mi;MSM=Nbr[S]/Mi;
    Tmp=CfL*LSM+CyL;Re2=Tmp/Mi;Re1=Tmp%Mi;
    Tmp=CfH*LSM+CfL*MSM+CyH+Re2;Re4=Tmp/Mi;Re1+=(Tmp%Mi)*Mi;
    Tmp=CfH*MSM+Re4-(Re1<0?1:0);CyH=Tmp/Mi;CyL=Tmp%Mi;
    Dest[S]=(Re1+Bi)%Bi;
    }
  Dest[0]=Nbr[0];
  AdjustSign(Dest);
  }

private void AddLargeNumbers(long Nbr1[],long Nbr2[],long Sum[]) {
  int S=(int)abs(Nbr1[0]);
  long T=(Nbr1[0]<0?Bi-1:0);
  long U=(Nbr2[0]<0?Bi-1:0);
  while(abs(Nbr2[0])>S) {
    Nbr1[++S]=T;
    }
  S=(int)abs(Nbr2[0]);
  while(S<abs(Nbr1[0])) {
    Nbr2[++S]=U;
    }
  Nbr1[S+1]=Nbr1[S+2]=T;
  Nbr2[S+1]=Nbr2[S+2]=U;
  Sum[0]=S;
  long Cy=0;
  for(S=1;S<=Sum[0]+2;S++) {
    Cy+=Nbr1[S]+Nbr2[S];Sum[S]=Cy%Bi;Cy/=Bi;
    }
  AdjustSign(Sum);
  }

private void AddLargeLong(long Src[],long Nbr,long Dest[]) {
  ExtendSign(Src);
  long Cy=Nbr;
  for (int i=1;i<abs(Src[0])+3;i++) {
    Dest[i]=Src[i]+(Cy<0?Cy+Bi:Cy);
    Cy=Dest[i]/Bi-(Cy<0?1:0);
    Dest[i]%=Bi;
    }
  Dest[0]=abs(Src[0]);
  AdjustSign(Dest);
  }

private void MultAddLargeNumbers(long CPrev,long Prev[],long CAct,long Act[],long Dest[]) {
  long Tmp, Re1, Re2, Re4, PL, PH, AL, AH;
  int S=(int)abs(Act[0]);
  long T=(Prev[0]<0?Bi-1:0);
  long U=(Act[0]<0?Bi-1:0);
  while(S<abs(Prev[0])) {
    Act[++S]=U;
    }
  S=(int)abs(Prev[0]);
  while(S<abs(Act[0])) {
    Prev[++S]=T;
    }
  Prev[S+1]=Prev[S+2]=T;
  Act[S+1]=Act[S+2]=U;
  Dest[0]=S;
  long CPH=CPrev/Mi;long CPL=CPrev%Mi;
  long CAH=CAct/Mi;long CAL=CAct%Mi;
  long CyL=0;long CyH=0;
  for(S=1;S<=Dest[0]+2;S++) {
    PL=Prev[S]%Mi;PH=Prev[S]/Mi;
    AL=Act[S]%Mi;AH=Act[S]/Mi;
    Tmp=PL*CPL+AL*CAL+CyL;Re2=Tmp/Mi;Re1=Tmp%Mi;
    Tmp=PH*CPL+PL*CPH+AH*CAL+AL*CAH+CyH+Re2;Re4=Tmp/Mi;Re1+=(Tmp%Mi)*Mi;
    Tmp=PH*CPH+AH*CAH+Re4-(Re1<0?1:0);CyH=Tmp/Mi;CyL=Tmp%Mi;
    Dest[S]=(Re1+Bi)%Bi;
    }
  AdjustSign(Dest);
  }

private void LongToDoublePrecLong(long Nbr, long Out[]) {
  Out[0] = Nbr & DosALa32_1;
  Out[1] = Nbr >>> 32;
  Out[2] = Out[3] = Out[4] = Out[5] = (Nbr<0? DosALa32_1: 0);
  }

private void ChSignDoublePrecLong(long Nbr[]) {
  long Cy=0;
  for(int i=0; i<6; i++) {
    Cy -= Nbr[i];
    Nbr[i]=(Cy>=0?Cy:Cy+DosALa32);
    Cy=(Cy>=0?0:-1);
    }
  }

private void AddDoublePrecLong(long Nbr1[], long Nbr2[], long Sum[]) {
  long Cy=0;
  for(int i=0; i<6; i++) {
    Cy += Nbr1[i]+Nbr2[i];
    Sum[i] = Cy & DosALa32_1;
    Cy >>= 32;
    }
  }

private void SubtDoublePrecLong(long Nbr1[], long Nbr2[], long Diff[]) {
  long Cy=0;
  for(int i=0; i<6; i++) {
    Cy += Nbr1[i]-Nbr2[i];
    Diff[i] = Cy & DosALa32_1;
    Cy = (Cy>=0?0:-1);
    }
  }

private void MultDoublePrecLong(long Nbr1[], long Nbr2[], long Prod[]) {
  long Cy, Pr;
  int i,j;
  Cy = Pr = 0;
  for (i=0; i<6; i++) {
    Pr = Cy & DosALa32_1;
    Cy >>>= 32;
    for (j=0; j<=i; j++) {
      Pr += Nbr1[j]*Nbr2[i-j];
      Cy += (Pr >>> 32);
      Pr &= DosALa32_1;
      }
    Prod[i] = Pr;
    }
  }

private long DivDoublePrecLong(long Dividend[], long Divisor[]) {
  boolean ChSignDividend = false;
  boolean ChSignDivisor = false;
  double dDividend, dDivisor;
  long QuotientH, QuotientL, Cy;
  long D0, D1, D2, D3, D4, D5, E0, E1, E2, E3, E4, E5;

  if (Dividend[5] >= DosALa31) {
    ChSignDividend = true;
    ChSignDoublePrecLong(Dividend);
    }
  if (Divisor[5] >= DosALa31) {
    ChSignDivisor = true;
    ChSignDoublePrecLong(Divisor);
    }

  dDividend = ((((((double)Dividend[5]*dDosALa32+
                   (double)Dividend[4])*dDosALa32+
                   (double)Dividend[3])*dDosALa32+
                   (double)Dividend[2])*dDosALa32+
                   (double)Dividend[1])*dDosALa32+
                   (double)Dividend[0]);

  dDivisor  = ((((((double)Divisor[5]*dDosALa32+
                   (double)Divisor[4])*dDosALa32+
                   (double)Divisor[3])*dDosALa32+
                   (double)Divisor[2])*dDosALa32+
                   (double)Divisor[1])*dDosALa32+
                   (double)Divisor[0]);

  QuotientH = (long)(dDividend/dDivisor/dDosALa32)+3;
  do {
    QuotientH--;
    D1 = Dividend[1] - Divisor[0] * QuotientH;
    if (D1>Dividend[1] || D1<0) {
      Cy = (D1>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    D2 = Dividend[2] - Divisor[1] * QuotientH + Cy;
    if (D2>Dividend[2] || D2<0) {
      Cy = (D2>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    D3 = Dividend[3] - Divisor[2] * QuotientH + Cy;
    if (D3>Dividend[3] || D3<0) {
      Cy = (D3>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    D4 = Dividend[4] - Divisor[3] * QuotientH + Cy;
    if (D4>Dividend[4] || D4<0) {
      Cy = (D4>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    D5 = Dividend[5] - Divisor[4] * QuotientH + Cy;
    } while (D5 < 0);

  D0 = Dividend[0];
  D1 &= DosALa32_1;
  D2 &= DosALa32_1;
  D3 &= DosALa32_1;
  D4 &= DosALa32_1;

  dDividend = ((((((double)D5*dDosALa32+
                   (double)D4)*dDosALa32+
                   (double)D3)*dDosALa32+
                   (double)D2)*dDosALa32+
                   (double)D1)*dDosALa32+
                   (double)D0);

  QuotientL = (long)(dDividend/dDivisor)+3;
  do {
    QuotientL--;
    E0 = D0 - Divisor[0] * QuotientL;
    if (E0>D0 || E0<0) {
      Cy = (E0>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E1 = D1 - Divisor[1] * QuotientL + Cy;
    if (E1>D1 || E1<0) {
      Cy = (E1>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E2 = D2 - Divisor[2] * QuotientL + Cy;
    if (E2>D2 || E2<0) {
      Cy = (E2>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E3 = D3 - Divisor[3] * QuotientL + Cy;
    if (E3>D3 || E3<0) {
      Cy = (E3>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E4 = D4 - Divisor[4] * QuotientL + Cy;
    if (E4>D4 || E4<0) {
      Cy = (E4>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E5 = D5 - Divisor[5] * QuotientL + Cy;
    } while (E5 < 0);

  if (ChSignDividend != ChSignDivisor) {
    if (((E0 | E1 | E2 | E3 | E4 | E5) & DosALa32_1) != 0) {
      QuotientL++;
      }
    }

  if (ChSignDividend) {
    ChSignDoublePrecLong(Dividend);
    }
  if (ChSignDivisor) {
    ChSignDoublePrecLong(Divisor);
    }

  if (ChSignDividend == ChSignDivisor) {
    return QuotientH << 32 | QuotientL;
    }
  return -(QuotientH << 32 | QuotientL);
  }

private void DivideDoublePrecLong(long Dividend[], long Divisor[], long Quotient[]) {
  boolean ChSignDividend = false;
  boolean ChSignDivisor = false;
  double dDividend, dDivisor, dAux;
  double dDosALa32 = (double)DosALa32;
  long QuotientH, D0, D1, D2, D3, D4, D5, E0, E1, E2, E3, E4, E5, Cy;

  if (Dividend[5] >= DosALa31) {
    ChSignDividend = true;
    ChSignDoublePrecLong(Dividend);
    }
  if (Divisor[5] >= DosALa31) {
    ChSignDivisor = true;
    ChSignDoublePrecLong(Divisor);
    }

  Quotient[0]=Quotient[1]=Quotient[2]=Quotient[3]=Quotient[4]=Quotient[5]=0;

  D0 = Dividend[0];
  D1 = Dividend[1];
  D2 = Dividend[2];
  D3 = Dividend[3];
  D4 = Dividend[4];
  D5 = Dividend[5];

  dDividend = ((((((double)D5*dDosALa32+
                   (double)D4)*dDosALa32+
                   (double)D3)*dDosALa32+
                   (double)D2)*dDosALa32+
                   (double)D1)*dDosALa32+
                   (double)D0);

  dDivisor  = ((((((double)Divisor[5]*dDosALa32+
                   (double)Divisor[4])*dDosALa32+
                   (double)Divisor[3])*dDosALa32+
                   (double)Divisor[2])*dDosALa32+
                   (double)Divisor[1])*dDosALa32+
                   (double)Divisor[0]);

  if (dDividend >= dDosALa32 * dDivisor) {
    if (dDividend >= dDosALa64 * dDivisor) {
      if (dDividend >= dDosALa96 * dDivisor) {
        if (dDividend >= dDosALa128 * dDivisor) {
          if (dDividend >= dDosALa160 * dDivisor) {
            QuotientH = (long)(dDividend/dDivisor/dDosALa160)+3;
            do {
              QuotientH--;
              E5 = D5 - Divisor[0] * QuotientH;
              } while (E5 < 0);
            Quotient[5] = QuotientH;
            D5 = E5;

            dDividend = ((((((double)D5*dDosALa32+
                             (double)D4)*dDosALa32+
                             (double)D3)*dDosALa32+
                             (double)D2)*dDosALa32+
                             (double)D1)*dDosALa32+
                             (double)D0);
            }
          QuotientH = (long)(dDividend/dDivisor/dDosALa128)+3;
          do {
            QuotientH--;
            E4 = D4 - Divisor[0] * QuotientH;
            if (E4>D4 || E4<0) {
              Cy = (E4>>>32) - DosALa32;
              }
            else {
              Cy = 0;
              }
            E5 = D5 - Divisor[1] * QuotientH + Cy;
            } while (E5 < 0);
          Quotient[4] = QuotientH;
          D4=E4 & DosALa32_1;
          D5=E5;
          dDividend = ((((((double)D5*dDosALa32+
                           (double)D4)*dDosALa32+
                           (double)D3)*dDosALa32+
                           (double)D2)*dDosALa32+
                           (double)D1)*dDosALa32+
                           (double)D0);
          }
        QuotientH = (long)(dDividend/dDivisor/dDosALa96)+3;
        do {
          QuotientH--;
          E3 = D3 - Divisor[0] * QuotientH;
          if (E3>D3 || E3<0) {
            Cy = (E3>>>32) - DosALa32;
            }
          else {
            Cy = 0;
            }
          E4 = D4 - Divisor[1] * QuotientH + Cy;
          if (E4>D4 || E4<0) {
            Cy = (E4>>>32) - DosALa32;
            }
          else {
            Cy = 0;
            }
          E5 = D5 - Divisor[2] * QuotientH + Cy;
          } while (E5 < 0);
        Quotient[3] = QuotientH;
        D3 = E3 & DosALa32_1;
        D4 = E4 & DosALa32_1;
        D5 = E5;
        dDividend = ((((((double)D5*dDosALa32+
                         (double)D4)*dDosALa32+
                         (double)D3)*dDosALa32+
                         (double)D2)*dDosALa32+
                         (double)D1)*dDosALa32+
                         (double)D0);
        }
      QuotientH = (long)(dDividend/dDivisor/dDosALa64)+3;
      do {
        QuotientH--;
        E2 = D2 - Divisor[0] * QuotientH;
        if (E2>D2 || E2<0) {
          Cy = (E2>>>32) - DosALa32;
          }
        else {
          Cy = 0;
          }
        E3 = D3 - Divisor[1] * QuotientH + Cy;
        if (E3>D3 || E3<0) {
          Cy = (E3>>>32) - DosALa32;
          }
        else {
          Cy = 0;
          }
        E4 = D4 - Divisor[2] * QuotientH + Cy;
        if (E4>D4 || E4<0) {
          Cy = (E4>>>32) - DosALa32;
          }
        else {
          Cy = 0;
          }
        E5 = D5 - Divisor[3] * QuotientH + Cy;
        } while (E5 < 0);
      Quotient[2] = QuotientH;
      D2 = E2 & DosALa32_1;
      D3 = E3 & DosALa32_1;
      D4 = E4 & DosALa32_1;
      D5 = E5;
      dDividend = ((((((double)D5*dDosALa32+
                       (double)D4)*dDosALa32+
                       (double)D3)*dDosALa32+
                       (double)D2)*dDosALa32+
                       (double)D1)*dDosALa32+
                       (double)D0);
      }
    QuotientH = (long)(dDividend/dDivisor/dDosALa32)+3;
    do {
      QuotientH--;
      E1 = D1 - Divisor[0] * QuotientH;
      if (E1>D1 || E1<0) {
        Cy = (E1>>>32) - DosALa32;
        }
      else {
        Cy = 0;
        }
      E2 = D2 - Divisor[1] * QuotientH + Cy;
      if (E2>D2 || E2<0) {
        Cy = (E2>>>32) - DosALa32;
        }
      else {
        Cy = 0;
        }
      E3 = D3 - Divisor[2] * QuotientH + Cy;
      if (E3>D3 || E3<0) {
        Cy = (E3>>>32) - DosALa32;
        }
      else {
        Cy = 0;
        }
      E4 = D4 - Divisor[3] * QuotientH + Cy;
      if (E4>D4 || E4<0) {
        Cy = (E4>>>32) - DosALa32;
        }
      else {
        Cy = 0;
        }
      E5 = D5 - Divisor[4] * QuotientH + Cy;
      } while (E5 < 0);
    Quotient[1] = QuotientH;
    D1 = E1 & DosALa32_1;
    D2 = E2 & DosALa32_1;
    D3 = E3 & DosALa32_1;
    D4 = E4 & DosALa32_1;
    D5 = E5;
    dDividend = ((((((double)D5*dDosALa32+
                     (double)D4)*dDosALa32+
                     (double)D3)*dDosALa32+
                     (double)D2)*dDosALa32+
                     (double)D1)*dDosALa32+
                     (double)D0);
    }
  QuotientH = (long)(dDividend/dDivisor)+3;
  do {
    QuotientH--;
    E0 = D0 - Divisor[0] * QuotientH;
    if (E0>D0 || E0<0) {
      Cy = (E0>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E1 = D1 - Divisor[1] * QuotientH + Cy;
    if (E1>D1 || E1<0) {
      Cy = (E1>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E2 = D2 - Divisor[2] * QuotientH + Cy;
    if (E2>D2 || E2<0) {
      Cy = (E2>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E3 = D3 - Divisor[3] * QuotientH + Cy;
    if (E3>D3 || E3<0) {
      Cy = (E3>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E4 = D4 - Divisor[4] * QuotientH + Cy;
    if (E4>D4 || E4<0) {
      Cy = (E4>>>32) - DosALa32;
      }
    else {
      Cy = 0;
      }
    E5 = D5 - Divisor[5] * QuotientH + Cy;
    } while (E5 < 0);

  if (ChSignDividend != ChSignDivisor) {
    if (((E0 | E1 | E2 | E3 | E4 | E5) & DosALa32_1) != 0) {
      if (++QuotientH == DosALa32) {
        QuotientH = 0;
        if (++Quotient[1] == DosALa32) {
          Quotient[1] = 0;
          if (++Quotient[2] == DosALa32) {
            Quotient[2] = 0;
            if (++Quotient[3] == DosALa32) {
              Quotient[3] = 0;
              if (++Quotient[4] == DosALa32) {
                Quotient[4] = 0;
                if (++Quotient[5] == DosALa32) {
                  Quotient[5] = 0;
                  }
                }
              }
            }
          }
        }
      }
    }

  Quotient[0] = QuotientH;

  if (ChSignDividend) {
    ChSignDoublePrecLong(Dividend);
    }
  if (ChSignDivisor) {
    ChSignDoublePrecLong(Divisor);
    }
  if (ChSignDividend != ChSignDivisor) {
    ChSignDoublePrecLong(Quotient);
    }
  }

private void DivDoublePrecLongByLong(long Dividend[], long Divisor, long Quotient[]) {
  int i;
  boolean ChSignDivisor=false;
  long Divid, Rem = 0;

  if (Divisor < 0) {
    ChSignDivisor = true;
    Divisor = -Divisor;
    }
  if (Dividend[5] >= DosALa31) {
    Rem = Divisor - 1;
    }
  for (i=5; i>=0; i--) {
    Divid = Dividend[i] + (Rem << 32);
    Rem = Divid % Divisor;
    Quotient[i] = Divid / Divisor;
    }
  if (ChSignDivisor) {
    ChSignDoublePrecLong(Quotient);
    }
  }

// Gcd calculation:
// Step 1: Set k<-0, and then repeatedly set k<-k+1, u<-u/2, v<-v/2
//         zero or more times until u and v are not both even.
// Step 2: If u is odd, set t<-(-v) and go to step 4. Otherwise set t<-u.
// Step 3: Set t<-t/2
// Step 4: If t is even, go back to step 3.
// Step 5: If t>0, set u<-t, otherwise set v<-(-t).
// Step 6: Set t<-u-v. If t!=0, go back to step 3.
// Step 7: The GCD is u*2^k.
private void GcdDoublePrecLong(long Nbr1[], long Nbr2[], long Gcd[]) {
  long u[]=new long[6];
  long v[]=new long[6];
  long t[]=new long[6];
  int k;

  System.arraycopy(Nbr1,0,u,0,6);
  System.arraycopy(Nbr2,0,v,0,6);
  if (u[0]==0 && u[1]==0 && u[2]==0 && u[3]==0 && u[4]==0 && u[5]==0) {
    System.arraycopy(v,0,Gcd,0,6);
    return;
    }
  if (v[0]==0 && v[1]==0 && v[2]==0 && v[3]==0 && v[4]==0 && v[5]==0) {
    System.arraycopy(u,0,Gcd,0,6);
    return;
    }
  if (u[5] >= DosALa31) {
    ChSignDoublePrecLong(u);
    }
  if (v[5] >= DosALa31) {
    ChSignDoublePrecLong(v);
    }
  k=0;
  while ((u[0] & 1)==0 && (v[0] & 1)==0) {   // Step 1
    k++;
    DivDoublePrecLongByLong(u,2,u);
    DivDoublePrecLongByLong(v,2,v);
    }
  if ((u[0] & 1)==1) {                       // Step 2
    System.arraycopy(v,0,t,0,6);
    ChSignDoublePrecLong(t);
    }
  else {
    System.arraycopy(u,0,t,0,6);
    }
  do {
    while ((t[0] & 1)==0) {                  // Step 4
      DivDoublePrecLongByLong(t,2,t);        // Step 3
      }
    if (t[5] < DosALa31) {                   // Step 5
      System.arraycopy(t,0,u,0,6);
      }
    else {
      System.arraycopy(t,0,v,0,6);
      ChSignDoublePrecLong(v);
      }
    SubtDoublePrecLong(u,v,t);               // Step 6

    } while (t[0]!=0 || t[1]!=0 || t[2]!=0 || t[3]!=0 || t[4]!=0 || t[5]!=0);
  System.arraycopy(u,0,Gcd,0,6);             // Step 7
  while (k>0) {
    AddDoublePrecLong(Gcd,Gcd,Gcd);
    k--;
    }
  }

private String DoublePrecLongToString(long Nbr[]) {
  long N0, N1, N2, N3, N4, N5, Rem;
  boolean ChSign = false;
  String nbrOutput="";

  if (Nbr[5] >= DosALa31) {
    ChSignDoublePrecLong(Nbr);
    ChSign = true;
    }
  N0 = Nbr[0];
  N1 = Nbr[1];
  N2 = Nbr[2];
  N3 = Nbr[3];
  N4 = Nbr[4];
  N5 = Nbr[5];
  do {
    Rem = N5 % Mi;
    N5 /= Mi;
    N4 += Rem << 32;
    Rem = N4 % Mi;
    N4 /= Mi;
    N3 += Rem << 32;
    Rem = N3 % Mi;
    N3 /= Mi;
    N2 += Rem << 32;
    Rem = N2 % Mi;
    N2 /= Mi;
    N1 += Rem << 32;
    Rem = N1 % Mi;
    N1 /= Mi;
    N0 += Rem << 32;
    Rem = N0 % Mi;
    N0 /= Mi;
    nbrOutput = ((Rem+Mi)+"").substring(1) + nbrOutput;
    } while (N0 != 0 || N1 != 0 || N2 != 0 || N3 != 0 || N4 != 0 || N5 != 0);
  while (nbrOutput.charAt(0)=='0' && nbrOutput.length() > 1) {
    nbrOutput = nbrOutput.substring(1);
    }
  if (ChSign) {
    ChSignDoublePrecLong(Nbr);
    }
  return nbrOutput;
  }
/*******************************************/
/* NextConv:                               */
/*  BigInteger Tmp = Prev * A1 + Act * B1; */
/*  Act = Prev * A2 + Act * B2;            */
/*  Prev = Tmp;                            */
/*******************************************/

private void NextConv(long Prev[],long Act[]) {
  long CP1H, CP1L, CP2H, CP2L;
  long CA1H, CA1L, CA2H, CA2L;
  long Cy1H, Cy1L, Cy2H, Cy2L;
  long PH, PL, AH, AL;
  long Tmp, Re1, Re2, Re4;
  int S=(int)abs(Act[0]);
  long T=(Prev[0]<0?Bi-1:0);
  long U=(Act[0]<0?Bi-1:0);
  while(S<abs(Prev[0])) {
    Act[++S]=U;
    }
  S=(int)abs(Prev[0]);
  while(S<abs(Act[0])) {
    Prev[++S]=T;
    }
  Prev[S+1]=Prev[S+2]=T;
  Act[S+1]=Act[S+2]=U;
  Act[0]=Prev[0]=S;
  CP1H=A1/Mi;CP1L=A1%Mi;
  CP2H=A2/Mi;CP2L=A2%Mi;
  CA1H=B1/Mi;CA1L=B1%Mi;
  CA2H=B2/Mi;CA2L=B2%Mi;
  Cy1L=Cy1H=Cy2L=Cy2H=0;
  for(S=1;S<=Prev[0]+2;S++) {
    PL=Prev[S]%Mi;PH=Prev[S]/Mi;
    AL=Act[S]%Mi;AH=Act[S]/Mi;
    Tmp=PL*CP1L+AL*CA1L+Cy1L;Re2=Tmp/Mi;Re1=Tmp%Mi;
    Tmp=PH*CP1L+PL*CP1H+AH*CA1L+AL*CA1H+Cy1H+Re2;Re4=Tmp/Mi;Re1+=(Tmp%Mi)*Mi;
    Tmp=PH*CP1H+AH*CA1H+Re4-(Re1<0?1:0);Cy1H=Tmp/Mi;Cy1L=Tmp%Mi;
    Prev[S]=(Re1+Bi)%Bi;
    Tmp=PL*CP2L+AL*CA2L+Cy2L;Re2=Tmp/Mi;Re1=Tmp%Mi;
    Tmp=PH*CP2L+PL*CP2H+AH*CA2L+AL*CA2H+Cy2H+Re2;Re4=Tmp/Mi;Re1+=(Tmp%Mi)*Mi;
    Tmp=PH*CP2H+AH*CA2H+Re4-(Re1<0?1:0);Cy2H=Tmp/Mi;Cy2L=Tmp%Mi;
    Act[S]=(Re1+Bi)%Bi;
    }
  AdjustSign(Prev);
  AdjustSign(Act);
  }

private void ShowStatus() {
  byte t=Show(CX2," x^2",(byte)2);
  t=Show(CXY," xy",t);
  t=Show(CY2," y^2",t);
  t=Show(CX," x",t);
  t=Show(CY," y",t);
  Show1(C1,t);
  txt+=" = 0 (x and y integer)";
  showStatus("Calculate "+txt);
  }

private void ShowBigEq(long biA[], long biB[], long biC[], String x, String y) {
  if (biA[5] >= DosALa31) {w("-");}
  w(DoublePrecLongToString(biA)+" "+x+sq+" ");
  if (biB[0]!=0 || biB[1]!=0 || biB[2]!=0 || biB[3]!=0 || biB[4]!=0 && biB[5]!=0) {
    if (biB[5] < DosALa31) {w("+ ");} else {w("- ");}
    w(DoublePrecLongToString(biB)+" "+x+y+" ");
    }
  if (biC[0]!=0 || biC[1]!=0 || biC[2]!=0 || biC[3]!=0 || biC[4]!=0 && biC[5]!=0) {
    if (biC[5] < DosALa31) {w("+ ");} else {w("- ");}
    w(DoublePrecLongToString(biC)+(y.equals("")?" ":" "+y+sq+" "));
    }
  }

private void GetRoot(long biA[],long biB[],long biC[]) {
  long G, K, L, M, P, S, T, Z;
  long biP[]=new long[6];
  long biM[]=new long[6];
  long biZ[]=new long[6];
  long biG[]=new long[6];
  long biK[]=new long[6];
  boolean DENis1;
  boolean NUMis0;

  G = biA[1]*DosALa32+biA[0];
  K = biB[1]*DosALa32+biB[0];
  L = biC[1]*DosALa32+biC[0];
  DET = K*K-4*G*L;
  System.arraycopy(biB,0,NUM,0,6);
  ChSignDoublePrecLong(NUM);
  NUMis0 = (NUM[0]==0 && NUM[1]==0 && NUM[2]==0 && NUM[3]==0 && NUM[4]==0 && NUM[5]==0);
  AddDoublePrecLong(biA,biA,DEN);
  if (teach) {
    w("We have to find the continued fraction expansion of the roots of <B>");
    ShowBigEq(biA, biB, biC, "t", "");
    w("= 0</B>, that is, ");
    if (NUMis0 == false) {
      w("(sqrt("+DET+") ");
      if (NUM[5] < DosALa31) {w("+ ");} else {w("- ");}
      w(DoublePrecLongToString(NUM)+") / ");
      }
    else {
      w("sqrt("+DET+") /");
      }
    if (DEN[5] >= DosALa31) {
      w("(-"+DoublePrecLongToString(DEN)+")");
      }
    else {
      w(DoublePrecLongToString(DEN));
      }
    w("<P>");
    }

  GcdDoublePrecLong(NUM,DEN,biG);
  MultDoublePrecLong(biG,biG,biZ);
  LongToDoublePrecLong(DET,biG);
  GcdDoublePrecLong(biZ,biG,biG);
  G = biG[0] | biG[1] << 32;
  K=1;T=3;
  while(G%4==0) {
    G/=4;K*=2;
    }
  while(G>=T*T) {
    while(G%(T*T)==0) {
      G/=T*T;K*=T;
      }
    T+=2;
    }
  DET/=K*K;
  DivDoublePrecLongByLong(NUM,K,NUM);
  DivDoublePrecLongByLong(DEN,K,DEN);
  DENis1 = (DEN[0]==1 && DEN[1]==0 && DEN[2]==0 && DEN[3]==0 && DEN[4]==0 && DEN[5]==0);
  if (teach && K!=1) {
    w("Simplifying, ");
    if (DENis1 == false && NUMis0 == false) {
      w("(");
      }
    w("sqrt("+DET+")");
    if (NUMis0 == false) {
      if (NUM[5] < DosALa31) {w(" + ");} else {w(" - ");}
      w(DoublePrecLongToString(NUM));
      }
    if (DENis1 == false) {
      if (NUMis0) {
        w(" / ");
        }
      else {
        w(") / ");
        }
      }
    if (DEN[5] >= DosALa31) {
      w("(-"+DoublePrecLongToString(DEN)+")");
      }
    else {
      if (DENis1 == false) {
        w(DoublePrecLongToString(DEN));
        }
      }
    w("<P>");
    }
  DET *= K*K;
  System.arraycopy(biB,0,NUM,0,6);
  ChSignDoublePrecLong(NUM);
  AddDoublePrecLong(biA,biA,DEN);
  SQD=sqrt(DET);
  if (teach) {
    w("The continued fraction expansion is:<BR>");
    System.arraycopy(DEN,0,biP,0,6);
    LongToDoublePrecLong(SQD+(DEN[5] >= DosALa31?1:0), biK);
    AddDoublePrecLong(biK,NUM,biK);
    Z=DivDoublePrecLong(biK,DEN);
    LongToDoublePrecLong(Z,biM);
    MultDoublePrecLong(biM,DEN,biK);
    SubtDoublePrecLong(biK,NUM,biM);
    w(""+Z);
    String sep="+ //";int cont=-1;
    K=P=L=M=-1;
    while(cont<0 || K!=P || L!=M) {
      w(sep);sep=", ";
      if (cont<0 &&
          biP[0]>0 && biP[0]<=SQD+biM[0] && biP[1]==0 && biP[2]==0 &&
          biP[3]==0 && biP[4]==0 && biP[5]==0 &&
          biM[0]>0 && biM[0]<=SQD && biM[1]==0 && biM[2]==0 &&
          biM[3]==0 && biM[4]==0 && biM[5]==0) {
        w("<B>");K=P=biP[0];L=M=biM[0];cont=0;
        }
      if (cont>=0) {
        P=(DET-M*M)/P;  /* both numerator and denominator are positive */
        Z=(SQD+M)/P;
        M=Z*P-M;
        cont++;
        }
      else {
        LongToDoublePrecLong(DET,biZ);
        MultDoublePrecLong(biM,biM,biG);
        SubtDoublePrecLong(biZ,biG,biG);
        DivideDoublePrecLong(biG,biP,biK);
        System.arraycopy(biK,0,biP,0,6);
        LongToDoublePrecLong(SQD+(biP[5]>=DosALa31?1:0),biZ);
        AddDoublePrecLong(biZ,biM,biK);
        Z = DivDoublePrecLong(biK,biP);
        LongToDoublePrecLong(Z,biG);
        MultDoublePrecLong(biG,biP,biZ);
        SubtDoublePrecLong(biZ,biM,biM);
        }
      w(""+Z);
      }
    w("//</B><BR>where the periodic part is marked in bold");
    if (cont>1) {
      w(" (the period has "+cont+" coefficients)");
      }
    w(".<P>");
    }
  }

/* type = 1: Find convergents */
/* type = 2: Find convergents for x^2 + Bxy + ACy^2 = 1 (recursion) */
/* type = 3: Find convergents for modified equation in homogeneous equation */
/* type = 4: Find convergents for modified equation in complete solution */
/* type = 5: Find convergents for x^2 + Bxy + ACy^2 = 1 (mod B^2-4AC) */
private boolean ContFrac(long biA[],byte type,byte SqrtSign,long s,long T,long MagnifY,long A) {
  long P, Z, M, P1, M1, Tmp, K, L, Mu;
  long biP[] = new long[6];
  long biM[] = new long[6];
  long biZ[] = new long[6];
  long biG[] = new long[6];
  long biMu[] = new long[6];
  long biK[] = new long[6];
  long biL[] = new long[6];
  long biM1[] = new long[6];
  long biP1[] = new long[6];
  long H1ModCY1=1,H2ModCY1=0,K1ModCY1=0,K2ModCY1=1;
  boolean Sols, secondDo = true;
  int i;
  int Conv;
  int Co=-1;
  String U=(type==4?"'":"");
  String X1="X";
  String Y1=(MagnifY==1?"Y":"Y'");
  Sols = true;
  if (biA[0]==1 && biA[1]==0 && biA[2]==0 && biA[3]==0 && biA[4]==0 && biA[5]==0) {
    H1[0]=SqrtSign;K1[0]=1;H1[1]=(SqrtSign+Bi)%Bi;K1[1]=0;
    if (type==1) {
      ShowLargeXY(X1,Y1,H1,K1,true,"","");
      return true;            /* Indicate there are solutions */
      }
    if ((type==3 || type==4) && (DET!=5 || A*F<0)) {
      if (ShowHomoSols(type,H1,K1,s,T,MagnifY,"","")) {
        return true;          /* Indicate there are solutions */
        }
      }
    }
  /* Paso = 1: Quick scan for solutions */
  /* Paso = 2: Show actual solutions */
  for (int Paso=(type==2 || DET==5 && A*F>0 && (type==3 || type==4)? 2: 1);
       Sols && Paso<=2; Paso++) {
    Conv = 0;
    Sols = false;
    System.arraycopy(DEN,0,biP,0,6);
    if (SqrtSign < 0) {
      ChSignDoublePrecLong(biP);
      }
    LongToDoublePrecLong(SQD+(biP[5] >= DosALa31?1:0), biK);
    if (SqrtSign < 0) {
      SubtDoublePrecLong(biK,NUM,biK);
      }
    else {
      AddDoublePrecLong(biK,NUM,biK);
      }
    Z=DivDoublePrecLong(biK,biP);
    LongToDoublePrecLong(Z,biM);
    MultDoublePrecLong(biM,DEN,biK);
    SubtDoublePrecLong(biK,NUM,biM);
    if (SqrtSign < 0) {
      ChSignDoublePrecLong(biM);
      }
    if (type == 4) {
      H2ModCY1 = Z%CY1;
      }
                                       // biM=(NUM+Z*DEN)*SqrtSign;
    if (type == 5) {
      A1=B2=1;A2=Z%T;B1=0;
      }
    else {
      H1[0]=SqrtSign;H2[0]=(Z*SqrtSign<0)?-1:1;K1[0]=1;K2[0]=H1[0];
      H1[1]=(SqrtSign+Bi)%Bi;H2[1]=(Z*SqrtSign+Bi)%Bi;K1[1]=0;K2[1]=H1[1];
      A1=B2=1; A2=B1=0;
      }
    Co=-1;
    for (i=0; i<6; i++) {
      biK[i] = biL[i] = DosALa32_1;
      }
    switch (type) {
      case 1:
        LongToDoublePrecLong(-2*F*SqrtSign,biMu);
        break;
      case 3:
      case 4:
        LongToDoublePrecLong(-2*SqrtSign,biMu);
        break;
      default:
        System.arraycopy(DEN,0,biMu,0,6);
        if (SqrtSign > 0) {
          ChSignDoublePrecLong(biMu);
          }
      }
    do {
      LongToDoublePrecLong(DET,biZ);       // biP1 = (DET-biM*biM)/biP
      MultDoublePrecLong(biM,biM,biG);
      SubtDoublePrecLong(biZ,biG,biG);
      DivideDoublePrecLong(biG,biP,biP1);
      LongToDoublePrecLong(SQD+(biP1[5]>=DosALa31?1:0),biZ);
      AddDoublePrecLong(biZ,biM,biK);
      Z = DivDoublePrecLong(biK,biP1);
      LongToDoublePrecLong(Z,biG);         // biM1 = Z * biP1 - biM
      MultDoublePrecLong(biG,biP1,biZ);
      SubtDoublePrecLong(biZ,biM,biM1);
      if (Co<0 &&
          biP[0]>0 && biP[0]<=SQD+biM[0] && biP[1]==0 && biP[2]==0 &&
          biP[3]==0 && biP[4]==0 && biP[5]==0 &&
          biM[0]>0 && biM[0]<=SQD && biM[1]==0 && biM[2]==0 &&
          biM[3]==0 && biM[4]==0 && biM[5]==0) {
        Co=0;
        System.arraycopy(biP,0,biK,0,6);
        System.arraycopy(biM,0,biL,0,6);
        }
      if (type==1 && biP[0]==biMu[0] && biP[1]==biMu[1] && biP[2]==biMu[2] &&
          biP[3]==biMu[3] && biP[4]==biMu[4] && biP[5]==biMu[5]) { // Solution found
        if (Co%2==0 ||
             biK[0]!=biP1[0] || biK[1]!=biP1[1] || biK[2]!=biP1[2] ||
             biK[3]!=biP1[3] || biK[4]!=biP1[4] || biK[5]!=biP1[5] ||
             biL[0]!=biM1[0] || biL[1]!=biM1[1] || biL[2]!=biM1[2] ||
             biL[3]!=biM1[3] || biL[4]!=biM1[4] || biL[5]!=biM1[5]) {
          if (Paso==2) {
            if (A2!=0) {
              NextConv(H1,H2);NextConv(K1,K2);
              A1=B2=1;A2=B1=0;
              }
            ShowLargeXY(X1,Y1,H2,K2,true,"NUM("+Conv+") = ","DEN("+Conv+") = ");
            }
          Sols = true;
          }
        secondDo = false;
        break;
        }
      if (type==3 || type==4) {
        if (Co==0 && A*F>0 && DET==5) {  /* Solution found */
          if (Paso==1) {
            secondDo = false;
            Sols = true;
            break;
            }
          else {
            NextConv(H1,H2);NextConv(K1,K2);
            A1=B2=1;A2=B1=0;
            ChangeSign(H2,(byte)1);
            AddLargeNumbers(H1,H2,H2);
            ChangeSign(H2,(byte)1);
            ChangeSign(K2,(byte)1);
            AddLargeNumbers(K1,K2,K2);
            ChangeSign(K2,(byte)1);
            if (ShowHomoSols(type,H2,K2,s,T,MagnifY,"NUM("+Conv+") - NUM("+(Conv-1)+") = ","DEN("+Conv+") - DEN("+(Conv-1)+") = ")) {
              secondDo = false;
              Sols = true;
              break;
              }
            AddLargeNumbers(H1,H2,H2);
            AddLargeNumbers(K1,K2,K2);
            }
          }
        if (biP1[0]==biMu[0] && biP1[1]==biMu[1] && biP1[2]==biMu[2] &&
            biP1[3]==biMu[3] && biP1[4]==biMu[4] && biP1[5]==biMu[5]) {  // Solution found
          if (Co%2==0 ||
               biK[0]!=biP1[0] || biK[1]!=biP1[1] || biK[2]!=biP1[2] ||
               biK[3]!=biP1[3] || biK[4]!=biP1[4] || biK[5]!=biP1[5] ||
               biL[0]!=biM1[0] || biL[1]!=biM1[1] || biL[2]!=biM1[2] ||
               biL[3]!=biM1[3] || biL[4]!=biM1[4] || biL[5]!=biM1[5]) {
            if (Paso==2) {
              if (A2!=0) {
                NextConv(H1,H2);NextConv(K1,K2);
                A1=B2=1;A2=B1=0;
                }
              if (ShowHomoSols(type,H2,K2,s,T,MagnifY,"NUM("+Conv+") = ","DEN("+Conv+") = ")) {
                secondDo = false;
                Sols = true;
                break;
                }
              }
            else {
              if (type==4) {
                Tmp = H2ModCY1*T;
                if ((Tmp-CY0)%CY1==0 || (Tmp+CY0)%CY1==0) {
                  secondDo = false;
                  Sols = true;
                  break;
                  }
                }
              else {
                secondDo = false;
                Sols = true;
                break;
                }
              }
            }
          }

        if (Paso==1 && type==4) {
          Tmp = (H1ModCY1+Z*H2ModCY1)%CY1;
          H1ModCY1 = H2ModCY1;
          H2ModCY1 = Tmp;
          Tmp = (K1ModCY1+Z*K2ModCY1)%CY1;
          K1ModCY1 = K2ModCY1;
          K2ModCY1 = Tmp;
          }
        }
      System.arraycopy(biM1,0,biM,0,6);
      System.arraycopy(biP1,0,biP,0,6);
      if (Co==0) {Co=1;}
      if (type==5) {
        Tmp=(A1+Z*A2)%T;
        A1=A2;A2=Tmp;
        Tmp=(B1+Z*B2)%T;
        B1=B2;B2=Tmp;
        }
      ChSignDoublePrecLong(biMu);
      if (Paso==2) {
        if (A2!=0 && Z>(Bi/10-A1)/A2 || B2!=0 && Z>(Bi/10-B1)/B2) {
          NextConv(H1,H2);NextConv(K1,K2);A1=B2=1;A2=B1=0;
          }
        A1+=Z*A2;B1+=Z*B2;
        Tmp=A1;A1=A2;A2=Tmp;Tmp=B1;B1=B2;B2=Tmp;
        }
      Conv++;
      } while (Co<0);

    if (secondDo == false) {
      continue;
      }
    Mu = biMu[0] | biMu[1] << 32;
    L = biL[0] | biL[1] << 32;
    K = biK[0] | biK[1] << 32;
    M = biM[0] | biM[1] << 32;
    P = biP[0] | biP[1] << 32;

    do {
      P1 = (DET-M*M)/P;    /* P & Q should be > 0 (See Knuth Ex 4.5.3-12) */
      Z = (SQD+M)/P1;
      M1=Z*P1-M;
      if (type==1 && P==Mu) {    /* Solution found */
        if (Co%2==0 || K!=P1 || L!=M1) {
          if (Paso==2) {
            if (A2!=0) {
              NextConv(H1,H2);NextConv(K1,K2);
              A1=B2=1;A2=B1=0;
              }
            ShowLargeXY(X1,Y1,H1,K1,true,"NUM("+Conv+") = ","DEN("+Conv+") = ");
            }
          Sols = true;
          }
        break;
        }
      if (type==3 || type==4) {
        if ((Co&1)==0 && A*F>0 && DET==5) {   /* Solution found */
          if (Paso==1) {
            Sols = true;
            break;
            }
          else {
            NextConv(H1,H2);NextConv(K1,K2);
            A1=B2=1;A2=B1=0;
            ChangeSign(H2,(byte)1);
            AddLargeNumbers(H1,H2,H2);
            ChangeSign(H2,(byte)1);
            ChangeSign(K2,(byte)1);
            AddLargeNumbers(K1,K2,K2);
            ChangeSign(K2,(byte)1);
            if (ShowHomoSols(type,H2,K2,s,T,MagnifY,"NUM("+Conv+") - NUM("+(Conv-1)+") = ","DEN("+Conv+") - DEN("+(Conv-1)+") = ")) {
              Sols = true;
              break;
              }
            AddLargeNumbers(H1,H2,H2);
            AddLargeNumbers(K1,K2,K2);
            }
          }   
        if (P1==Mu) {   /* Solution found */
          if (Co%2==0 || K!=P1 || L!=M1) {
            if (Paso==2) {
              if (A2!=0) {
                NextConv(H1,H2);NextConv(K1,K2);
                A1=B2=1;A2=B1=0;
                }
              if (ShowHomoSols(type,H2,K2,s,T,MagnifY,"NUM("+Conv+") = ","DEN("+Conv+") = ")) {
                Sols = true;
                break;
                }
              }
            else {
              if (type==4) {
                Tmp = H2ModCY1*T;
                if ((Tmp-CY0)%CY1==0 || (Tmp+CY0)%CY1==0) {
                  Sols = true;
                  break;
                  }
                }
              else {
                Sols = true;
                break;
                }
              }
            }
          }

        if (Paso==1 && type==4) {
          Tmp = (H1ModCY1+Z*H2ModCY1)%CY1;
          H1ModCY1 = H2ModCY1;
          H2ModCY1 = Tmp;
          Tmp = (K1ModCY1+Z*K2ModCY1)%CY1;
          K1ModCY1 = K2ModCY1;
          K2ModCY1 = Tmp;
          }
        }
      Co++;
      if (Co%5000 == 0) {showStatus("Conv: "+Co+" (Eq "+EqNbr+" of "+NbrEqs+")");}
      if (Co%5000 == 2500) {showStatus(NbrSols+" solution"+(NbrSols==1?"":"s"));}
      M=M1;P=P1;
      if (type==2 && P1==Mu) {
        NextConv(H1,H2);
        NextConv(K1,K2);
        Sols = true;
        break;
        }
      if (type==5) {
        if (P1==Mu) {
          NbrCo = Co;
          Sols = true;
          break;
          }
        else {
          Tmp=(A1+Z*A2)%T;
          A1=A2;A2=Tmp;
          Tmp=(B1+Z*B2)%T;
          B1=B2;B2=Tmp;
          }
        }
      Mu=-Mu;
      if (Paso==2) {
        if (A2!=0 && Z>(Bi/10-A1)/A2 || B2!=0 && Z>(Bi/10-B1)/B2) {
          NextConv(H1,H2);NextConv(K1,K2);A1=B2=1;A2=B1=0;
          }
        A1+=Z*A2;B1+=Z*B2;
        Tmp=A1;A1=A2;A2=Tmp;Tmp=B1;B1=B2;B2=Tmp;
        }
      Conv++;
      } while (NbrCo>0?Co != NbrCo: Co%2!=0 || K!=P || L!=M);
    if (type==5) {break;}
    }                            /* end for */
  return Sols;
  }

private boolean ShowHomoSols(byte type, long H1[], long K1[], long s, long T, long MagnifY, String eqX, String eqY) {
  int i;
  String U=(type==4?"'":"");
  String X1="X";
  String Y1=(MagnifY==1?"Y":"Y'");

  if (teach) {
    ShowLargeXY(Y1,"Z",H1,K1,true,eqX,eqY);
    w("Since "+X1+U+"<SUB>0</SUB> = ");
    if (T>1) {
      w(T+" "+UU+U+"<SUB>0</SUB> = "+T+" (");
      }
    ShowLin(s,-F,0,VU+"<SUB>0</SUB>","Z<SUB>0</SUB>");
    if (T>1) {
      w(")<BR>and "+Y1+"<SUB>0</SUB> = "+T+" V<SUB>0</SUB>");
      }
    w(":<BR>");
    }
  MultAddLargeNumbers(s*T,H1,-F*T,K1,L1);
  MultLargeNumber(T,H1,L2);
  if (type==4) {
    if (teach) {
      ShowLargeXY(X1+U,Y1+U,L1,L2,false,"","");
      }
    for (i=(LargeNumberIsZero(L1) && LargeNumberIsZero(L2)?1:0);i<2;i++) {
      AddLargeLong(L2,-CY0,L2);
      if (teach) {
        w(Y1+"<SUB>0</SUB> = (");
        ShowLin(0,i==0?1:-1,-CY0,"",Y1+"'<SUB>0</SUB>");
        w(")/"+par(CY1)+"<BR>");
        }
      if (DivLargeNumber(L2,CY1,L2) !=0) {
        if (teach) {
          w("It is not an integer number.<P>");
          }
        }
      else {
        if (teach) {
          w(X1+"<SUB>0</SUB> = (");
          ShowLin(i==0?1:-1,-B,-D,X1+"'<SUB>0</SUB>",Y1+"<SUB>0</SUB>");
          w(") / "+par(2*A)+"<BR>");
          }
        AddLargeLong(L1,-D,L1);
        MultAddLargeNumbers(1,L1,-B,L2,L1);
        if (DivLargeNumber(L1,2*A,L1) !=0) {
          if (teach) {
            w("It is not an integer number.<P>");
            }
          }
        else {
          if (teach) {
            if (MagnifY!=1) { 
              ShowLargeXY(X1,Y1,L1,L2,false,"","");
              w("Since <B>Y = "+MagnifY+Y1+"</B>");
              }
            w("<P>");
            }
          MultLargeNumber(MagnifY,L2,L2);
          ShowLargeXY("X","Y",L1,L2,false,"","");
          return true;
          }
        }
      MultAddLargeNumbers(-s*T,H1,F*T,K1,L1);
      MultLargeNumber(-T,H1,L2);
      }
    }
  else {
    if (teach) {
      if (MagnifY!=1) { 
        ShowLargeXY(X1,Y1,L1,L2,false,"","");
        w("Since <B>"+Y1+" = "+MagnifY+"Y</B>");
        }
      w("<P>");
      }
    MultLargeNumber(MagnifY,L2,L2);
    if (teach) {
      ShowLargeXY("X","Y",L1,L2,false,"","");
      ChangeSign(L1,(byte)1);
      ChangeSign(L2,(byte)1);
      ShowLargeXY("X","Y",L1,L2,false,"","");
      ChangeSign(L1,(byte)1);
      ChangeSign(L2,(byte)1);
      }
    return true;
    }
  return false;
  }

private void ShowAllLargeSolutions() {
  int i;

  allSolsFound = true;
  also = false;
  for (i=0; i<sortedSolsY.size(); i++) {
    ShowLargeXY("X","Y",((ArrayLongs)(sortedSolsX.elementAt(i))).element,
                        ((ArrayLongs)(sortedSolsY.elementAt(i))).element,
                        false,"","");
    w("<P>");
    }
  }

private void ShowLargeXY(String x,String y,long H1[],long K1[],boolean sol,String eqX,String eqY) {
  if (x.equals("X") && y.equals("Y")) {
    showStatus(NbrSols+" solutions");
    if (teach == false && allSolsFound == false) {
      also=true;
      if (sol && LargeNumberIsZero(K1)) {
        InsertNewSolution(H1,K1);
        ChangeSign(H1,(byte)1);
        InsertNewSolution(H1,K1);
        ChangeSign(H1,(byte)1);
        }
      else {
        if (sol && K1[0]<0) {
          ChangeSign(H1,(byte)1);
          ChangeSign(K1,(byte)1);
          InsertNewSolution(H1,K1);
          ChangeSign(H1,(byte)1);
          ChangeSign(K1,(byte)1);
          }
        else {
          InsertNewSolution(H1,K1);
          }
        }
      return;
      }
    }
  if (teach && y.equals("Y")) {
    w("<B>");
    }
  if (y.equals("Y")) {
    showAlso();
    }
  w(x+"<SUB>0</SUB> = ");
  if (teach) {w(eqX);}
  ShowLargeNumber(H1);
  w("<BR>"+y+"<SUB>0</SUB> = ");
  if (teach) {w(eqY);}
  ShowLargeNumber(K1);
  w("<BR>");
  if (y.equals("Y") && sol) {
    ChangeSign(H1,(byte)1);
    ChangeSign(K1,(byte)1);
    w("<BR>and also:<BR>"+x+"<SUB>0</SUB> = ");
    if (teach) {w(eqX.equals("")?"":"-"+eqX);}
    ShowLargeNumber(H1);
    w("<BR>"+y+"<SUB>0</SUB> = ");
    if (teach) {w(eqY.equals("")?"":"-"+eqY);}
    ShowLargeNumber(K1);w("<BR>");
    ChangeSign(H1,(byte)1);
    ChangeSign(K1,(byte)1);
    }
  if (teach && y.equals("Y")) {
    w("</B>");
    }
  }

     /* Insert new number into sorted solutions vector */
     /* performing a binary search and then using the  */
     /* insertElementAt method of the vector           */
private void InsertNewSolution(long H1[],long K1[]) {
  int sizeVector, indexVector, increment, compare, i;

  sizeVector = sortedSolsY.size();
  indexVector = 0;
  if (sizeVector > 0) {
    increment = 1;
    while (increment * 2 <= sizeVector) {
      increment *= 2;
      }
    while (increment > 0) {  /* Perform binary search */
      if (indexVector + increment <= sizeVector) {
        compare = Compare(((ArrayLongs)(sortedSolsY.elementAt(indexVector + increment-1))).element, K1);
        if (compare == 0) {
          compare = Compare(((ArrayLongs)(sortedSolsX.elementAt(indexVector + increment-1))).element, H1);
          }
        if (compare == 0) {
          return;
          }
        if (compare < 0) {
          indexVector += increment;
          }
        }
      increment /= 2;
      }
    }

  NbrSols++;
  sortedSolsX.insertElementAt(new ArrayLongs((int)(abs(H1[0])+1)), indexVector);
  System.arraycopy(H1,
                   0,
                   ((ArrayLongs)(sortedSolsX.elementAt(indexVector))).element,
                   0,
                   (int)(abs(H1[0])+1));
  sortedSolsY.insertElementAt(new ArrayLongs((int)(abs(K1[0])+1)), indexVector);
  System.arraycopy(K1,
                   0,
                   ((ArrayLongs)(sortedSolsY.elementAt(indexVector))).element,
                   0,
                   (int)(abs(K1[0])+1));
  }

private int Compare(long longarray[], long K1[]) {
  int length, i, retcode;
  boolean signlongarray=false, signK1=false;

  if (longarray[0] < 0) {
    signlongarray = true;
    ChangeSign(longarray,(byte)1);
    }
  if (K1[0] < 0) {
    signK1 = true;
    ChangeSign(K1,(byte)1);
    }

  length = (int)longarray[0];
  if (length < (int)K1[0]) {
    retcode = -1;
    }
  else {
    if (length > (int)K1[0]) {
      retcode = 1;
      }
    else {
      for (i=length; i>0; i--) {
        if (longarray[i] != K1[i]) {
          break;
          }
        }
      if (i==0) {    /* Absolute values are equal */
        if (signlongarray == signK1) {
          retcode = 0;
          }
        else {
          if (signlongarray == false) {
            retcode = 1;
            }
          else {
            retcode = -1;
            }
          }
        }
      else {           /* Different absolute values */
        if (longarray[i] < K1[i]) {
          retcode = -1;
          }
        else {
          retcode = 1;
          }
        }
      }
    }
  if (signlongarray == true) {
    ChangeSign(longarray,(byte)1);
    }
  if (signK1 == true) {
    ChangeSign(K1,(byte)1);
    }
  return retcode;
  }

/*******************************************************/
/* H = Constant term                                   */
/* T = Divisor of the square part of the constant term */
/* A = X^2 coefficient                                 */
/* B = XY coefficient                                  */
/* C = Y^2 coefficient                                 */
/* S = Nothing or apostrophe (complete quad equation)  */
/*******************************************************/
private void SolContFrac(long H,long T,long A,long B,long C,String S) {
  long factor[] = new long[64];
  long P[]=new long[64];
  long Q[]=new long[64];
  long Dif[]=new long[64];   /* Holds difference */
  long mod[]=new long[64];
  long pos[]=new long[64];
  long Tmp, q, r, s, t, u, v, Pp, dif, Sol1, Sol2, Modulo;
  long Tmp1=SQD;
  long Tmp2[] = new long[6];
  long Tmp3[] = new long[6];
  long Tmp4=DET;
  long ValA,ValB,ValC,ValF,ValAM,ValBM,ValCM;
  long VarD,VarK,VarQ,VarR,VarT,VarV,VarW,VarX,VarY,VarY1;
  long biA[] = new long[6];
  long biB[] = new long[6];
  long biC[] = new long[6];
  long biR[] = new long[6];
  long biS[] = new long[6];
  long biT[] = new long[6];
  ArrayLongs LongArrayX, LongArrayY;
  int index,index2,cont;
  int NbrFactors;
  long gcdAF, MagnifY;
  int cuenta=0;
  long OrigA, OrigC;
  boolean ShowHR = false;

  System.arraycopy(NUM,0,Tmp2,0,6);
  System.arraycopy(DEN,0,Tmp3,0,6);
  F=H/T/T;
  if (teach && T>1) {
    w("Since "+T+" * "+T+" is a divisor of the constant term (<B>"+H+"</B>), the solutions should be "+T+" times the solutions of ");
    ShowEq(A,B,C,0,0,F,"u","v");
    w(" = 0.");
    if (abs(F)!=1) {
      w(" Let <B>F</B> be the constant term.");
      }
    w("<P>");
    UU="U";VU="V";UL="u";VL="v";FP="F";
    }
  if (teach && T==1) {
    UU="X"+S;VU="Y"+S;UL="x"+S;VL="y"+S;FP="f"+S;
    }
  gcdAF = gcd(A,F);
  OrigA = A;
  OrigC = C;
  if (teach && gcdAF > 1) {
    w("Since <B>gcd(A,F) = gcd("+A+","+F+") = "+gcdAF+" > 1</B>, we have to replace <B>y = ny'</B> where <B>n</B> is a divisor of <B>gcd(A,F)</B>.<P>");
    }
  for (MagnifY=1; MagnifY*MagnifY <= gcdAF; MagnifY++) {
   do {
    if (gcdAF/MagnifY*MagnifY != gcdAF) {continue;}
    if (teach) {
      if (ShowHR) {w("<HR>");} else {ShowHR = true;}
      }
    MagnifY = gcdAF/MagnifY;
    F=H/T/T/MagnifY;
    ValF=abs(F);
    A = OrigA/MagnifY;
    C = OrigC*MagnifY;
    ValA=(A+ValF)%ValF;
    ValB=(B+ValF)%ValF;
    ValC=(C+ValF)%ValF;
    if (teach) {
      if (MagnifY != 1) {
        w("Let <B>y = "+MagnifY+"y'</B><P>We obtain: ");
        }
      w("<B>");
      ShowEq(A,B,C,0,0,F,"x","y'");
      w(" = 0</B><P>");
      }
           /* Find factors of F */
    NbrFactors = 0;
    Tmp = ValF;
    if (Tmp==1) {
      factor[NbrFactors++] = 1;
      }
    else {
      while ((Tmp%2)==0) {
        factor[NbrFactors++] = 2;
        Tmp/=2;
        }
      while ((Tmp%3)==0) {
        factor[NbrFactors++] = 3;
        Tmp/=3;
        }
      s=5;        /* Sequence of divisors 5, 7, 11, 13, 17, 19,... */
      do {
        while ((Tmp%s)==0) {
          factor[NbrFactors++] = s;
          Tmp/=s;
          }
        s+=2;
        while ((Tmp%s)==0) {
          factor[NbrFactors++] = s;
          Tmp/=s;
          }
        s+=4;
        } while (s*s<=Tmp);
      if (Tmp != 1) {
        factor[NbrFactors++] = Tmp;
        }
      }
    mod[NbrFactors] = Tmp = 1;
    Pp = (2*ValA)%ValF;
    for (index=NbrFactors-1;index>=0;index--) {
      P[index] = Pp;
      Tmp *= factor[index];
      mod[index] = Tmp;
      Pp = MultMod(MultMod(Pp,factor[index],ValF),factor[index],ValF);
      }
    Modulo = factor[NbrFactors-1];
    ValAM=(ValA+Modulo)%Modulo;
    ValBM=(ValB+Modulo)%Modulo;
    ValCM=(ValC+Modulo)%Modulo;
    if (ValAM == 0) {  /* Linear equation: sol=-C/B */
      Sol1=Sol2=MultMod(Modulo-ValCM,ModInv(ValBM,Modulo),Modulo);
      }
    else {    /* Quadratic equation Ax^2+Bx+C=0 (mod F) */
      if (Modulo>2) {
        Sol1=MultMod(ValBM,ValBM,Modulo)-MultMod(4*ValAM,ValCM,Modulo);
        if (Sol1<0) {Sol1+=Modulo;}
               /* Find square root of Sol1 mod Modulo */
        if (Sol1 == 0) {                 /* if double root: sol = -b/2a */
          Sol1 = Sol2 = MultMod(ModInv((2*ValAM+Modulo)%Modulo,Modulo),((-ValBM)+Modulo)%Modulo,Modulo);
          }
        else {
          if (ModPow(Sol1,(Modulo-1)/2,Modulo) == 1) { /* if sols exist */
            if (Modulo%8 == 5) {
              s=ModPow(2*Sol1,(Modulo-5)/8,Modulo);
              Sol1=MultMod(MultMod(MultMod(MultMod(2*Sol1,s,Modulo),s,Modulo)-1,Sol1,Modulo),s,Modulo);
              }
            else {
              if (Modulo%8 != 1) {
                Sol1=ModPow(Sol1,(Modulo+1)/4,Modulo);
                }
              else {
                VarR=1;
                VarQ=Modulo-1;
                while (VarQ%2==0) {
                  VarQ/=2;
                  VarR*=2;
                  }
                VarX=2;
                while (true) {
                  VarY=ModPow(VarX,VarQ,Modulo);
                  if (ModPow(VarY,VarR/2,Modulo)!=1) {break;}
                  VarX++;
                  }
                VarX=ModPow(Sol1,(VarQ-1)/2,Modulo);
                VarV=MultMod(Sol1,VarX,Modulo);
                VarW=MultMod(VarV,VarX,Modulo);
                while (VarW!=1) {
                  VarK=1;VarD=VarW;
                  while (VarD!=1) {
                    VarD=MultMod(VarD,VarD,Modulo);
                    VarK*=2;
                    }
                  VarD=ModPow(VarY,VarR/VarK/2,Modulo);
                  VarY1=MultMod(VarD,VarD,Modulo);
                  VarR=VarK;
                  VarV=MultMod(VarV,VarD,Modulo);
                  VarW=MultMod(VarW,VarY1,Modulo);
                  VarY=VarY1;
                  }   /* end while */
                Sol1=VarV;
                }     /* end modulo 8 = 1 */
              }       
            s=ModInv((2*ValAM)%Modulo,Modulo);
            Sol2=MultMod((Modulo+Sol1-ValBM)%Modulo,s,Modulo);
            Sol1=MultMod((2*Modulo-Sol1-ValBM)%Modulo,s,Modulo);
            }
          else {   /* No solution exists */
            Sol1=Sol2=-1;
            }
          }
        }
      else {         /* Modulo <= 2 */
        if (Modulo==2) {
          switch ((int)ValBM * 2 + (int)ValCM) {
            case 0:           /* A = 1, B = 0, C = 0 */
              Sol1=Sol2=0;    /* Solution only for s=0 */
              break;
            case 1:           /* A = 1, B = 0, C = 1 */
              Sol1=Sol2=1;    /* Solution only for s=1 */
              break;
            case 2:           /* A = 1, B = 1, C = 0 */
              Sol1=0;         /* Solution for s=0 and s=1 */
              Sol2=1;
              break;
            default:          /* A = 1, B = 1, C = 1 */
              Sol1=Sol2=-1;   /* No solutions */
              break;
            }
          }                   /* End Modulo = 2 */
        else {                /* Modulo = 1 */
          Sol1=Sol2=0;
          }
        }
      }               /* End Quadratic Equation */
    ValAM=(ValA+ValF)%ValF;
    ValBM=(ValB+ValF)%ValF;
    ValCM=(ValC+ValF)%ValF;
    if (teach) {
      UL1 = UL;
      VL1 = VL+(MagnifY==1?"":"'");
      w("Let <B>"+UL1+" = s"+VL1+" - "+FP+"z</B>, so <B>[-(as"+sq+" + bs + c)/"+FP+"]"+VL1+sq+" + (2as + b)"+VL1+"z - a"+FP+"z"+sq+" = 1</B>.<P>So <B>");
      ShowEq(A,0,0,B,0,C,"s","");
      w("</B> should be multiple of <B>"+ValF+"</B>.<P>");
      }
    NbrEqs=EqNbr=0;
    byte sol=0;
    t=Sol1;
    for (cont=(Sol1<0?2:0); cont<2; cont++) {
      index=NbrFactors-1;v=mod[index];
      dif=0;
      q=(MultMod((MultMod(ValAM,t,ValF)+ValBM)%ValF,t,ValF)+ValCM)%ValF;  /* q%v = 0 */
      while (true) {
        if (q%v==0) {
          if (index==0) {          /* Solution found */
            NbrEqs++;
            if (teach) {
              s=t*mod[1];
              for (index2=1;index2<NbrFactors;index2++) {
                s+=pos[index2]*mod[index2+1];
                }
              s=s%ValF;
              if (sol==0) {
                sol=1;
                w("This holds for <B>s</B> = "+s);
                }
              else {
                w(", "+s);
                }
              }
            }
          else {
            pos[index]=t;
            t=0;
            for (index2=index;index2<NbrFactors;index2++) {
              t+=pos[index2]*mod[index2+1];
              }
            t=t%ValF;
            Dif[index]=dif;
            Q[index]=q;
            dif=MultMod((MultMod((2*t+v)%ValF,ValAM,ValF)+ValBM)%ValF,v,ValF);
            Pp=P[--index];
            t=0;v=mod[index];
            continue;
            }
          }
        if (index != NbrFactors-1 && ++t < factor[index]) {
          q = (q+dif)%ValF;
          dif = (dif+Pp)%ValF;
          continue;
          }
        else {
          while (++index < NbrFactors) {
            t=pos[index];v=mod[index];   /* Restore previous values */
            if (index < NbrFactors-1 && ++t < factor[index]) {
              Pp=P[index];
              dif=Dif[index];
              q=(Q[index]+dif)%ValF;
              dif = (dif+Pp)%ValF;
              break;
              }
            }
          if (index>=NbrFactors) {break;}
          }
        }
      if (Sol1==Sol2) {break;}  /* Do not process double root */
      t=Sol2;
      }
    if (teach) {
      if (sol==0) {
        w("No values of <B>s</B> makes the previous assertion true.<P>");
        SQD=Tmp1;
        System.arraycopy(Tmp2,0,NUM,0,6);
        System.arraycopy(Tmp3,0,DEN,0,6);
        DET=Tmp4;
        return;
        }
      w(".<P><UL>");
      }     /* end if (teach) */
    t=Sol1;
    for (cont=(Sol1<0?2:0); cont<2; cont++) {
      index=NbrFactors-1;v=mod[index];
      dif=0;
      q=(MultMod((MultMod(ValAM,t,ValF)+ValBM)%ValF,t,ValF)+ValCM)%ValF;  /* q%v = 0 */
      while (true) {
        if (q%v==0) {
          if (index==0) {          /* Solution found */
            EqNbr++;
            s=t*mod[1];
            for (index2=1;index2<NbrFactors;index2++) {
              s+=pos[index2]*mod[index2+1];
              }
            s=s%ValF;
                  //  Calculate biA = -((As+B)s+C)/F
                  //            biB = 2As+B
                  //            biC = -AF
                  //            biR = gcd(biA, biB, biC)
            LongToDoublePrecLong(A,biR);         // biR = A
            LongToDoublePrecLong(s,biS);         // biS = s
            LongToDoublePrecLong(B,biB);         // biB = B
            MultDoublePrecLong(biR,biS,biC);     // biC = As
            AddDoublePrecLong(biC,biB,biT);      // biT = As+B
            MultDoublePrecLong(biT,biS,biA);     // biA = (As+B)s
            LongToDoublePrecLong(C,biR);         // biR = C
            AddDoublePrecLong(biA,biR,biT);      // biT = (As+B)s+C
            LongToDoublePrecLong(-F,biR);        // biR = -F
            DivideDoublePrecLong(biT,biR,biA);   // biA = -((As+B)s+C)/F
            AddDoublePrecLong(biC,biC,biC);      // biC = 2As
            AddDoublePrecLong(biC,biB,biB);      // biB = 2As+B
            LongToDoublePrecLong(A,biT);         // biT = A
            MultDoublePrecLong(biR,biT,biC);     // biC = -AF
            GcdDoublePrecLong(biA,biB,biT);      // biT = gcd(biA, biB)
            GcdDoublePrecLong(biT,biC,biR);      // biR = gcd(biA, biB, biC)
            if (teach) {
              w("<LI>Let s = "+s+". Replacing in the above equation:<BR>");
              ShowBigEq(biA,biB,biC,VL1,"z");
              w(" = 1<P>");
              }
            if (biR[0]>1 || biR[1]>0 || biR[2]>0 || biR[3]>0 || biR[4]>0 ||
                biR[5]>0 && biR[5]<DosALa31) {  // if biR > 1 ...
              if (teach) {
                w("Since the <B>gcd</B> of the three coefficients is "+DoublePrecLongToString(biR)+" there are no integer solutions.<P>");
                }
              }
            else {
              LongArrayX = LongArrayY = null;
              GetRoot(biA,biB,biC);
              if (S.equals("")) {
                if (ContFrac(biA,(byte)3,(byte)1,s,T,MagnifY,A)) {
                  if (L2[0] < 0) {
                    ChangeSign(L1,(byte)1);
                    ChangeSign(L2,(byte)1);
                    }
                  LongArrayX = new ArrayLongs((int)(abs(L1[0])+1));
                  System.arraycopy(L1,
                                   0,
                                   LongArrayX.element,
                                   0,
                                   (int)(abs(L1[0])+1));
                  LongArrayY = new ArrayLongs((int)(abs(L2[0])+1));
                  System.arraycopy(L2,
                                   0,
                                   LongArrayY.element,
                                   0,
                                   (int)(abs(L2[0])+1));
                  }

                if (ContFrac(biA,(byte)3,(byte)(-1),s,T,MagnifY,A)) {
                  if (L2[0] < 0) {
                    ChangeSign(L1,(byte)1);
                    ChangeSign(L2,(byte)1);
                    }
                  if (teach) {
                    w("Choosing the solution with minimum <B>y</B>:<P><TABLE BORDER=1><TR><TH>");
                    NbrSols++;
                    }
                  if (LongArrayY != null && Compare(L2,LongArrayY.element) > 0) {
                    ShowLargeXY("X","Y",LongArrayX.element,LongArrayY.element,true,"","");
                    }
                  else {
                    ShowLargeXY("X","Y",L1,L2,true,"","");
                    }
                  w("</TABLE>");
                  }
                else {
                  if (LongArrayY != null) {
                    ShowLargeXY("X","Y",LongArrayX.element,LongArrayY.element,true,"","");
                    }
                  }
                }
              else {
                ContFrac(biA,(byte)4,(byte)1,s,T,MagnifY,A);
                ContFrac(biA,(byte)4,(byte)(-1),s,T,MagnifY,A);
                }
              }
            } 
          else {
            pos[index]=t;
            t=0;
            for (index2=index;index2<NbrFactors;index2++) {
              t+=pos[index2]*mod[index2+1];
              }
            t=t%ValF;
            Dif[index]=dif;
            Q[index]=q;
            dif=MultMod((MultMod((2*t+v)%ValF,ValAM,ValF)+ValBM)%ValF,v,ValF);
            Pp=P[--index];
            t=0;v=mod[index];
            continue;
            }
          }
        if (index < NbrFactors-1 && ++t < factor[index]) {
          q = (q+dif)%ValF;
          dif = (dif+Pp)%ValF;
          continue;
          }
        else {
          while (++index < NbrFactors) {
            t=pos[index];v=mod[index];   /* Restore previous values */
            if (index < NbrFactors-1 && ++t < factor[index]) {
              Pp=P[index];
              dif=Dif[index];
              q=(Q[index]+dif)%ValF;
              dif = (dif+Pp)%ValF;
              break;
              }
            }
          if (index>=NbrFactors) {break;}
          }
        }
      if (Sol1==Sol2) {break;}  /* Do not process double root */
      t=Sol2;
      }
    if (teach) {
      w("</UL>");
      }
    SQD=Tmp1;
    System.arraycopy(Tmp2,0,NUM,0,6);
    System.arraycopy(Tmp3,0,DEN,0,6);
    DET=Tmp4;
    } while (MagnifY*MagnifY > gcdAF);
   }
  }

private byte ShowSols(byte type) {
  if (type == 1) {
    if (CalculateKandL() != 0) {     /* if K or L not integers */
      return 1;                      /* bye */
      }
    } 
  if (teach) {
    w("m = ");
    ShowLargeNumber(H2);
    w("<BR>n = ");
    ShowLargeNumber(K2);
    w("<P>Using the formulas:<BR><B>P =");
    w(" m<BR>Q = -Cn<BR>");
    if (type == 1) {
      w("<TABLE><TR><TH>K =<TH>CD(P+S-2) + E(B-Bm-2ACn)<HR>4AC-B<SUP>2</SUP></TABLE><BR>");
      }
    w("R = An<BR>S = m + Bn<BR>");
    if (type == 1) {
      w("<TABLE><TR><TH>L =<TH>D(B-Bm-2ACn) + AE(P+S-2)<HR>4AC - B<SUP>2</SUP><TH> + Dn</TABLE><BR>");
      }
    w("</B>we obtain:<P><TABLE BORDER=1><TR><TH>");
    }
  w("P = ");
  ShowLargeNumber(H2);
  w("<BR>Q = ");
  MultLargeNumber(-C,K2,H1);
  ShowLargeNumber(H1);
  if (type==1) {
    w("<BR>K = ");
    ShowLargeNumber(L1);
    }
  w("<BR>R = ");
  MultLargeNumber(A,K2,H1);
  ShowLargeNumber(H1);
  w("<BR>S = ");
  MultAddLargeNumbers(1,H2,B,K2,H1);
  ShowLargeNumber(H1);
  if (type==1) {
    w("<BR>L = ");
    ShowLargeNumber(L2);
    }
  w(teach?"</TABLE><P>":"<P>");
  return 0;
  }

/* It returns zero if K and L are valid */
private byte CalculateKandL() {
  MultAddLargeNumbers(2,H2,B,K2,L1);  /* 2r + Bs */
  AddLargeLong(L1,-2,H1);  /* Kd = 2r + Bs - 2 */
  AddLargeLong(H2,-1,K1);  /* r - 1 */
  MultAddLargeNumbers(-B,K1,-2*A*C,K2,K1);  /* Ke */
  MultAddLargeNumbers(C*D,H1,E,K1,L1);  /* K(4AC - BB) */
  if (DivLargeNumber(L1,4*A*C-B*B,L1) != 0) {  /* K */
    return 1;               /* K not integer */
    }
  MultAddLargeNumbers(D,K1,A*E,H1,L2);
  if (DivLargeNumber(L2,4*A*C-B*B,L2) != 0) {
    return 1;               /* L not integer */
    }
  MultAddLargeNumbers(1,L2,D,K2,L2);    /* L */
  return (byte)0;
  }

private void ShowRecursionRoot(byte type) {
  byte t;
  t=ShowSols(type);
  if (type == 1) {
    ChangeSign(H2,(byte)1);
    ChangeSign(K2,(byte)1);
    if (t == 1 && ShowSols((byte)1) == 1) {
      if (teach) {
        w("m = ");
        ShowLargeNumber(H2);
        w("<BR>n = ");
        ShowLargeNumber(K2);
        w("<P>Using the formulas:<BR><B>P = m"+sq+" - ACn"+sq+"<BR>Q = -Cn(2m + Bn)<BR>K = -n(Em + CDn)<BR>R = An(2m + Bn)<BR>S = m"+sq+" + 2Bmn + (B"+sq+" - AC)n"+sq+"<BR>L = n(Dm + (BD-AE)n)<BR></B>we obtain:<TABLE BORDER=1><TR><TH>");
        }
      Mult2LargeNumbers(H2,H2,H1);   /* m^2 */
      Mult2LargeNumbers(H2,K2,K1);   /* mn */
      Mult2LargeNumbers(K2,K2,L1);   /* n^2 */
      MultAddLargeNumbers(1,H1,-A*C,L1,L2);
      w("P = ");
      ShowLargeNumber(L2);
      MultAddLargeNumbers(-2*C,K1,-B*C,L1,L2);
      w("<BR>Q = ");
      ShowLargeNumber(L2);
      MultAddLargeNumbers(-E,K1,-C*D,L1,L2);
      w("<BR>K = ");
      ShowLargeNumber(L2);
      MultAddLargeNumbers(2*A,K1,A*B,L1,L2);
      w("<BR>R = ");
      ShowLargeNumber(L2);
      MultAddLargeNumbers(B*B-A*C,L1,2*B,K1,L2);
      AddLargeNumbers(L2,H1,L2);
      w("<BR>S = ");
      ShowLargeNumber(L2);
      MultAddLargeNumbers(D,K1,B*D-A*E,L1,L2);
      w("<BR>L = ");
      ShowLargeNumber(L2);
      w(teach?"</TABLE><P>":"<P>");
      }
    }
  }

private void ShowRecursion(byte type) {
  String t1;
  long biA[] = new long[6];
  long biB[] = new long[6];
  long biC[] = new long[6];
  w("X<SUB>n+1</SUB> = P X<SUB>n</SUB> + Q Y<SUB>n</SUB>"+(type==1? " + K": ""));
  w("<BR>Y<SUB>n+1</SUB> = R X<SUB>n</SUB> + S Y<SUB>n</SUB>"+(type==1? " + L": ""));
  t1=" integer solution of the equation <B>m"+sq+" + bmn + acn"+sq+" = ";
  if (teach) {
    w("</TABLE><P>In order to find the values of P, Q, R, S we have to find first an"+t1);
    ShowEq(1,B,A*C,0,0,0,"m","n");
    w(" = 1</B>.<P>");
    }
  else {
    w("<P>");
    }
  LongToDoublePrecLong(A, biA);
  LongToDoublePrecLong(C, biB);
  MultDoublePrecLong(biA, biB, biC);
  LongToDoublePrecLong(1, biA);
  LongToDoublePrecLong(B, biB);
  GetRoot(biA, biB, biC);
  ContFrac(biA,(byte)2,(byte)1,0,0,1,A);
  if (teach) {
    w("An"+t1);
    ShowEq(1,B,A*C,0,0,0,"m","n");
    w(" = 1</B> is:<P>");
    }
  ShowRecursionRoot(type);
  if (B!=0) {
    if (teach) {
      w("<P>Another"+t1);
      ShowEq(1,B,A*C,0,0,0,"m","n");
      w(" = 1</B> is:<P>");
      }
    else {
      w("<P>as well as<P>");
      }
    MultAddLargeNumbers(1,H2,B,K2,H2); /* r <- r + Bs */
    ChangeSign(K2,(byte)1);            /* s <- -s */
    ShowRecursionRoot(type);
    }
  }

private boolean CheckMod(long R,long S,long X2,long X1,long X0) {
  long Y2=gcd(R,S);
  if (teach) {
    if (Y2!=1) {
      w("We try to check the equation modulo the prime divisors of ");
      if (R!=0 && S!=0) {
        w("gcd("+abs(R)+","+abs(S)+") = "+Y2+".<BR>");
        }
      else {
        w(Y2+".<BR>");
        }
      }
    }
  long Y1=2;
  int indH=1;
  long D1=abs(Y2);
  while(D1>=Y1*Y1) {
    int T=0;
    while(D1%Y1==0) {
      T++;
      if (T==1) {
        H1[indH++]=Y1;
        }
      D1/=Y1;
      }
    Y1++;
    if (Y1>3) {
      Y1++;
      }
    }
  if (D1>1) {
    H1[indH++]=D1;
    }
  for (int T=1;T<indH;T++) {
    if (H1[T]>1) {
      long Z=((X1*X1-4*X0*X2)%H1[T])+H1[T];
      long N=(H1[T]-1)/2;
      long Y=1;
      while (N!=0) {
        if (N%2!=0) {
          Y=(Y*Z)%H1[T];
          }
        N/=2;
        Z=(Z*Z)%H1[T];
        }
      if (Y>1) {
        if (teach) {
          w("There are no solutions modulo "+H1[T]+",");
          NoSol();
          }
        return true;
        }
      }
    }
  return false;
  }
}

class ArrayLongs {
  public long element[];

  public ArrayLongs(int NbrElements) {
    this.element = new long[NbrElements];
    }
  }

// </XMP>