quadraticsieve.java

factorization.zip

      }

      // Compute the leading coefficient in biS.

      BigNumOps.LongToBigNbr(afact[0], biS);
      for (index = 1; index < s; index++) {
        BigNumOps.MultBigNbrByLong(biS, afact[index], biS);
      }
      for (index = 0; index < s; index++) {
        D = 1;
        E = afact[index];
        for (index2 = 0; index2 < s; index2++) {
          if (index != index2) {
            D = D * afact[index2] % E;
          }
        }
        amodq[index] = (int)D;
        Q = modsqrt[aindex[index]] * BigNumOps.modInv(D, E) % E;
        if (Q > E/2) {
          Q = E - Q;
        }
        BigNumOps.DivBigNbrByLong(biS, E, biR);
        BigNumOps.MultBigNbrByLong(biR, Q, aiJS[index]);
      }
      for (index = 0; index < Factorizer.NumberLength; index++) {
        biN[index] = aiJS[0][index];
      }
      for (index2 = 1; index2 < s; index2++) {
        BigNumOps.AddBigNbr(biN, aiJS[index2], biN);
      }
      for (index = 1; index < NbrPrimes; index++) {
        D = 1;
        E = prime[index];
        for (index2 = 0; index2 < s; index2++) {
          D = D * afact[index2] % E;
        }
        ainv[index] = BigNumOps.modInv(D, E);
        for (index2 = 0; index2 < s; index2++) {
          Bainv2[index2][index] = (int)(BigNumOps.RemDivBigNbrByLong(aiJS[index2], E) * 2 *
              ainv[index] % E);
        }
        D = BigNumOps.RemDivBigNbrByLong(biN, E);
        soln1[index] = (int)((SieveLimit + ainv[index]*(E + modsqrt[index] - D)) % E);
        soln2[index] = (int)((SieveLimit + ainv[index]*(E + E - modsqrt[index] - D)) % E);
      }
      do {
//        if (onlyFactoring) {
          polynomialsSieved++;
//        }

  //
  // Sieve stage
  //
      D = PolynomialIndex;
      index2 = 0;
      while ((D & 1) == 0) {
        D /= 2;
        index2++;
        }
      if (polyadd = ((D & 2) != 0)) {
        BigNumOps.AddBigNbr(biN, aiJS[index2], biN);
        BigNumOps.AddBigNbr(biN, aiJS[index2], biN);
        }
      else {
        BigNumOps.SubtractBigNbr(biN, aiJS[index2], biN);
        BigNumOps.SubtractBigNbr(biN, aiJS[index2], biN);
        }
      int [] rowBainv2 = Bainv2[index2];
     // Compute solutions for divisors of the leading coefficients
      for (index2 = 0; index2 < s; index2++) {
        index = aindex[index2];
        E = prime[index];
        D = BigNumOps.RemDivBigNbrByLong(Factorizer.TestNbr, E*E);
        Q = BigNumOps.RemDivBigNbrByLong(biN, E*E);
        soln1[index] = (int)(((D - Q*Q)/E * BigNumOps.modInv(amodq[index2]*Q%E, E)%E+E+SieveLimit)%E);
        soln2[index] = 2*SieveLimit+1;
        }
      X1 = 2 * SieveLimit;
      F1 = polyadd? 2 - rowBainv2[1]: rowBainv2[1];
      if ((soln1[1] += F1)%2 == 0) {
        SieveArray[0] = (byte)(logar[1] - G);
        SieveArray[1] = (byte)(-G);
        }
      else {
        SieveArray[0] = (byte)(-G);
        SieveArray[1] = (byte)(logar[1] - G);
        }
      F2 = 2;
      index = 2;
      while (true) {
        F = (int)prime[index];
        F3 = F2 * F;
        if (X1+1 < F3) {
          F3 = X1+1;
          }
        F4 = F2;
        while (F4*2 <= F3) {
          System.arraycopy(SieveArray, 0, SieveArray, F4, F4);
          F4 *= 2;
          }
        System.arraycopy(SieveArray, 0, SieveArray, F4, F3-F4);
        if (F3 == X1+1) {
          break;
          }
        Q1 = logar[index];
        F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
        for (index2=(soln1[index] += F1) % F; index2 < F3; index2 += F) {
          SieveArray[index2] += Q1;
          }
        for (index2=(soln2[index] += F1) % F; index2 < F3; index2 += F) {
          SieveArray[index2] += Q1;
          }
        index++;
        F2 *= F;
        }

      for (; index<smallPrimeUpperLimit; index++) {
        F = (int)prime[index];
        F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
        soln1[index] += F1;
        soln2[index] += F1;
        }

      int primesmall = (int)prime[index-1];
      int soln1small = soln1[index-1];
      int soln2small = soln2[index-1];
      for (; index<firstLimit; index++) {
        F = (int)prime[index];
        F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
        F2 = F+F; F3=F2+F; F4=F3+F;
        S1 = (soln1[index] += F1) % F;
        S2 = (soln2[index] += F1) % F;
        if ((G0 = S2 - S1) < 0) {
          S1 = S2;
          G0 = -G0;
          }
        G1 = G0+F; G2 = G1+F; G3 = G2+F;
        Q1 = logar[index];
        index2=X1/F4*F4+S1;
        do {
          SieveArray[index2] += Q1;
          SieveArray[index2+F] += Q1;
          SieveArray[index2+F2] += Q1;
          SieveArray[index2+F3] += Q1;
          SieveArray[index2+G0] += Q1;
          SieveArray[index2+G1] += Q1;
          SieveArray[index2+G2] += Q1;
          SieveArray[index2+G3] += Q1;
          } while ((index2 -= F4) >= 0);
        }
      for (; index<secondLimit; index++) {
        F = (int)prime[index];
        F1 = (polyadd? F - rowBainv2[index]: rowBainv2[index]);
        X2 = X1 - 4*F;
        F2 = F+F; F3=F2+F; F4=F2+F2;
        Q1 = logar[index];
        for (index2=(soln1[index] += F1)%F; index2 <= X2; index2 += F4) {
          SieveArray[index2] += Q1;
          SieveArray[index2+F] += Q1;
          SieveArray[index2+F2] += Q1;
          SieveArray[index2+F3] += Q1;
          }
        for (; index2 <= X1; index2 += F) {
          SieveArray[index2] += Q1;
          }
        for (index2=(soln2[index] += F1)%F; index2 <= X2; index2 += F4) {
          SieveArray[index2] += Q1;
          SieveArray[index2+F] += Q1;
          SieveArray[index2+F2] += Q1;
          SieveArray[index2+F3] += Q1;
          }
        for (; index2 <= X1; index2 += F) {
          SieveArray[index2] += Q1;
          }
        }
      for (; index<thirdLimit; index++) {
        F = (int)prime[index];
        Q1 = logar[index];
        F1 = (polyadd? F - rowBainv2[index]: rowBainv2[index]);
        for (index2=(soln1[index] += F1)%F; index2 <= X1; index2 += F) {
          SieveArray[index2] += Q1;
          }
        for (index2=(soln2[index] += F1)%F; index2 <= X1; index2 += F) {
          SieveArray[index2] += Q1;
          }
        }
      if (polyadd) {
        for (; index<NbrPrimes2; index++) {
          Q1 = logar[index];
          F1 = (F = (int)prime[index]) - rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F1 = (F = (int)prime[++index]) - rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F1 = (F = (int)prime[++index]) - rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F1 = (F = (int)prime[++index]) - rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          }
        for (; index<NbrPrimes; index++) {
          Q1 = logar[index];
          F1 = (F = (int)prime[index]) - rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          }
        }
      else {
        for (; index<NbrPrimes2; index++) {
          Q1 = logar[index];
          F = (int)prime[index];
          F1 = rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F = (int)prime[++index];
          F1 = rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F = (int)prime[++index];
          F1 = rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          F = (int)prime[++index];
          F1 = rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          }
        for (; index<NbrPrimes; index++) {
          Q1 = logar[index];
          F = (int)prime[index];
          F1 = rowBainv2[index];
          if ((F2 = (soln1[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          if ((F2 = (soln2[index] += F1) % F) < X1) {
            SieveArray[F2] += Q1;
            }
          }
        }
  //
  // Trial division stage
  //
      index2 = 2*SieveLimit+1;
      do {
        if (SieveArray[--index2] >= 0 && SieveArray[index2]<64) {
          if (SieveArray[index2] < logarsmall) {
            if ((index2-soln1small)%primesmall != 0) {
              if ((index2-soln2small)%primesmall != 0) {continue;}
              }
            }
          trialDivisions++;
          BigNumOps.MultBigNbrByLong(biS, index2 - SieveLimit, biT);
          BigNumOps.AddBigNbr(biT, biN, biT);
          BigNumOps.MultBigNbr(biT, biT, biR);
          BigNumOps.SubtractBigNbr(biR, Factorizer.TestNbr, biR);  // Number to factor: (Ax+B)^2-N
          for (i=0; i<Factorizer.NumberLength; i++) {
            biT[i] = biR[i];
            }
          // factor biR

          if (kk>1) {
            while (BigNumOps.RemDivBigNbrByLong(biR, kk*kk)==0) {
              BigNumOps.DivBigNbrByLong(biR, kk*kk, biR);
              }
            }
          F = Factorizer.NumberLength;              // Back up Factorizer.NumberLength
          for (index=0; index<s; index++) {
            BigNumOps.DivBigNbrByLong(biR, afact[index], biR);
            if (biR[Factorizer.NumberLength-1] == 0 && biR[Factorizer.NumberLength-2] < Constants.DosALa31 ||
                biR[Factorizer.NumberLength-1] == Constants.MaxUInt && biR[Factorizer.NumberLength-2] >= Constants.DosALa31) {
                Factorizer.NumberLength--;
              }
            }
          switch (Factorizer.NumberLength) {
            case 6:
              biR5 = biR[5];
            case 5:
              biR4 = biR[4];
            case 4:
              biR3 = biR[3];
            case 3:
              biR2 = biR[2];
            case 1:
            case 2:
              biR1 = biR[1];
              biR0 = biR[0];
            }
          nbrColumns = 0;
          if (Factorizer.NumberLength <= 2) {
            Divid = (biR1 << 32) | biR0;
            for (index=1; index<NbrPrimes; index++) {
              Divisor = prime[index];
              while (Divid % Divisor == 0) {
                Divid /= Divisor;
                }
              }
            }