import java.io.*;
import java.math.*;
//==========================================================C
//  OS(Quadratic Sieve) Function Program                    C
//  2次ふるい法javaプログラム ふるい改良3(改良2+(4)追加)    C
//                         +  行列改良2(行列をbitで表示)    C
// 　(1) ふるいは倍精度浮動小数点を使用し、本来のふるい処理 C
//   (2) 素数基底に-1を追加                                 C
//   (3) 素数基底のh倍の素数も2個存在すればふるいで採用     C
//   (4) f(x)=x^2-Nの基底での分解方法の改良(ふるいを応用)   C
//----------------------------------------------------------C
//   1. M = (int)sqrt(N)                                    C
//        using Bsqrt                                       C
//   2. Seting the Factor Base                              C
//        Prime Table:  PTB[NP]                             C
//        Factor Base:  FTB[NB]                             C
//   3. Sieving by QS                                       C
//     (1) Setting sieving base value                       C
//        Sieving base value: SBV[NS][3]                    C
//     (2) Siving for over factor base                      C
//        Sparase sieving Matrix: SQS[NT][NG]               C
//     (3) Sieving                                          C
//        Sieved Sparse Matrix: MQS[ND][spM][2]             C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/11/23              C
//        ( Tokyo Kougei University )                       C
//        後 保範： 東京工芸大学                            C
//==========================================================C
public class QS5{
//==========================================================C
// Global data  (Constant and Work area)                    C
//==========================================================C
// テーブル(素数、基底、ふるい)
 public final static int bno=32, sbM=250;
 public final static int spM=41, szM=7900, szN=8300, szwk=1000010; 
 public final static int szNP=150010, szNB=75010, szdt=300010;
 public static int PTB[] = new int[szNP];  // 素数テーブル
 public static int FTB[] = new int[szNB];  // 分解基底テーブル(-1は含まず)
 public static int TTB[] = new int[szNB];  // 分解基底の利用数テーブル(同上)
 public static int UTB[] = new int[szM];   // 利用基底テーブル(-1を含む)
 public static int SBV[][] = new int[szNB][3];  // ふるいテーブル(P,Q,S)
// f(x)=x^2-N分解テーブル (Add 2008.11.23)
 public int SBF[][] = new int[szN][2];  // 基底分解テーブル(S1,S2),なし:-1
 public int FCT[] = new int[spM];       // 分解基底要因、先頭は要因数
 public int FCN[][] = new int[spM][2];  // 同上+ベキ数、先頭は要因数
// ふるい実行テーブル、STB[*][0]は|X*X-N|, STB[*][1]はふるい累積値
 public static double STB[][] = new double[2][szwk]; // 基本基底+利用基底で分解されたふるい行列
 public static int MQS[][][] = new int[szN][spM][2]; // *[0]: 基底番号, *[1]: 基底のベキ数
 public static int MNO[] = new int[szN];     // ふるい行列の各行の数
 public static BigInteger MX[] = new BigInteger[szN];  // 対応する値X=M+x
// 基本基底以外の分解基底で分解されたふるいデータ 
 public static int MSN[][] = new int[szdt][2];  // L:分解基底の番号(>=NBS) (Change 2008.11.23)
 public static BigInteger MSV[] = new BigInteger[szdt]; // 対応する値X=M+x
// 行列計算用
 public static int VC[] = new int[szN];         // 軸交換ベクトル
 public static int VP[] = new int[szM];         // 枢軸ベクトル
 public static int AM[][] = new int[szN][sbM];  // 0-1行列(ビット圧縮)
 public static int BAM[][] = new int[szN][bno]; // ビット展開ワーク
 public static int BWK[] = new int[szN];        // ビット圧縮ワーク
// 因数分解用
 public static int LN[] = new int[szM];         // 従属する行番号
 public static int PW[] = new int[szM];         // 利用基底のベキ数
 public static int PK[] = new int[szNB];        // 分解ワーク
 public static BigInteger FactV[][] = new BigInteger[700][4]; // X*X-Y*Y=0 mod N
                            // [0]:X, [1]:Y; [2]:P=GCD(X+Y,N), [3];N/P
 public int NOS = 0, NH = 0;   // ふるいの総数、分解基底ふるいで得たデータ数
 public int NSD, NBU;       // 利用基底で得たデータ数、利用基底の数
 public int NNZ=0;          // 分解基底ふるいで取得に失敗したデータ数
 public long tims=0;
//==========================================================C
public BigInteger QS(BigInteger N, int NP, int LP, int h, int NV[])
//==========================================================C
//  OS(Quadratic Sieve) Program                             C
//  2次ふるい法javaプログラム ふるい改良１　                C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    NP      int,  In,  Size of Prime Table                C
//    LP      int,  In,  Size of one sieve deta             C
//    h       int,  In,  NBS=NB/h                           C
//    NV[5]   int,  Out, NV[0]=NB, NV[1]=NDD, NV[2]=NS      C
//                       NV[3]=NBS, NV[4]=ND                C
//    return  Bint, Out, M=sqrt(N)                          C
//----------------------------------------------------------C
//    NP  : 素数基底の総数                                  C
//    NB  : 分解基底の総数(-1は含まず)                      C
//    NBS : ふるいに使用する基本基底の数(=NB/h)             C
//    NBU : 利用基底の数(利用した分解基底で-1を含む)        C
//    NS  : ふるいテーブルの数                              C
//    ND  : 基本基底だけで得られたデータの数                C
//    NDD : ふるいで得られたデータの総数                    C
//    NH  : 全分解基底のふるいで得られたデータの数          C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/15              C
//    Ver.1  2008.11.23  Change from Bfact to Divide        C
//        ( Tokyo Kougei University )                       C
//        後 保範： 東京工芸大学                            C
//==========================================================C
 { int IWK[] = new int[100000];
   int NIWK = 100000, BV;
   int NBS, NBU=0;
   long t1, t2; 
   int Debug = 0; 
// M = (int)sqrt(N)
   t1 = System.currentTimeMillis();
   BigInteger M = Bsqrt(N);
   if(Debug==1) { System.out.println("M="+M.toString()); }
// Seting the Factor Base
   TPRIM(NP);
   int NB = BASE(N, NP);
   NBS = NP/(2*h);  if(h == 1) NBS = NB;
   BV = 11*FTB[NBS-1] + 1 ;
   int NS = SBASE(N, M, NBS, BV);
   t2 = System.currentTimeMillis();
   double T1 = (t2 - t1)*1.0e-3;
//   System.out.println("基底計算時間="+T1);
   if(Debug == 1) {
     for (int k=0; k<NB; k++)
       { System.out.println("k="+k+" "+FTB[k]); }
     for (int k=0; k<NS; k++)
       { System.out.println("SBV="+SBV[k][0]+" "+SBV[k][1]
                           +" "+SBV[k][2]); }
     } 
// Sieving by QS
   int ND = Qsieve(N, M, NB, NBS, NS, LP, h);
//  Set MSN, MSV and UTB
   int NDD = SetMQX(N, ND, NB, NBS);
   NV[0] = NB;   NV[1] = NDD;
   NV[2] = NS;   NV[3] = NBS;
   NV[4] = ND;  
   return M;}
//==========================================================C
public int Qsieve(BigInteger N, BigInteger M, int NB, 
                  int NBS, int NS,int LP, int h) 
//==========================================================C
//  Sieving by QS                                           C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    M       Bint, In,  M=sqrt(N)                          C
//    NB      int,  In,  Size of all Factor Base (FTB[NB])  C
//    NBS     int,  In,  Size of basic Factor Base          C
//    NS      int,  In,  Size of sieve base value (SBV)     C
//    LP      int,  In,  Size of one sieve deta             C
//    h       int,  In,  NBS=NB/h                           C
//    return  int,  out, Number of sieved data( > NB )      C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/15              C
//    Ver.1  2008.11.23  Change from Bfact to Divide        C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int k, SV, ID, NTD, NTS, ND=0;
  int out=0;
  double DV[] = new double[3];
  double cpu, TL = 5.0;
  BigInteger MS, BLP, M2N;
  String s1;
  long t1, t2, tt1, tt2, tt3, Tim1, Tim2, Tim3;
  int Debug = 0; 
// Initial set
  s1 = Integer.toString(LP);
  BLP = new BigInteger(s1);
  s1 = N.toString();
  DV[0] = Double.parseDouble(s1);
  s1 = M.toString();
  DV[1] = Double.parseDouble(s1);
  M2N = M.multiply(M);  M2N = M2N.subtract(N);
  s1 = M2N.toString();
  DV[2] = Double.parseDouble(s1);
//  Clear TTB
  for (k=0; k<NB; k++)
   { TTB[k] = 0; }
// Main Loop
  SV = 0;  ID = 1;
  Tim1 = Tim2 = Tim3 = 0;
  t1 = System.currentTimeMillis();
  while (ND+NSD < NBS+NBU+h+1){
     tt1 = System.currentTimeMillis();
     if(Debug > 0) System.out.println(ND+" "+NSD+" "+NBS+" "+NBU);
//  Sieving by QS
     Dsieve(DV, NS, SV*ID, LP);
//  Select and set matrix
     s1 = Integer.toString(SV*ID);
     MS = new BigInteger(s1); 
     MS = M.add(MS);
     tt2 = System.currentTimeMillis();
     Tim1 += tt2 - tt1;
//  Seving by QS 
     ND = Ssieve(LP, ND, N, MS, NB, NBS, SV*ID); 
     tt3 = System.currentTimeMillis();
     Tim2 += tt3 - tt2;
     if(Debug > 0)
      { for (int i=0; i<NB; i++)
         { System.out.print(TTB[i]+" "); } 
        System.out.println(); }
//  Setting NSD,NBU
     NSD = NBU = 0;
     for (k=NBS; k<NB; k++)
      { if(TTB[k] >= 2) 
         { NBU ++;
           NSD += TTB[k]; }
      }
     if(ID > 0) SV += LP; 
     ID = -ID;
     NOS += LP;
     t2 = System.currentTimeMillis();
     cpu = (t2 - t1)*1.0e-3;
     if(cpu > TL) {
        NTD = ND + NSD;  NTS = NBS + NBU;
        System.out.println("途中経過： 取得="+ NTD+", 基底="+NTS); 
        t1 = t2;  }
     Tim3 += t2 - tt3;
   }
 if(out > 0) {
   float CPU1 = (float)(Tim1*1.0e-3);
   float CPU2 = (float)(Tim2*1.0e-3);
   float CPU3 = (float)(Tim3*1.0e-3);
   float CPU4 = (float)(tims*1.0e-3);
   System.out.println("Sive Time="+CPU1+" "+CPU2+" "+CPU3);
   System.out.println("Ssive Time="+CPU4); }   
 return ND;}
