mersenneforum.org MPQS Trial Division Optimisation
 User Name Remember Me? Password
 Register FAQ Search Today's Posts Mark Forums Read

 2020-06-08, 21:50 #1 Sam Kennedy     Oct 2012 1228 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?
2020-06-09, 01:24   #2
R.D. Silverman

Nov 2003

26·113 Posts

Quote:
 Originally Posted by Sam Kennedy 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,
Don't do this. Instead, reconstruct the actual factorizations by resieving.
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:
 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.
Then there is a bug somewhere. Either q(x) mod p = 0 was solved incorrectly,
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:
 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.
You shouldn't even attempt the LA until you know you have enough smooth relations.

All I can say is that there is a bug somewhere.... It is impossible to debug from the
info given.

 2020-06-09, 10:11 #3 Sam Kennedy     Oct 2012 5216 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 2020-06-09 at 10:40
2020-06-09, 11:40   #4
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

17×83 Posts

Quote:
 Originally Posted by Sam Kennedy however if I calculate -1064 mod 109 I get 79 which doesn't equal soln1 or soln2.
My Pari gives:
Code:
? (-1064)%109
%12 = 26
?
It is a quite antique installation, but I believe in its result.

2020-06-09, 11:56   #5
Sam Kennedy

Oct 2012

2×41 Posts

Quote:
 Originally Posted by R. Gerbicz My Pari gives: Code: ? (-1064)%109 %12 = 26 ? It is a quite antique installation, but I believe in its result.
This is weird, I could have sworn I got 79 a bunch of times using the calculator that comes with windows, but now I'm getting 26, maybe I imagined it, maybe I was temporarily teleported to an alternate universe. My code on the other hand *is* getting 79 for some bizarre reason, once I figure that out everything should hopefully come together.

Thank you!

 2020-06-09, 15:58 #6 chris2be8     Sep 2009 5×383 Posts And using perl: Code: \$ perl -e 'print -1064%109,"\n"' 26 Just belt and braces. Chris
 2020-06-09, 17:10 #7 henryzz Just call me Henry     "David" Sep 2007 Cambridge (GMT/BST) 573210 Posts
 2020-06-09, 18:43 #8 jasonp Tribal Bullet     Oct 2004 348910 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.
2020-06-09, 21:01   #9
Till

"Tilman Neumann"
Jan 2016
Germany

1101000012 Posts

Quote:
 Originally Posted by jasonp 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.

Hey ho, I hope that I'm not the only Java friend over here^^ (despite certain shortcomings of the platform in high-performance computations)

2020-06-09, 21:38   #10
Sam Kennedy

Oct 2012

2×41 Posts

Quote:
 Originally Posted by henryzz
That was part of it, I managed to fix that but it turns out the issue was a mix of my factor base indices being off-by-one and my partial-combining routine interacting incorrectly with my trial division routine. It's all fixed now

Quote:
 Originally Posted by jasonp;547551 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.
jasonp, is that how msieve works? I was absolutely blown away by how fast it is! I know you've put a lot of effort into optimising it and it's very impressive. My MPQS code can factor a 60 digit semiprime in a couple of minutes, but yours does it in a couple of seconds, so now my goal is to figure out how to optimise my implementation as much as possible. How much faster is SIQS than MPQS? (I know it will be implementation dependent, I read somewhere it's about a factor of 2, does that sound about right?)

I would be very interested in discussing optimisations, would it be better for me to start a new thread on that subject?

Quote:
 Originally Posted by Till Hey ho, I hope that I'm not the only Java friend over here^^ (despite certain shortcomings of the platform in high-performance computations)
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 2020-06-09 at 21:47

 2020-06-10, 04:30 #11 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 2×3×151 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/en-us/rese...ote-letter.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).

 Similar Threads Thread Thread Starter Forum Replies Last Post yih117 Math 5 2018-02-02 02:49 mathPuzzles Math 8 2017-04-21 07:21 Peter Hackman Factoring 7 2009-10-26 18:27 SPWorley Math 8 2009-08-24 23:26 ewmayer Factoring 7 2008-12-11 22:12

All times are UTC. The time now is 23:24.

Tue Oct 27 23:24:37 UTC 2020 up 47 days, 20:35, 2 users, load averages: 1.72, 1.61, 1.67