I been reading a lot about quadratic sieve:

In all those places they are talking about factoring 100 digits numbers, as the the upper bounds of this algorithm. And so I been trying to achieve for months.

I came up with my own implementation, based on the knowledge of the mentioned above. But I can't bit the 60 digits limit, which takes me several hours to solve. Below is my algorithm, the bottleneck is in the sieving process. Am I doing it wrong? What else can I improve?

The Algorithm

  • First I start by building the prime base. I am using a predefined B_BOUNDS value to determinate how many primes I want to include in it. I iterate over all the numbers starting from 2, and using the primality test to check if they are prime(BigInteger.nextProbablePrime). I am using Legendre symbol to include only those primes that are in quadratic residue mod N.
  • Next I am building wheels, each wheel holds one prime, it tells what is the next number that the prime will divide. I use them during the sieving loop. In order to find the first two values of each prime I use Tonelli–Shanks algorithm, he called it ressol. Also the wheel calculates the prime log. This is in order to avoid unnecessary division during the sieve loop.
  • Now I start the sieve loop, I start checking values from the ceil root of N(the values that I am trying to factor) and define sieveVectorBound as the size of the biggest prime in the factory base, this is how many numbers I will process over each loop. I choose this size so each wheel will be used attlist once.
    • First I calculate the base log, it's the size ceil root of N plus current position, position starts from 0 and increased by sieveVectorBound each loop.
    • I create 2 arrays of doubles, logs and trueLogs with the size of sieveVectorBound.
    • I create another array of VectorData with the size of sieveVectorBound. VectorData contains information about the position where it been found, and the b-smooth vector already transformed to boolean, I use BitSet for it. Position is used to calculate $x$ and $y$ such that $x^2 - N = y$, it's important for the last step of the algorithm.
    • Now I iterate over all the wheels and sum all their logs in the logs array. I do it with the help of two important methods that wheel got:
      • testMove(long limit) it checks if the wheel already reached a limit position
      • nextLog() it increase the wheel current position and return the log for that position(At this point it's the same log all the time, but in early implementations I been returning the log multiplied by the biggest power of that prime that could divide that position)
    • Next I return the wheels to the position they been at when this loop started and make another iteration over them all
    • For each wheel I check if it's corresponding log[index] is bigger than the log of biggest prime in the prime base power two. If it's bigger I skip to the next wheel but if it's not, this number is either a prime or divided by one of the prime base primes.
    • So if not, I find the actual log of that position by taking it from the trueLogs array or calculating and saving it there if it's the first time. I compare the actual log to the log I calculated if it's equal within a small error of $0.0000001$, I know that this position is fully factored by the prime base and I add its VectorData to bSmoothVectors VectorData array to be used latter.
    • if it's not equal then it's a prime, so I add its VectorData to bigPrimesList.
    • We are still in the if case of not ignoring the wheel, we added the VectorData to the right list and now I turn on the bit of VectorData b-smooth vector at the current index of the wheels array
  • When enough big primes and b-smooth vectors been found, I construct and solve the matrix using gaussian elimination, I know it can be improved by using algorithms to find the null space, but I am using a very efficient implementation with bit calculation, so it's not the bottleneck of the algorithm
  • I make small optimization in the bigPrimeList, I only take those big primes that been found more than once, and if the big prime been found just twice, I xor, those two vectors to one, as it's the only way it will be part of the solution.

After solving the matrix, I extract the solution and it's pretty much it


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.