//==========================================================C
public int Ssieve(int LP, int NO, BigInteger N,
       BigInteger MS, int NB, int NBS, int IVI)
//==========================================================C
//  Sieving by QS                                           C
//   select STB and set MQS(sieved sparce matrix)           C
//    STB[0][*] = X*X - N                                   C
//    STB[1][*] = multiply for base factor                  C
//----------------------------------------------------------C
//    LP      int,  In,  Size of one sieve deta             C
//    NO      int,  In,  Number of old selected sieve data  C
//    N       Bint, In,  Number for Factorizing             C
//    MS      Bint, In,  Start value (M+SV)                 C
//    NB      int,  In,  Size of Factor Base (FTB[NB])      C
//    NBS     int,  In,  Size of basic Factor Base          C
//    IVI     int,  In,  Start value (M+SV)                 C
//    return  int,  out, Number of new selected sieve data  C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/21              C
//    Ver.1  2008.11.23  Change from Bfact to Divide        C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int k, j, NN, IP, IV, VB, NZ, IVP; 
  String s1;
  BigInteger Zero = new BigInteger("0");
  BigInteger One = new BigInteger("1");
  BigInteger X, X2, NF;
  long t1, t2;
// Initial set
  t1 = System.currentTimeMillis();
  X = MS;  NN = NO;
  VB = FTB[NB-1]; 
