quadraticsieve.java

factorization.zip

    int [] matrixXmY = new int[k];
    int [] matrixCalc3 = new int[k];     // Matrix that holds temporary data
    int [] matrixTemp;
    int [] matrixCalc1 = new int[32];    // Matrix that holds temporary data
    int [] matrixCalc2 = new int[32];    // Matrix that holds temporary data
    int [] matr;
    long seed;
    int Temp, Temp1;
    int dimension = 0;
    int stepNbr = 0;
    int currentOrder, currentMask;
    int dim, maxdim;
    int minind, min, minanz;
    int [] rowMatrixB;

    newDiagonalSSt = oldDiagonalSSt = -1;

    // Initialize matrix X-Y with random data
    seed = 123456789L;
    for (i=matrixXmY.length-1; i>=0; i--) {
      matrixXmY[i] = (int)seed;
      seed = (seed * 62089911L + 54325442L) % DosALa31_1;
      matrixXmY[i] += (int)(seed * 6543265L);
      seed = (seed * 62089911L + 54325442L) % DosALa31_1;
    }
    // Compute matrix V(0) = A(X-Y)
    MatrixOps.MultiplyAByMatrix(matrixB, matrixXmY, matrixCalc3, matrixV);
    // Compute matrix Vt(0) * V(0)
    MatrixOps.MatrTranspMult(matrixV, matrixV, matrixVtV0);
    while (true) {
      oldDiagonalSSt = newDiagonalSSt;
      stepNbr++;
      // Compute matrix A * V(i)
      MatrixOps.MultiplyAByMatrix(matrixB, matrixV, matrixCalc3, matrixAV);
      // Compute matrix Vt(i) * A * V(i)
      MatrixOps.MatrTranspMult(matrixV, matrixAV, matrixVtAV);

      // If Vt(i) * A * V(i) = 0, end of loop
      for (i=matrixVtAV.length-1; i>=0; i--) {
        if (matrixVtAV[i] != 0) {break;}
      }
      if (i<0) {break;}               // End X-Y calculation loop

      // Selection of S(i) and W(i)

      matrixTemp = matrixWinv2;
      matrixWinv2 = matrixWinv1;
      matrixWinv1 = matrixWinv;
      matrixWinv = matrixTemp;

      mask = 1;
      for (j=31; j>=0; j--) {
        matrixD[j] = matrixVtAV[j];        //  D = VtAV
        matrixWinv[j] = mask;              //  Winv = I
        mask *= 2;
      }

      index = 31;
      indexC = 31;
      for (mask = 1; mask != 0; mask *= 2) {
        if ((oldDiagonalSSt & mask) != 0) {
          matrixE[index] = indexC;
          matrixF[index] = mask;
          index--;
        }
        indexC--;
      }
      indexC = 31;
      for (mask = 1; mask != 0; mask *= 2) {
        if ((oldDiagonalSSt & mask) == 0) {
          matrixE[index] = indexC;
          matrixF[index] = mask;
          index--;
        }
        indexC--;
      }
      newDiagonalSSt = 0;
      for (j=0; j<32; j++) {
        currentOrder = matrixE[j];
        currentMask = matrixF[j];
        for (k=j; k<32; k++) {
          if ((matrixD[matrixE[k]] & currentMask) != 0) {
            break;
          }
        }
        if (k<32) {
          i = matrixE[k];
          Temp = matrixWinv[i];
          matrixWinv[i] = matrixWinv[currentOrder];
          matrixWinv[currentOrder] = Temp;
          Temp1 = matrixD[i];
          matrixD[i] = matrixD[currentOrder];
          matrixD[currentOrder] = Temp1;
          newDiagonalSSt |= currentMask;
          dimension++;
          for (k=31; k>=0; k--) {
            if (k != currentOrder) {
              if ((matrixD[k] & currentMask) != 0) {
                matrixWinv[k] ^= Temp;
                matrixD[k] ^= Temp1;
              }
            }
          }                      // end for k
        }
        else {
          for (k=j; k<32; k++) {
            if ((matrixWinv[matrixE[k]] & currentMask) != 0) {
              break;
            }
          }
          i = matrixE[k];
          Temp = matrixWinv[i];
          matrixWinv[i] = matrixWinv[currentOrder];
          matrixWinv[currentOrder] = Temp;
          Temp1 = matrixD[i];
          matrixD[i] = matrixD[currentOrder];
          matrixD[currentOrder] = Temp1;
          for (k=31; k>=0; k--) {
            if ((matrixWinv[k] & currentMask) != 0) {
              matrixWinv[k] ^= Temp;
              matrixD[k] ^= Temp1;
            }
          }                      // end for k
        }                        // end if
      }                          // end for j
      // Compute D(i), E(i) and F(i)
      if (oldDiagonalSSt != -1) {    // S=I => F=0
        // F = -Winv(i-2) * (I - Vt(i-1)*A*V(i-1)*Winv(i-1)) * ParenD * S*St
        MatrixOps.MatrixMultiplication(matrixVt1AV1, matrixWinv1, matrixCalc2);
        index = 31;           // Add identity matrix
        for (mask = 1; mask != 0; mask *= 2) {
          matrixCalc2[index] ^= mask;
          index--;
        }
        MatrixOps.MatrixMultiplication(matrixWinv2, matrixCalc2, matrixCalc1);
        MatrixOps.MatrixMultiplication(matrixCalc1, matrixCalcParenD, matrixF);
        MatrixOps.MatrMultBySSt(matrixF, newDiagonalSSt, matrixF);
      }
      // E = -Winv(i-1) * Vt(i)*A*V(i) * S*St
      MatrixOps.MatrixMultiplication(matrixWinv1, matrixVtAV, matrixE);
      MatrixOps.MatrMultBySSt(matrixE, newDiagonalSSt, matrixE);

      // ParenD = Vt(i)*A*A*V(i) * S*St + Vt(i)*A*V(i)
      // D = I - Winv(i) * ParenD
      MatrixOps.MatrTranspMult(matrixAV, matrixAV, matrixCalc1); // Vt(i)*A*A*V(i)
      MatrixOps.MatrMultBySSt(matrixCalc1, newDiagonalSSt, matrixCalc1);
      MatrixOps.MatrixAddition(matrixCalc1, matrixVtAV, matrixCalcParenD);
      MatrixOps.MatrixMultiplication(matrixWinv, matrixCalcParenD, matrixD);
      index = 31;           // Add identity matrix
      for (mask = 1; mask != 0; mask *= 2) {
        matrixD[index] ^= mask;
        index--;
      }

      // Update value of X - Y
      MatrixOps.MatrixMultiplication(matrixWinv, matrixVtV0, matrixCalc1);
      MatrixOps.MatrixMultAdd(matrixV, matrixCalc1, matrixXmY);

      // Compute value of new matrix V(i)
      // V(i+1) = A * V(i) * S * St + V(i) * D + V(i-1) * E + V(i-2) * F
      MatrixOps.MatrMultBySSt(matrixAV, newDiagonalSSt, matrixCalc3);
      MatrixOps.MatrixMultAdd(matrixV, matrixD, matrixCalc3);
      MatrixOps.MatrixMultAdd(matrixV1, matrixE, matrixCalc3);
      if (oldDiagonalSSt != -1) {         // F != 0
        MatrixOps.MatrixMultAdd(matrixV2, matrixF, matrixCalc3);
      }

      // Compute value of new matrix Vt(i)V0
      if (stepNbr > 3) {
        // Vt(i+1)V(0) = Dt * Vt(i)V(0) + Et * Vt(i-1)V(0) + Ft * Vt(i-2)V(0)
        MatrixOps.MatrTranspMult(matrixD, matrixVtV0, matrixCalc1);
        MatrixOps.MatrTranspMult(matrixE, matrixVt1V0, matrixCalc2);
        MatrixOps.MatrixAddition(matrixCalc1, matrixCalc2, matrixCalc2);
        if (oldDiagonalSSt != -1) {          // F != 0
          MatrixOps.MatrTranspMult(matrixF, matrixVt2V0, matrixCalc1);
          MatrixOps.MatrixAddition(matrixCalc1, matrixCalc2, matrixCalc2);
        }
      }
      else {
        if (stepNbr == 1) {
          MatrixOps.MatrTranspMult(matrixCalc3, matrixV, matrixCalc2);   // V(1)t * V(0)
        }
        else {    // if stepNbr == 2 ...
          MatrixOps.MatrTranspMult(matrixCalc3, matrixV1, matrixCalc2);  // V(2)t * V(0)
        }
      }
      matrixTemp = matrixV2;
      matrixV2 = matrixV1;
      matrixV1 = matrixV;
      matrixV = matrixCalc3;
      matrixCalc3 = matrixTemp;
      matrixTemp = matrixVt2V0;
      matrixVt2V0 = matrixVt1V0;
      matrixVt1V0 = matrixVtV0;
      matrixVtV0 = matrixCalc2;
      matrixCalc2 = matrixTemp;
      matrixTemp = matrixVt1AV1;
      matrixVt1AV1 = matrixVtAV;
      matrixVtAV = matrixTemp;
    }          // end while

    // Find matrix V1:V2 = B * (X-Y:V)
    {
      int row;
      for (row = matrixB.length-1; row >= 0; row--) {
        matrixV1[row] = matrixV2[row] = 0;
      }
      for (row = matrixB.length-1; row >= 0; row--) {
        rowMatrixB = matrixB[row];
        for (index = rowMatrixB.length-1; index>=0; index--) {
          matrixV1[rowMatrixB[index]] ^= matrixXmY[row];
          matrixV2[rowMatrixB[index]] ^= matrixV[row];
        }
      }
    }
    maxdim=64;
    dim=0;
    while (dim<maxdim) {
      for (i=dim; i<maxdim; i++) {
        matr = (i>=32?matrixV1: matrixV2);
        mask = 1 << (31-(i & 31));
        vectorIndex[i] = -1;
        for (j=0; j<matrixV1.length; j++) {
          if ((matr[j] & mask) != 0) {
            vectorIndex[i] = j;
            break;
          }
        }
      }
      for (i=dim; i<maxdim; i++) {
        if (vectorIndex[i] < 0) {
          MatrixOps.colexchange(matrixXmY, matrixV, matrixV1, matrixV2, dim, i);
          vectorIndex[i] = vectorIndex[dim];
          vectorIndex[dim]=-1;
          dim++;
        }
      }
      if (dim==maxdim) {break;}
      min=vectorIndex[dim];
      minind=dim;
      for (i=dim+1; i<maxdim; i++) {
        if (vectorIndex[i] < min) {
          min=vectorIndex[i];
          minind=i;
        }
      }
      minanz=0;
      for (i=dim; i<maxdim; i++) {
        if (vectorIndex[i] == min) {
          minanz++;
        }
      }
      if (minanz>1) {
        for (i=minind+1; i<maxdim; i++) {
          if (vectorIndex[i]==min) {
            MatrixOps.coladd(matrixXmY, matrixV, matrixV1, matrixV2, minind, i);
          }
        }
      }
      else {
        maxdim--;
        MatrixOps.colexchange(matrixXmY, matrixV, matrixV1, matrixV2, minind, maxdim);
      }
    }
    dim=0;          // find linear independent solutions
    while (dim<maxdim) {
      for (i=dim; i<maxdim; i++) {
        matr = (i>=32?matrixXmY: matrixV);
        mask = 1 << (31-(i & 31));
        vectorIndex[i] = -1;
        for (j=0; j<matrixV1.length; j++) {
          if ((matr[j] & mask) != 0) {
            vectorIndex[i] = j;
            break;
          }
        }
      }
      for (i=dim; i<maxdim; i++) {
        if (vectorIndex[i] < 0) {
          maxdim--;
          MatrixOps.colexchange(matrixXmY, matrixV, matrixV1, matrixV2, maxdim, i);
          vectorIndex[i] = vectorIndex[maxdim];
          vectorIndex[maxdim]=-1;
        }
      }
      if (dim==maxdim) {break;}
      min=vectorIndex[dim];
      minind=dim;
      for (i=dim+1; i<maxdim; i++) {
        if (vectorIndex[i] < min) {
          min=vectorIndex[i];
          minind=i;
        }
      }
      minanz=0;
      for (i=dim; i<maxdim; i++) {
        if (vectorIndex[i] == min) {
          minanz++;
        }
      }
      if (minanz>1) {
        for (i=minind+1; i<maxdim; i++) {
          if (vectorIndex[i] == min) {
            MatrixOps.coladd(matrixXmY, matrixV, matrixV1, matrixV2, minind, i);
          }
        }
      }
      else {
        MatrixOps.colexchange(matrixXmY, matrixV, matrixV1, matrixV2, minind, dim);
        dim++;
      }
    }
    return matrixV;
  }

  /* */
  static String getTotals() {
      String s = "";
      if (smoothsFound > 0) {
            s = "SIQS: " + QuadraticSieve.polynomialsSieved +
                  " polynomials sieved\n      " +
                  QuadraticSieve.trialDivisions +
                  " sets of trial divisions\n      " +
                  smoothsFound +
                  " smooth congruences found\n      " +
                  totalPartials +
                  " partial congruences found\n      " +
                  partialsFound + " useful partial congruences";
      }
      return s;
  }

}