quadraticsieve.java

factorization.zip

                          BigNumOps.AddBigNbrModN(biR, biU, biR);
                          }
                        }
                      }
                    if (expParity != 0) {
                      // Check if the index is already in the row.
                      for (j=0; j<nbrColumns; j++) {
                        if (index <= rowMatrixB[j]) {
                          break;
                          }
                        }
                      if (j<nbrColumns && index == rowMatrixB[j]) {
                                              // Index already in row.
                        D = prime[rowMatrixB[j]];
                        Factorizer.NumberLength--;
                        BigNumOps.MultBigNbrByLongModN(biR, D, biR);
                        Factorizer.NumberLength++;
                        if (BigNumOps.RemDivBigNbrByLong(Factorizer.TestNbr, D) == 0) {
                          BigNumOps.DivBigNbrByLong(Factorizer.TestNbr, D, biU);
                          BigNumOps.AddBigNbrModN(biR, biU, biR);
                          }
                        // Delete entry from row.
                        for (k=j+1; k<nbrColumns; k++) {
                          rowMatrixB[k-1] = rowMatrixB[k];
                          }
                        nbrColumns--;
                        }
                      else {               // Index not in row.
                        // Insert entry in row.
                        for (k=nbrColumns; k>j; k--) {
                          rowMatrixB[k] = rowMatrixB[k-1];
                          }
                        rowMatrixB[j] = index;
                        nbrColumns++;
                        }
                      }
                    }
                 if ((biV[Factorizer.NumberLength-1] & 0x80000000L) != 0) {
                   BigNumOps.AddBigNbr(biV, Factorizer.TestNbr, biV);
                   }
                 if ((biW[Factorizer.NumberLength-1] & 0x80000000L) != 0) {
                   BigNumOps.AddBigNbr(biW, Factorizer.TestNbr, biW);
                   }
                 Factorizer.NumberLength--;
                 BigNumOps.MultBigNbrModN(biV, biW, biU);
                 Factorizer.NumberLength++;
                 if (InsertNewRelation(biR, biT, biU, nbrColumns,
                          matrixB, rowMatrixB, pp, vectLeftHandSide)) {
                    partialsFound++;
                    pp++;
                    ShowSIQSStatus(pp, matrixB.length, startTime);
                    if (pp == matrixB.length) {
                      i = EraseSingletons(matrixB, vectLeftHandSide, vectExpParity);
                      if (i!=0) {
                        Util.verbose(i + " singletons discarded");
                        pp -= i;
                      }
                      else {
                        Util.verbose("Solving congruence matrix using Block Lanczos algorithm");
                        if (LinearAlgebraPhase(NbrPrimes, matrixB, prime,
                              biT, biR, biU, vectExpParity, vectLeftHandSide)) {
                          return BigNumOps.BigIntToBigNbr(biT);        // Factor found
                        }
                        else {
                          Util.verbose("Linear dependences were found. Discarding 50 congruences...");
                          pp -= 50;               // Factor not found
                          }
                        }
                      }
                    break;
                    }
                  i = -2;   // Do not execute next if.
                  break;
                  }
                prev = i;
                i = matrixPartial[i][1];    // Get next index for same hash.
                }       // end while
              if (i==-1 && nbrPartials<matrixPartial.length) { // No match.
                // Add partial to table of partials.
                if (prev >= 0) {
                  matrixPartial[prev][1] = nbrPartials;
                  }
                else {
                  matrixPartialHashIndex[(int)(Divid & 0x7FE)/2] = nbrPartials;
                  }
                matrixPartial[nbrPartials][0] = (int)Divid;
                matrixPartial[nbrPartials][1] = -1;  // Indicate last index
                                                     // with this hash.
                BigNumOps.MultBigNbrByLong(biS, index2 - SieveLimit, biT);
                BigNumOps.AddBigNbr(biT, biN, biT);    // biT = Ax+B
                for (index=0; index<Factorizer.NumberLength; index++) {
                  matrixPartial[nbrPartials][index+2] = (int)biT[index];
                  }
                nbrPartials++;
                }
              }
            }
          }
        } while (index2 > 0);

  // Next polynomial
      PolynomialIndex++;
      } while (PolynomialIndex < NbrPolynomials);
    PolynomialSetNbr++;
    } while (true);
  }

  static void ShowSIQSStatus(int pp, int matrixBLength, long startTime) {
    long New,u;

//    if (TerminateThread) {
//      throw new ArithmeticException();
//    }
//    calcThread.yield();
    New=System.currentTimeMillis();
    if (Factorizer.OldTimeElapsed >= 0 &&
        Factorizer.OldTimeElapsed/1000 != (Factorizer.OldTimeElapsed + New - Factorizer.Old)/1000) {
      Factorizer.OldTimeElapsed += New-Factorizer.Old;
      Factorizer.Old = New;
      long t = Factorizer.OldTimeElapsed/1000;
      if (New - startTime > 5 && pp>10) {
        u = (New - startTime) * (matrixBLength-pp) / pp / 1000;
        Util.verbose(t/86400+"d "+(t%86400)/3600+"h "+((t%3600)/60)+"m "+(t%60)+
            "s     Congruences found: "+pp+" ("+(int)((float)(pp*100)/(float)matrixBLength)+
            "%)   End sieve in "+u/86400+"d "+(u%86400)/3600+"h "+((u%3600)/60)+"m "+(u%60)+"s");
      }
      else {
        Util.verbose(t/86400+"d "+(t%86400)/3600+"h "+((t%3600)/60)+"m "+(t%60)+
            "s     Congruences found: "+pp+" ("+(int)((float)(pp*100)/(float)matrixBLength)+"%)");
      }
    }
  }

  static boolean InsertNewRelation(long [] biR, long [] biT, long [] biU,
      int nbrColumns, int [][] matrixB, int [] rowMatrixB, int pp,
      long [][] vectLeftHandSide) {
    /* Convert negative numbers to the range 0 <= n < Factorizer.TestNbr */

    int i,k,index;

    if ((Factorizer.TestNbr[0] & 1) == 0) {
      BigNumOps.DivBigNbrByLong(Factorizer.TestNbr,2,Factorizer.TestNbr);
      for (i=0; i<Factorizer.NumberLength; i++) {
        biT[i] = 0;
      }
      BigNumOps.AddBigNbrModN(biR, biT, biR);
      Factorizer.NumberLength--;
      GCD.ModInvBigNbr(biR, biT, Factorizer.TestNbr);
      Factorizer.NumberLength++;
      BigNumOps.MultBigNbrByLong(Factorizer.TestNbr,2,Factorizer.TestNbr);
    }
    else {
      GCD.ModInvBigNbr(biR, biT, Factorizer.TestNbr);
    }
    if ((biU[Factorizer.NumberLength-1] & 0x80000000L) != 0) {
      BigNumOps.AddBigNbr(biU, Factorizer.TestNbr, biU);
    }

    // Compute biU / biR  (mod Factorizer.TestNbr)

    Factorizer.NumberLength--;
    BigNumOps.MultBigNbrModN(biU, biT, biR);
    Factorizer.NumberLength++;

    // Insert it only if it is different from previous relations.

    if (pp == 0) {     // First relation. Insert it always.
      matrixB[pp] = new int [nbrColumns]; /* Add relation to matrix B */
      for (k=0; k<nbrColumns; k++) {
        matrixB[0][k] = rowMatrixB[k];
      }
      vectLeftHandSide[0] = new long[Factorizer.NumberLength];
      for (index=0; index<Factorizer.NumberLength; index++) {
        vectLeftHandSide[0][index] = biR[index];
      }
      return true;
    }
    for (i=0; i<pp; i++) {
      if (nbrColumns > matrixB[i].length) {
        continue;
      }
      if (nbrColumns == matrixB[i].length) {
        for (k=0; k<nbrColumns; k++) {
          if (rowMatrixB[k] != matrixB[i][k]) {
            break;
          }
        }
        if (k==nbrColumns) {
          return false;              // Do not insert same relation.
        }
        if (rowMatrixB[k] > matrixB[i][k]) {
          continue;
        }
      }
      for (k=pp-1; k>=i; k--) {
        matrixB[k+1] = matrixB[k];
        vectLeftHandSide[k+1] = vectLeftHandSide[k];
      }
      break;
    }
    matrixB[i] = new int [nbrColumns]; // Add relation to matrix B.
    for (k=0; k<nbrColumns; k++) {
      matrixB[i][k] = rowMatrixB[k];
    }
    vectLeftHandSide[i] = new long[Factorizer.NumberLength];
    for (index=0; index<Factorizer.NumberLength; index++) {
      vectLeftHandSide[i][index] = biR[index];
    }
    return true;
  }

  static int EraseSingletons(int [][] matrixB, long [][] vectLeftHandSide,
      int [] vectExpParity) {
    int i,j,delta;
    int [] rowMatrixB;
    // Find singletons in matrixB.
    for (i=matrixB.length-1; i>=0; i--) {
      vectExpParity[i] = 0;
    }
    for (i=matrixB.length-1; i>=0; i--) {
      rowMatrixB = matrixB[i];
      for (j=rowMatrixB.length-1; j>=0; j--) {
        vectExpParity[rowMatrixB[j]]++;
      }
    }
    delta = 0;
    // Erase singletons from matrixB.
    for (i=0; i<matrixB.length; i++) {
      rowMatrixB = matrixB[i];
      for (j=rowMatrixB.length-1; j>=0; j--) {
        if (vectExpParity[rowMatrixB[j]] == 1) {    // Singleton found
          delta++;
          break;
        }
      }
      if (j<0) {                                    // Singleton not found
        matrixB[i-delta] = matrixB[i];
        vectLeftHandSide[i-delta] = vectLeftHandSide[i];
      }
    }
    return delta;
  }

  /************************/
  /* Linear algebra phase */
  /************************/
 static  boolean LinearAlgebraPhase(int NbrPrimes, int [][] matrixB, long [] prime,
      long [] biT, long [] biR, long [] biU,
      int [] vectExpParity, long[][] vectLeftHandSide) {
    int mask, i, j, index;
    int [] rowMatrixB;

    int [] matrixV = BlockLanczos(matrixB);
    for (mask=1; mask != 0; mask*=2) {
      BigNumOps.LongToBigNbr(1, biT);
      BigNumOps.LongToBigNbr(1, biR);
      for (i=matrixB.length-1; i>=0; i--) {
        vectExpParity[i] = 0;
      }
      for (i=matrixB.length-1; i>=0; i--) {
        if ((matrixV[i] & mask) != 0) {
          Factorizer.NumberLength--;
          BigNumOps.MultBigNbrModN(vectLeftHandSide[i], biR, biU);
          Factorizer.NumberLength++;
          for (j=0; j<Factorizer.NumberLength; j++) {
            biR[j] = biU[j];
          }
          rowMatrixB = matrixB[i];
          Factorizer.NumberLength--;
          for (j=rowMatrixB.length-1; j>=0; j--) {
            vectExpParity[rowMatrixB[j]] ^= 1;
            if (vectExpParity[rowMatrixB[j]] == 0) {
              if (rowMatrixB[j] == 0) {
                BigNumOps.SubtractBigNbr(Factorizer.TestNbr, biT, biT);    // Multiply biT by -1.
              }
              else {
                BigNumOps.MultBigNbrByLongModN(biT, prime[rowMatrixB[j]], biT);
              }
            }
          }
          Factorizer.NumberLength++;
        }
      }
      for (i=NbrPrimes-1; i>=0; i--) {
        if (vectExpParity[i] != 0) {
          break;
        }
      }
      BigNumOps.SubtractBigNbrModN(biR, biT, biR);
      GCD.GcdBigNbr(biR, TestNbr2, biT);
      index = 0;
      if (biT[0] == 1) {
        for (index=1; index<Factorizer.NumberLength; index++) {
          if (biT[index] != 0) {
            break;
          }
        }
      }
      if (index<Factorizer.NumberLength) {     /* GCD not 1 */
        for (index=0; index<Factorizer.NumberLength; index++) {
          if (biT[index] != TestNbr2[index]) {
            break;
          }
        }
        if (index<Factorizer.NumberLength) {     /* GCD not 1 */
          return true;
        }
      }
    }
    return false;
  }

  static final long DosALa31_1 = Constants.DosALa31 - 1;

  static int [] BlockLanczos(int [][] matrixB) {
    int i, j, k=matrixB.length;
    int oldDiagonalSSt,newDiagonalSSt;
    int index, indexC, mask;
    int [] matrixD = new int[32];
    int [] matrixE = new int[32];
    int [] matrixF = new int[32];
    int [] matrixWinv = new int[32];
    int [] matrixWinv1 = new int[32];
    int [] matrixWinv2 = new int[32];
    int [] matrixVtV0 = new int[32];
    int [] matrixVt1V0 = new int[32];
    int [] matrixVt2V0 = new int[32];
    int [] matrixVtAV = new int[32];
    int [] matrixVt1AV1 = new int[32];
    int [] matrixAV = new int[k];
    int [] matrixCalcParenD = new int[32];
    int [] vectorIndex = new int[64];
    int [] matrixV = new int[k];
    int [] matrixV1 = new int[k];
    int [] matrixV2 = new int[k];