// Main Loop
  for (k=0; k<LP; k++)
//  Select
   { IVP = IVI + k;
     int NOM = 0;
//     int SV = (int)(STB[0][k]/STB[1][k] + 0.5);
//     if(SV <= VB)
     if(STB[0][k] <= STB[1][k]*VB)
//  Set MQS matrix and MX vector
      { X2 = X.multiply(X);
        X2 = X2.subtract(N);
        if(X2.compareTo(Zero) < 0) 
         { MQS[NN][NOM][0] = 0; MQS[NN][NOM][1] = 1;
           NOM ++;  }
        X2 = X2.abs();
//        NF = BFact(X2, NBS);
        NF = Divide(X2, NBS, IVP);
        if(NF.compareTo(One) == 0)
          { int NEN = FCN[0][0];
            for (j=0; j<NEN; j++)
              { MQS[NN][j+NOM][0] = FCN[j+1][0] + 1;
                MQS[NN][j+NOM][1] = FCN[j+1][1];  }
//  Change 2008.11.23
//          { for (j=0; j<NBS; j++)
//             { if(PK[j] > 0) {
//                  MQS[NN][NOM][0] = j + 1;
//                  MQS[NN][NOM][1] = PK[j];
//                  NOM++; }
//              } 
            MX[NN] = X; MNO[NN] = NEN + NOM; NN++; }
        else
//  Set MSN, MSV and TTB vector
          { s1 = NF.toString();
            IV = Integer.parseInt(s1);
//            double sd = Double.parseDouble(s1);
            if(IV <= VB) {
              NZ = Search(IV, NB, FTB);
              if(NZ != -1)
                { MSN[NH][0] = NZ;  MSN[NH][1] = IVP;
                  MSV[NH] = X;
                  TTB[NZ] += 1;  NH ++; }
              else { NNZ++; }
              }
            else { NNZ++; } 
           }
       }
     X = X.add(One);
   }
  t2 = System.currentTimeMillis();
  tims += t2 - t1;
  return NN; }
//==========================================================C
public void Dsieve(double DV[], int NS, int SV, int LP)
//==========================================================C
//  Sieving by QS                                           C
//   set STB(Sieving table) using                           C
//    STB[0][*] = X*X - N                                   C
//    STB[1][*] = multiply for base factor                  C
//----------------------------------------------------------C
//    DV[3]   dbl,  In,  DV[0]=N, DV[1]=M, DV[2]=M^2-N      C
//    NS      int,  In,  Size of SBV(Sieving base value)    C
//    SV      int,  In,  Start Value for sieving            C
//    LP      int,  In,  Size of one sieve deta             C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/15              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int k, j, Q, ST;
  double DN, DM, DM2N, X, P;
// Initial set
  DN = DV[0];   DM = DV[1];
  DM2N = DV[2]; X = SV;
  for (j=0; j<LP; j++)
   { STB[0][j] = Math.abs(DM2N + X*(2.0*DM + X));
     STB[1][j] = 1.0;  X += 1.0; }
