mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2013-01-18, 16:18   #1
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default Quick smoothness test

I've started coding the MPQS from scratch, trying to modify my original code just didn't work.

I was wondering if there was a quicker way of testing if a number is smooth than I'm currently doing it.

Here's my algorithm:

1. For each prime p in the factor base:
1.1 Calculate X mod p
1.2 While the result is 0:
1.2.1 X = X / p
1.2.2 Recalculate X mod p, update result
1.2.3 If X = 1, break out of loop

I've profiled my code, 91.7% of the CPU time is spent on line 1.1, the division on line 1.2.1 accounts for only 0.2% of the CPU time.

Are there any clever ways of doing this? How could I improve this algorithm?
Sam Kennedy is offline   Reply With Quote
Old 2013-01-18, 16:56   #2
bsquared
 
bsquared's Avatar
 
"Ben"
Feb 2007

D4A16 Posts
Default

What is X? The full precision (ax + b)^2 - N value?
bsquared is offline   Reply With Quote
Old 2013-01-18, 17:47   #3
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

Quote:
Originally Posted by bsquared View Post
What is X? The full precision (ax + b)^2 - N value?
Yes, sorry should have specified that.
Sam Kennedy is offline   Reply With Quote
Old 2013-01-18, 18:13   #4
bsquared
 
bsquared's Avatar
 
"Ben"
Feb 2007

2×35×7 Posts
Default

Quote:
Originally Posted by Sam Kennedy View Post
Yes, sorry should have specified that.
One approach to make this step a lot faster is the following:

When you are done sieving, carry the final sieve step of each prime into the trial division stage. I.e., if the blocksize is B and prime p exits the sieve interval at B+k, take the value B+k into the trial division routine.

Then, rather than checking if X % p is 0, you can check if (B+k - x) % p is 0, where x is the sieve index you are currently looking at (that is used to compute X = (ax + b)^2 - N). This changes the operation from a multiple-precision division involving X to a single precision division that will be a lot faster.
bsquared is offline   Reply With Quote
Old 2013-01-18, 19:30   #5
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

As a small example, let's say the block size is 10, with the elements running 0 - 9. Starting at index 2, each value is divisible by 3, so the prime p exits the array at element 11.

B+k = 11

Then say I'm checking to see if the 5th element is divisible by 3 (which it should be!), I calculate (11 - 5) % 3 = 6 % 3 = 0.

The 7th element shouldn't be, just to check: (11 - 7) % 3 = 4 % 3 = 1

Have I understood correctly, or just got lucky with the numbers I picked?
Sam Kennedy is offline   Reply With Quote
Old 2013-01-18, 19:54   #6
bsquared
 
bsquared's Avatar
 
"Ben"
Feb 2007

2·35·7 Posts
Default

You got it.

Also handy is the fact that k is the starting offset of the prime in the next sieve block. So once you are done sieving and trial dividing block B1, subtract B from all of the values in your (B+k) array and the results are where you start sieving block B2.

Last fiddled with by bsquared on 2013-01-18 at 19:57
bsquared is offline   Reply With Quote
Old 2013-01-18, 20:41   #7
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

How much speed is to be gained from splitting into blocks of 32kb? In my old implementation larger arrays tended to be faster than smaller arrays, but my old implementation was pretty slow anyway.
Sam Kennedy is offline   Reply With Quote
Old 2013-01-18, 21:22   #8
bsquared
 
bsquared's Avatar
 
"Ben"
Feb 2007

65128 Posts
Default

Quote:
Originally Posted by Sam Kennedy View Post
How much speed is to be gained from splitting into blocks of 32kb? In my old implementation larger arrays tended to be faster than smaller arrays, but my old implementation was pretty slow anyway.
Difficult to say - it will be very implementation dependent. However, if you are frequently reading/modifying/writing a small data of data (i.e., a sieve block), it generally pays to keep that information close to the processor. You can't get much closer than L1 cache, and that's why 32kB is a popular choice for the size of the sieve block.
bsquared is offline   Reply With Quote
Old 2013-01-18, 22:36   #9
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

2·41 Posts
Default

I was thinking, I have two arrays (soln1 and soln2), which hold the start index for the sieving, instead, could I just compute (x - soln1) % p, where x is the sieve array index, and if the result is 0, do the division?

Is creating the blocks as simple as calculating 32768/sizeof(size_t), and using that many per block?
Sam Kennedy is offline   Reply With Quote
Old 2013-01-18, 22:43   #10
chalsall
If I May
 
chalsall's Avatar
 
"Chris Halsall"
Sep 2002
Barbados

24·5·7·17 Posts
Default

Quote:
Originally Posted by Sam Kennedy View Post
Is creating the blocks as simple as calculating 32768/sizeof(size_t), and using that many per block?
Why don't you try it, and see?
chalsall is offline   Reply With Quote
Old 2013-01-18, 23:27   #11
Sam Kennedy
 
Sam Kennedy's Avatar
 
Oct 2012

10100102 Posts
Default

Quote:
Originally Posted by chalsall View Post
Why don't you try it, and see?
I've still got quite a bit to code before I start with that
Sam Kennedy is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Quick primality test? guptadeva Miscellaneous Math 12 2017-12-28 00:09
Inverse of Smoothness Probability paul0 Math 6 2017-07-25 16:41
A quick question Pegos Information & Answers 6 2016-08-11 14:39
Quick & Dirty storm5510 Programming 37 2009-09-08 06:52
Smoothness test Washuu Math 12 2005-06-27 12:19

All times are UTC. The time now is 05:04.

Sun Apr 11 05:04:29 UTC 2021 up 2 days, 23:45, 1 user, load averages: 1.82, 1.88, 1.75

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, 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.