mersenneforum.org  

Go Back   mersenneforum.org > Factoring Projects > Factoring

Reply
 
Thread Tools
Old 2020-06-08, 21:50   #1
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default 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?
Sam Kennedy is offline   Reply With Quote
Old 2020-06-09, 01:24   #2
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

22·5·373 Posts
Default

Quote:
Originally Posted by Sam Kennedy View Post
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.
R.D. Silverman is offline   Reply With Quote
Old 2020-06-09, 10:11   #3
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

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
Sam Kennedy is offline   Reply With Quote
Old 2020-06-09, 11:40   #4
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

7·199 Posts
Default

Quote:
Originally Posted by Sam Kennedy View Post
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.
R. Gerbicz is offline   Reply With Quote
Old 2020-06-09, 11:56   #5
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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!

Sam Kennedy is offline   Reply With Quote
Old 2020-06-09, 15:58   #6
chris2be8
 
chris2be8's Avatar
 
Sep 2009

2×33×5×7 Posts
Default

And using perl:
Code:
$ perl -e 'print -1064%109,"\n"'
26
Just belt and braces.

Chris
chris2be8 is offline   Reply With Quote
Old 2020-06-09, 17:10   #7
henryzz
Just call me Henry
 
henryzz's Avatar
 
"David"
Sep 2007
Cambridge (GMT/BST)

165A16 Posts
Default

Is https://stackoverflow.com/questions/...gative-numbers the issue?
henryzz is offline   Reply With Quote
Old 2020-06-09, 18:43   #8
jasonp
Tribal Bullet
 
jasonp's Avatar
 
Oct 2004

3,529 Posts
Default

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.
jasonp is offline   Reply With Quote
Old 2020-06-09, 21:01   #9
Till
 
Till's Avatar
 
"Tilman Neumann"
Jan 2016
Germany

2·32·23 Posts
Cool

Quote:
Originally Posted by jasonp View Post
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)
Till is offline   Reply With Quote
Old 2020-06-09, 21:38   #10
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

10100102 Posts
Default

Quote:
Originally Posted by henryzz View Post
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 View Post
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
Sam Kennedy is offline   Reply With Quote
Old 2020-06-10, 04:30   #11
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

5×181 Posts
Default

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).
danaj is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Sublinear complexity of trial division? yih117 Math 5 2018-02-02 02:49
Mersenne trial division implementation mathPuzzles Math 8 2017-04-21 07:21
trial division over a factor base Peter Hackman Factoring 7 2009-10-26 18:27
P95 trial division strategy SPWorley Math 8 2009-08-24 23:26
Need GMP trial-division timings ewmayer Factoring 7 2008-12-11 22:12

All times are UTC. The time now is 13:07.

Sat Sep 19 13:07:09 UTC 2020 up 9 days, 10:18, 1 user, load averages: 1.98, 1.61, 1.50

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.