// Main Loop
  for (k=0; k<NS; k++)
    { P = SBV[k][0];
      Q = SBV[k][1];
      ST = SBV[k][2];
      ST = (ST - SV)%Q;
      if(ST < 0) ST += Q;
      for (j=ST; j<LP; j=j+Q)
       { STB[1][j] *= P;  }
    }
 }
//==========================================================C
public int SetMQX(BigInteger N, int ND, int NB, int NBS)
//==========================================================C
//  Setting MQX, MX and UTB by sparse matrix                C
//    V = TBL[L],  L is a output                            C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    ND      int,  In,  Number of sieved data in NBS       C
//    NB      int,  In,  Size of Factor Base (FTB[NB])      C
//    NBS     int,  In,  Size of basic Factor Base          C
//    return  int, out, Number of all sieved data           C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/21              C
//    Ver.1  2008.11.23  Change from Bfact to Divide        C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int k, j, l, NDD, LV, IVP;
  BigInteger X, X2;
  BigInteger Zero = new BigInteger("0");
//  Initial Set
  NDD = ND;
//  Set UTB
    UTB[0] = -1;
    for (k=0; k<NBS; k++)
     { UTB[k+1] = FTB[k]; }
    NBU = NBS + 1;
    for (k=NBS; k<NB; k++)
     { if(TTB[k] >= 2)
        { UTB[NBU] = FTB[k];
          NBU++; }
     }
//  Set MQX, MX
//   Set
    for (k=0; k<NH; k++)
     { l = MSN[k][0];  IVP = MSN[k][1];
       LV = FTB[l];
       if(TTB[l] >= 2)
        { int NOM = 0;
          X = MSV[k];  X2 = X.multiply(X);
          X2 = X2.subtract(N);
          if(X2.compareTo(Zero) < 0)
           { MQS[NDD][NOM][0] = 0; MQS[NDD][NOM][1] = 1;
             NOM ++;  }
          X2 = X2.abs();
//          BigInteger NF = BFact(X2, NBS);
          BigInteger NF = Divide(X2, NBS, IVP);
          if(Bmodi(NF, LV) == 0)
           { int NZ = Search(LV, NBU, UTB);
             int NEN = FCN[0][0];
             for (j=0; j<NEN; j++)
               { MQS[NDD][j+NOM][0] = FCN[j+1][0] + 1;
                 MQS[NDD][j+NOM][1] = FCN[j+1][1];  }
//  2008.11.23
//             for (j=0; j<NBS; j++)
//              { if(PK[j] > 0) {
//                  MQS[NDD][NOM][0] = j + 1;
//                  MQS[NDD][NOM][1] = PK[j];
//                  NOM++; }  }
             NOM = NEN + NOM;
             MQS[NDD][NOM][0] = NZ; 
             MQS[NDD][NOM][1] = 1;  NOM++; 
             MX[NDD] = X;  MNO[NDD] = NOM; NDD++; }
        } 
     }
  return NDD; }
//==========================================================C
public int Search(int V, int N, int TBL[])
//==========================================================C
//  Table search                                            C
//    V = TBL[L],  L is a output                            C
//----------------------------------------------------------C
//    V       int, in,  Searching Value                     C
//    N       int, in,  Size of TBL[]                       C
//    TBL[N]  int, in,  Searching Table                     C
//    return  int, out, Searched Number                     C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/08/30              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int l, m, h;
  int L = -1;
//  Search
  l = 0;  h = N;  m = h/2;
  while (h - l > 0)
    { if(h-l == 1)
       { if(V == TBL[l]) L = l;
         if(V == TBL[h]) L = h;
         break; }  
      if(V == TBL[m] ) { L = m;  break; }
      if(V > TBL[m] ) { l = m;  m = (h+m)/2; }
         else         { h = m;  m = (l+m)/2; }
     }
  return L; }
//==========================================================C
public void TPRIM(int NP)
//==========================================================C
//  Setting the Table of Prime Numbers                      C
//   PTB(k) = k's Prime Number                              C
//     k=0,1,...,NP-1,  PTB(1)=2, PTB(2)=3, PTB(3)=5,...    C
//----------------------------------------------------------C
//    NP     int, In,  Size of Prime Table                  C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/01/12              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int i, j, j1,k;
  int IWK[] = new int[10000];
  int NIWK=10000;
  int NO, LPM, NE, NH, NS, NS1; 
  int NIT[] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43};
//  Initial Set
  NO = 14;  LPM = 43;
  for (i=0; i<Math.min(NO,NP); i++)
    { PTB[i] = NIT[i]; }
