20200608, 21:50  #1 
Oct 2012
2·41 Posts 
MPQS Trial Division Optimisation
I've been adding debugging statements to my MPQS implementation and I can't figure out why this is happening:
I'm currently trial dividing by every prime in the factor base, and as far as I understand if you are at index x in the sieve array, then x % p should equal either soln1 or soln2. However I am finding numbers which are divisible by p at indices which do not yield soln1 or soln2 as calculated above. If I restrict the trial division to primes in the factor base which satisfy x % p == soln1 or soln2, I appear to find fewer smooths and the linear algebra stage always fails. I did find some weird behaviour with C++ and taking the modulus of negative numbers, I've fixed that but I still have the same problem. I'm just wondering if anyone else came across a similar problem when implementing MPQS or if anybody has any idea what could be going wrong here? 
20200609, 01:24  #2  
Nov 2003
2^{2}×5×373 Posts 
Quote:
Instead of accumulating log(p) for all the primes as you resieve, you should SAVE the primes that hit a location known to be smooth. Quote:
or you are trial dividing an incorrect value > which means you have not computed q(x) correctly. Perhaps there is a sign error somewhere. Quote:
All I can say is that there is a bug somewhere.... It is impossible to debug from the info given. 

20200609, 10:11  #3 
Oct 2012
2×41 Posts 
Thank you very much for your response, I'm looking at things a lot closer and I have a specific example to show:
N = 523022617466601111760007224100074291200000001 A = 539044587247440481 B = 347870396507058651 At index 1064 (indices in the range: [30000,30000]): (((A * 1064) + B)^2  N) / A = 969667587191076007674875352, which factors into primes < 6991 which is smooth with F=10000. (The full factorisation just for reference is 1 * 2^3 * 3 * 41^2 * 109 * 379 * 1667 * 3041 * 3833 * 4283 * 6991) One prime where the solution doesn't appear to be correct is 109, I've checked the calculated modular square root: 7^2 = N mod 109, which is correct. The inverse of A mod 109 is: 16 soln1 = 16 * (7  B) mod 109 = 32 soln2 = 16 * (7  B) mod 109 = 26 I've calculated soln1 and soln2 by hand and my program calculates the same solutions, however if I calculate 1064 mod 109 I get 79 which doesn't equal soln1 or soln2. Is my method for calculating the solutions correct? I've taken it from Scott Contini's thesis. Last fiddled with by Sam Kennedy on 20200609 at 10:40 
20200609, 11:40  #4 
"Robert Gerbicz"
Oct 2005
Hungary
2·3^{2}·83 Posts 

20200609, 11:56  #5  
Oct 2012
2·41 Posts 
Quote:
Thank you! 

20200609, 15:58  #6 
Sep 2009
882_{16} Posts 
And using perl:
Code:
$ perl e 'print 1064%109,"\n"' 26 Chris 
20200609, 17:10  #7 
Just call me Henry
"David"
Sep 2007
Cambridge (GMT/BST)
2^{5}×5×37 Posts 
Is https://stackoverflow.com/questions/...gativenumbers the issue?

20200609, 18:43  #8 
Tribal Bullet
Oct 2004
3·1,181 Posts 
You should be prepared to manually force a remainder whose sign you want when the numerator of a modulo operation is negative, because a case can be made for both a positive and negative remainder. It might even be CPU dependent in C and C++.
If the sieve interval is [M,M] then you can sidestep the issue by translating the sieve roots mod p backwards by (M%p) and then sieving from 0 to 2*M. Numbers are always positive in this formulation and you can use unsigned arithmetic throughout (unless you are Till, since Java doesn't have unsigned anything :) I also find it conceptually simpler to iterate through one sieve rather than two sieves where the second one uses negated roots. 
20200609, 21:01  #9  
"Tilman Neumann"
Jan 2016
Germany
11×43 Posts 
Quote:
Hey ho, I hope that I'm not the only Java friend over here^^ (despite certain shortcomings of the platform in highperformance computations) 

20200609, 21:38  #10  
Oct 2012
2·41 Posts 
Quote:
Quote:
I would be very interested in discussing optimisations, would it be better for me to start a new thread on that subject? I feel like coding in C++ is like giving Manuel from Fawlty Towers instructions: you can never be 100% certain what's going to happen next. Java definitely has its benefits! Last fiddled with by Sam Kennedy on 20200609 at 21:47 

20200610, 04:30  #11 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
3^{2}·101 Posts 
https://en.wikipedia.org/wiki/Modulo_operation
For even more interesting reading. It's implementation defined in C (this bit me a few times while doing printer firmware, because our HP/UX machines gave different results than x86 and MIPS). This paper is referenced and is a nice read about 5 different consistent methods, how they work, and how to convert results from one into the others.. https://www.microsoft.com/enus/rese...oteletter.pdf So that's 5, but then we get the inconsistent, implementation defined, or just weird ones. Pari's divrem uses the Euclidian version, which seems like a good choice for its purpose. It's also pretty darn consistent, as we would like. Perl is ... interesting. It even changes behavior with "use integer" to whatever C does (i.e. implementation defined). 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Sublinear complexity of trial division?  yih117  Math  5  20180202 02:49 
Mersenne trial division implementation  mathPuzzles  Math  8  20170421 07:21 
trial division over a factor base  Peter Hackman  Factoring  7  20091026 18:27 
P95 trial division strategy  SPWorley  Math  8  20090824 23:26 
Need GMP trialdivision timings  ewmayer  Factoring  7  20081211 22:12 