I been reading a lot about quadratic sieve:
- here in math.stackexchange
- at wiki
- A Tale of Two Sieves by Carl Pomerance, the invented of the algorithm
- Factoring of large primes by Stefan Buttcher, he also got Java implementation
- And many others...
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?
- 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