//  Main Loop of prime numbers
  for (k=1; k<=10000; k++)
    { if(NP <= NO) { break; }
      NE = NIWK;
      if(LPM <= 30000) NE = Math.min(LPM*LPM/2, NIWK);
      for (i=0; i<NE; i++)
        { IWK[i] = 0; }
//  Set odd composit number
      NH = LPM + 2*NE;
      for (j1=1; j1<NO; j1++)
        { j = PTB[j1];
          if(j*j > NH) { break; }
          NS1 = j - (LPM+2)%j;
          if(NS1 == j) { NS1 = 0; }
          NS  = (j*(NS1%2) + NS1)/2;
          for (i=NS; i<NE; i+=j)
            { IWK[i] = 1; }
        }
//  Search prime number
      for (i=0; i<NE; i++)
        { LPM = LPM + 2;
          if(IWK[i] == 0)
            { PTB[NO] = LPM;
              NO += 1; }
          if(NO >= NP) { break; }
        }
     }
  }
//==========================================================C
public BigInteger Divide(BigInteger X, int NB, int IV)
//==========================================================C
//  Factorization with Primes by SBF                        C
//    X=FTB[FCN[1][0]]^FCN[1][1]*...*FTB[FCN[N-1][0]]^      C
//          FCN[N-1][1]*Y                                   C
//      N=FCN[0][0]                                         C
//----------------------------------------------------------C
//   X         Big, In,  Target value                       C
//   NB        int, In,  Size of Factor Base                C
//   IV        int, In,  Relative value; x=M+IV             C
//   FCN[][]   int, Out, Power of Factor Base               C
//   return    Big, Out, No factorized value (Y)            C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/11/23              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ BigInteger Y, T, S;
  int k, IW, NO;  String s1;
  BigInteger Zero = new BigInteger("0");
  Y = X;
// set FCT[k]; X=x^2-Nが分解される基底番号、FCT[0]は個数
  NO = 0;
  for (k=0; k<NB; k++)
    { IW = IV%FTB[k];
      if(IW < 0) { IW += FTB[k]; }
      if(IW==SBF[k][0] || IW==SBF[k][1])
        { FCT[NO+1] = k;  NO++; }
     }
  FCT[0] = NO;
// set T; TはFTB[FCT[k]]の乗算
  T = new BigInteger("1");
  for (k=1; k<=NO; k++)
    { s1 = Integer.toString(FTB[FCT[k]]);
      S = new BigInteger(s1);
      T = T.multiply(S); }
// set FCN[k][]; FCN[k][0]=基底番号、FCN[k][1]=ベキ数、FCN[0][0]は個数
  FCN[0][0] = NO;
  Y = Y.divide(T);
  for (k=1; k<=NO; k++)
    { IW = FTB[FCT[k]];
      FCN[k][0] = FCT[k];  FCN[k][1] = 1;
      s1 = Integer.toString(IW);
      S = new BigInteger(s1);
      for ( int j=0; j<30; j++)
        { T = Y.mod(S);
          int IT = T.compareTo(Zero);
          if(IT == 0) 
            { FCN[k][1] += 1;
              Y = Y.divide(S); }
          else { break; } 
   }   }
   return Y; }
//==========================================================C
public BigInteger BFact(BigInteger X, int NB)
//==========================================================C
//  Factorization with Primes                               C
//    X=FTB[0]^PW[0]*...*FTB[N-1]^PW[N-1]*Y                 C
//----------------------------------------------------------C
//   X         Big, In,  Target value                       C
//   NB        int, In,  Size of Factor Base                C
//   PW[NP]    int, Out, Power of Factor Base               C
//   return    Big, Out, No factorized value (Y)            C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/12/31              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ BigInteger Y, PS, P0;
  BigInteger Zero = new BigInteger("0");
  Y = X;
// Main Loop
  for (int k=0; k<NB; k++)
    { int IP = FTB[k];  PK[k] = 0;
      String s1 = Integer.toString(IP);
      PS = new BigInteger(s1);
      for ( int j=0; j<30; j++)
        { P0 = Y.mod(PS);
          int IT = P0.compareTo(Zero);
          if(IT == 0) 
            { PK[k] += 1;
              Y = Y.divide(PS); }
          else { break; } 
    }   }
  return Y; }
//==========================================================C
public int SBASE(BigInteger N, BigInteger M, int NB, int MV) 
//==========================================================C
//  Setting SBV: Sieving base vale                          C
//    if(S^2 - N (mod P^l) == 0) --> set SBV                C
//      SBV[*][0] = P in FTB                                C
//      SBV[*][1] = Q = P^l (l=1,2,...)                     C
//      SBV[*][2] = M-S (mod Q) (0 =< S < Q)                C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    M       Bint,In,  M=sqrt(N)                           C
//    NB      int, In,  Size of Factor Base                 C
//    MV      int, In,  Maximum value for sieving value     C
//    return  int, Out, Size of SBV(Sieving base value)     C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/27              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int NS = 0;
  int k, j, S2, Nin;
  int P, Q, IM, IP;
  double D, D2;
//  main Loop
  for (k=0; k<NB; k++)
    { P = FTB[k];  
      Q = P;  IP = 0;
      SBF[k][0] = SBF[k][1] = -1;
      while ((Q <= MV) & (Q > 0))
       { Nin = Bmodi(N, Q);
         IM  = Bmodi(M, Q);
         for (j=0; j<Q; j++)
          { D  = j + IM;
            D2 = D*D;
            S2 = (int)(D2 - Nin);
//            S2 = (j+IM)*(j+IM) - Nin;
            S2 = S2%Q;
            if(S2 < 0) S2 += Q;
            if(S2 == 0) 
             { SBV[NS][0] = P;
               SBV[NS][1] = Q;
               SBV[NS][2] = j;
               NS++;
               if(P == Q) 
                { SBF[k][IP] = j; IP++; } 
             }
          }
         Q *= P;
       }
    }  
  return NS; }
//==========================================================C
public int BASE(BigInteger N, int NP) 
//==========================================================C
//  Setting FTB: Factor Base                                C
//    if(N^(P-1)/2) (mod P) == 1) --> set FTB               C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    NP      int, In,  Size of Prime Table                 C
//    return  int, Out, Size of Factor Base                 C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/12              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int NB = 0;
  int k, P, P2;
  String s1, s2;
  BigInteger BP, BP2, Q;
  BigInteger One = new BigInteger("1");
//  main Loop
  for (k=0; k<NP; k++)
    { P  = PTB[k];   P2 = (P - 1)/2;
      s1 = Integer.toString(P);
      s2 = Integer.toString(P2);
      BP = new BigInteger(s1);
      BP2 = new BigInteger(s2);
      Q = N.modPow(BP2, BP);
      if(Q.compareTo(One) == 0) 
        { FTB[NB] = P;  NB++; }
    }  
  return NB; }
//==========================================================C
public int BASE1(BigInteger N, int NP) 
//==========================================================C
//  Setting FTB: Factor Base                                C
//    if(S^2 - N (mod P) == 0) --> set FTB                  C
//----------------------------------------------------------C
//    N       Bint, In,  Input Value for Factorizing        C
//    NP      int, In,  Size of Prime Table                 C
//    return  int, Out, Size of Factor Base                 C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/01/12              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int NB = 0;
  int k, j, S2, Nin;
//  main Loop
  for (k=0; k<NP; k++)
    { int P = PTB[k];
      Nin = Bmodi(N, P);
      for (j=0; j<P; j++)
        { S2 = j*j - Nin;
          if(S2%P == 0) 
            { FTB[NB] = P;  NB++;
              break; }
        }
    }  
  return NB; }
//==========================================================C
public int Bmodi(BigInteger X, int Y) 
//==========================================================C
//  Z = X (mod Y)                                           C
//    X: BigInteger,  Y,Z: int                              C
//----------------------------------------------------------C
//    X       Bint, In,  First input value                  C
//    Y       int,  In,  Second input value                 C
//    return  int,  Out, Z = X (mod Y)                      C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/01/12              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
{ int Z;
  BigInteger BY, BZ;
  String s1;
  s1 = Integer.toString(Y);
  BY = new BigInteger(s1);
//  Z = X (mod Y)
  BZ = X.mod(BY);
  s1 = BZ.toString();
  Z  = Integer.parseInt(s1);
  if(Z < 0) Z += Y;
  return Z;  }
//==========================================================C
//  Bsqrt(x) : Sqrt for x of BigInteger                     C
//----------------------------------------------------------C
//   x        Bint, In,  Input integer value                C
//   return   Bint, Out, Integer sqrt(x)                    C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/07/27              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public BigInteger Bsqrt(BigInteger x)
  { BigInteger u, w, t, TEN, One, Zero;
    int  k, n, m;
    Zero = new BigInteger("0");
    One  = new BigInteger("1");
    TEN  = new BigInteger("10");
//  double xs=sqrt(x) 
    String x1 = x.toString();
    double xd = Double.parseDouble(x1);
    double xs = Math.sqrt(xd);
//  int 10**(9*n)*is=sqrt(x)
    n = 0;
    for (k=0; k<=150; k++) {
       if(xs <= 1.0e9) break;
       n = n + 1;
       xs = xs/10.0; }
    int ix = (int)xs;
    String ss = Integer.toString(ix);
//  BigInteger t = a start value of sqrt(x) 
    t = new BigInteger(ss);
    for (k=1; k<=n; k++)
      { t = t.multiply(TEN); }
//  Loop Number for Newton Method
    m = 4;
    if(n >=8) m = 6; 
    if(n >=30) m = 8;
//  BigInteger t = sqrt(x) by Newton Method
    if(t.compareTo(Zero) != 0) {
      for (k = 1; k <= m; k++)
        { w = t.multiply(t);
          u = w.add(x);
          w = t.add(t);
          t = u.divide(w); }  }
    t = t.add(One);
    return t; }
//==========================================================C
//  Gauss Elimination Program with 0-1                      C
//    Make matrix AM(0-1) from MQS                          C
//    Gauss Elimination A with bits                         C
//      A=AM[N][M]                                          C
//----------------------------------------------------------C
//   N         int,  in,  Number of rows for AM             C
//   M         int,  in,  Number of columns for AM          C
//   return int,  Out, ID >= 1; Number of dependency        C
//                         = 0; no-dependency               C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/30              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public int Gauss(int N, int M) throws Exception 
{ int ID, IP, IW;
  int k, kb, kl, j, i, ii, kk, kg;
  int IO=0;
// Setting Output File
/*  FileOutputStream fo = new FileOutputStream("Gauss.txt");
  PrintWriter out = new PrintWriter(fo,true); */
// Make matrix AM(0-1) with bits from MQS
  int MBT = SetAM(N, M);
//  Output Matrix
/*  if(IO >=1)
   { out.println("Input");
     for (i=0; i<N; i++)
      { out.print(i+"  ");
        for (j=0; j<MBT; j++)
         { out.print(AM[i][j]+" "); }
        out.println("");
      } 
   } */
// Gauss Elimination
//  Set VC,VP
  for (i=0; i<N; i++)
   { VC[i] = i; }
  for (j=0; j<M; j++)
   { VP[j] = 0; }
//  Main Loop
  k = 0; kk = 0;
  for (kb=0; kb<MBT; kb++){
//  対象区間のbitをintに展開
    AMtoB(N, kk, kb);  kg = kk;
/*    if(IO >=1){
      out.println("BAM Before");
      for (i=0; i<N; i++)
       { for (j=0; j<bno; j++)
          { out.print(BAM[i][j]+" "); }
         out.println("");
        }
     } */
//  対象区間の消去処理
    for (kl=0; kl<bno; kl++) {
//   bit内の処理
       if(k >= M) { break; }
       IP = 0;
       if(BAM[kk][kl] == 0) {
//   Search Pivot
         IP = -1;
         for (i=kk+1; i<N; i++)
          { if(BAM[i][kl] == 1) 
             { IP = i;  break; }
          }
//   Exchange row for kk <--> IP 
         if(IP != -1)
          { IW = VC[IP];  VC[IP] = VC[kk];
            VC[kk] = IW;
//   Exchange BAM
            for (j=0; j<bno; j++)
             { IW = BAM[IP][j];
               BAM[IP][j] = BAM[kk][j];
               BAM[kk][j] = IW; }
//   Exchange AM
            for (j=0; j<MBT; j++)
             { IW = AM[IP][j];
               AM[IP][j] = AM[kk][j];
               AM[kk][j] = IW; } 
          }
        }
//   if(k-column vector is all zero) --> skip 
       if(IP != -1) {
          VP[k] = VC[kk];
          for (i=kk+1; i<N; i++)
           { if (BAM[i][kl] == 1)
//   Elimination of BAM
              { for (j=0; j<kl; j++)
                 { BAM[i][j] = BAM[i][j]^BAM[kk][j]; } 
                for (j=kl+1; j<bno; j++)
                 { BAM[i][j] = BAM[i][j]^BAM[kk][j]; } 
//   Elimination of AM
                for (j=0; j<kb; j++)
                 { AM[i][j] = AM[i][j]^AM[kk][j]; }
                for (j=kb+1; j<MBT; j++)
                 { AM[i][j] = AM[i][j]^AM[kk][j]; }
              }
           }
//    Output
/*        if(IO >=1){
             out.println("AM After, i="+i);
             for (ii=0; ii<N; ii++)
              { for (j=0; j<bno; j++)
                 { out.print(BAM[ii][j]+" "); }
                for (j=0; j<MBT; j++)
                 { out.print(" "+AM[ii][j]); }
                out.println(""); 
              }
           } */
//  count up kk for N
          kk++;
        }
//  count up k for M
      k++;
     }
//  対象区間のintをbitにしてAMに戻す
    BtoAM(N, kg, kb);
/*    if(IO >=1)
     { out.println("kb="+kb);
       for (i=kg; i<N; i++)
        { out.print(i+"  ");
          for (j=0; j<MBT; j++)
           { out.print(AM[i][j]+" "); }
          out.println("");
        }
       out.println("VC");
       for (i=0; i<N; i++)
        { out.print(VC[i]+" "); }
       out.println(""); 
       out.println("VP");
       for (j=0; j<M; j++)
        { out.print(VP[j]+" "); }
       out.println("");
     }  */
  } 
// Number of Dependency
  ID = 0;
  if(N-kk > 0) { ID = N - kk; }
  return ID; }
//==========================================================C
//  Set 0-1 matrix AM with bits                             C
//    Make 0-1 matrix AM from MQS                           C
//----------------------------------------------------------C
//   N         int, in,  Number of rows for AM              C
//   M         int, in,  Number of columns for AM           C
//   return    int, out, Number of columns/bits             C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/30              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public int SetAM(int N, int M)
{ int i, j, k, jj, MBT, MBB, Bit;
  MBT = (M + bno - 1)/bno;
  MBB = MBT*bno;
// Main loop
  for (k=0; k<N; k++)
//  set k's row 
   { for (j=0; j<MBB; j++)
      { BWK[j] = 0; }
     int NOM = MNO[k];
     for (i=0; i<NOM; i++)
      { j = M - MQS[k][i][0] - 1;
        BWK[j] = MQS[k][i][1]%2; }
//  Set AM with bits
     jj = 0;
     for (j=0; j<MBT; j++)
      { Bit = 0;
        for (i=0; i<bno; i++)
         { Bit = Bit << 1;
           Bit = Bit | BWK[jj]; jj++; }
        AM[k][j] = Bit;
      }
   }
  return MBT; }
//==========================================================C
//  対象列区間のbitをintに展開                              C
//    Change from AM[kk:N][kb] to BAM[kk:N][*]              C
//----------------------------------------------------------C
//   N         int, in,  Number of rows for AM              C
//   kk        int, in,  current row for AM, BAM            C
//   kb        int, in,  current column for AM              C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/24              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public void AMtoB(int N, int kk, int kb)
{ int k, j, js, Bit, One=1;
  for (k=kk; k<N; k++){
     Bit = AM[k][kb];
     for (j=0; j<bno; j++){
        js = bno - j - 1;
        BAM[k][j] = (Bit >>> js) & One;
        }
     }  
  }
//==========================================================C
//  対象列区間のintをbitに戻す                              C
//    Change from BAM[kk:N][*] to AM[kk:N][kb]              C
//----------------------------------------------------------C
//   N         int, in,  Number of rows for AM              C
//   kk        int, in,  current row for AM, BAM            C
//   kb        int, in,  current column for AM              C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/24              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public void BtoAM(int N, int kk, int kb)
{ int k, j, js, Bit, One=1;
  for (k=kk; k<N; k++){
     Bit = 0;
     for (j=0; j<bno; j++){
        Bit = Bit << 1;
        Bit = Bit | BAM[k][j];
        }
     AM[k][kb] = Bit;
     }  
  }
//==========================================================C
//  対象行区間のbitをintに展開                              C
//    Change from AM[k][0:M/bno] to BWK[0:M]                C
//----------------------------------------------------------C
//   M         int, in,  Number of columns for AM           C
//   kt        int, in,  current row for AM, BAM            C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/26              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public void AMtoBL(int M, int kt)
{ int i, j, k, js, Bit, MBT, One=1;
  MBT = (M + bno - 1)/bno;
  k = 0;
  for (j=0; j<MBT; j++){
     Bit = AM[kt][j];
     for (i=0; i<bno; i++){
        if(k >= M) { break; }
        js = bno - i - 1;
        BWK[k] = (Bit >>> js) & One;
        k++;  }
     }  
  }
//==========================================================C
//  Factor : Factorizeing by the QS                         C
//    using:  MQS[N][spM], FTB[M], AM[N][M], VC[N], VP[M]   C
//    output: LN[NO][], PW[NO][M]                           C
//----------------------------------------------------------C
//   NV      Bint, In,  Input Value for Factorizing         C
//   MV      Bint, In,  MV=sqrt(NV), Number of QS starting  C
//   N       int,  In,  Number of rows for MQS, EM          C
//   M       int,  In,  Number of columns for EM            C
//   NO      int,  In,  Number of factorized value          C
//----------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2008/09/30              C
//        ( Tokyo Kougei University )                       C
//==========================================================C
 public void Factor(BigInteger NV, BigInteger MV,
             int N, int M, int NO)
 { int k, j, i, NN, NOP, CN;
   BigInteger X, Y, XV, YV, YI;
   BigInteger One = new BigInteger("1");
   String s1;
//  Set to dependet column
   for (k=0; k<Math.min(NO,500); k++)
    { NN = 1;  int RNO = N - NO + k;
      AMtoBL(N, RNO);
      LN[NN] = VC[RNO];
      for (j=0; j<M; j++)
       { if(BWK[j] == 1)
          { NN++;  LN[NN] = VP[j];  } }
      LN[0] = NN; 
//  Set to power of base-factorr
      NOP = LN[0];
//  Set to PW
      for (i=0; i<M; i++)
       { PW[i] = 0; }
      for (j=1; j<=NOP; j++)
       { CN = LN[j];
         int NOM = MNO[CN];
         for (i=0; i<NOM; i++)
          { int iNO = MQS[CN][i][0];
            PW[iNO] += MQS[CN][i][1]; }
       }
//  Compute X, Y with X^2=Y^2 (mod NV)
      X = One;  Y = One;
      for (j=1; j<=NOP; j++)
       { CN = LN[j];
//   Compute X
         XV = MX[CN];
         X  = X.multiply(XV);
         X  = X.mod(NV); }
//   Compute Y
      for (i=1; i<M; i++)
       { int IPW = PW[i]/2;
         if(IPW >= 1) 
          { s1 = Integer.toString(IPW);
            YI = new BigInteger(s1);
            s1 = Integer.toString(UTB[i]);
            YV = new BigInteger(s1);
            YV = YV.modPow(YI, NV);
            Y  = Y.multiply(YV);
            Y  = Y.mod(NV);  }
       }
      FactV[k][0] = X;
      FactV[k][1] = Y;
//   Comput GCD(X+Y,NV) and GCD(X-Y,NV)
      XV = X.add(Y);
      X  = XV.gcd(NV);
      Y  = NV.divide(X);
      FactV[k][2] = X;
      FactV[k][3] = Y;
    }
  }
 }