Fastest way to Factor (Prime-1)/2 for 64-bit Prime?
<p>I'm trying to gather some statistics on prime numbers, among which is the distribution of factors for the number (prime-1)/2. I know there are general formulas for the size of factors of uniformly selected numbers, but I haven't seen anything about the distribution of factors of one less than a prime.</p>

<p>I've written a program to iterate through primes starting at the first prime after 2^63, and then factor the (prime - 1)/2 using trial division by all primes up to 2^32. However, this is extremely slow because that is a lot of primes (and a lot of memory) to iterate through. I store the primes as a single byte each (by storing the increment from one prime to the next). I also use a deterministic variant of the Miller-Rabin primality test for numbers up to 2^64, so I can easily detect when the remaining value (after a successful division) is prime.</p>

<p>I've experimented using a variant of pollard-rho and elliptic curve factorization, but it is hard to find the right balance of between trial division and switching to these more complicated methods. Also I'm not sure I've implemented them correctly, because sometimes they seem to take a very lone time to find a factor, and based on their asymptotic behavior, I'd expect them to be quite quick for such small numbers.</p>

<p>I have not found any information on factoring many numbers (vs just trying to factor one), but it seems like there should be some way to speed up the task by taking advantage of this.</p>

<p>Any suggestions, pointers to alternate approaches, or other guidance on this problem is greatly appreciated.</p>

<hr>

<p>Edit:
The way I store the primes is by storing an 8-bit offset to the next prime, with the implicit first prime being 3. Thus, in my algorithms, I have a separate check for division by 2, then I start a loop:</p>

<pre class="lang-py prettyprint-override"><code>factorCounts = collections.Counter()
while N % 2 == 0:
    factorCounts[2] += 1
    N //= 2
pp = 3
for gg in smallPrimeGaps:
    if pp*pp &gt; N:
        break
    if N % pp == 0:
        while N % pp == 0:
            factorCounts[pp] += 1
            N //= pp
    pp += gg
</code></pre>

<p>Also, I used a wheel sieve to calculate the primes for trial division, and I use an algorithm based on the remainder by several primes to get the next prime after the given starting point.</p>

<hr>

<p>I use the following for testing if a given number is prime (porting code to c++ now):</p>

<pre class="lang-cpp prettyprint-override"><code>bool IsPrime(uint64_t n)
{
    if(n &lt; 341531)
        return MillerRabinMulti(n, {9345883071009581737ull});
    else if(n &lt; 1050535501)
        return MillerRabinMulti(n, {336781006125ull, 9639812373923155ull});
    else if(n &lt; 350269456337)
        return MillerRabinMulti(n, {4230279247111683200ull, 14694767155120705706ull, 1664113952636775035ull});
    else if(n &lt; 55245642489451)
        return MillerRabinMulti(n, {2ull, 141889084524735ull, 1199124725622454117, 11096072698276303650});
    else if(n &lt; 7999252175582851)
        return MillerRabinMulti(n, {2ull, 4130806001517ull, 149795463772692060ull, 186635894390467037ull, 3967304179347715805ull});
    else if(n &lt; 585226005592931977)
        return MillerRabinMulti(n, {2ull, 123635709730000ull, 9233062284813009ull, 43835965440333360ull, 761179012939631437ull, 1263739024124850375ull});
    else
        return MillerRabinMulti(n, {2ull, 325ull, 9375ull, 28178ull, 450775ull, 9780504ull, 1795265022ull});
}
</code></pre>

<p>This is how I store primes for later:
(I'm going to assume you want the factors of the number, and not just a primality test).</p>

<p>Copied from my website <a href="http://chemicaldevelopment.us/programming/2016/10/03/PGS.html" rel="nofollow">http://chemicaldevelopment.us/programming/2016/10/03/PGS.html</a></p>

<p>Iâm going to assume you know the binary number system for this part. If not, just think of 1 as a âyesâ and 0 as a ânoâ.</p>

<p>So, there are plenty of algorithms to generate the first few primes. I use the Sieve of Eratosthenes to compute a list.</p>

<p>But, if we stored the primes as an array, like [2, 3, 5, 7] this would take up too much space. How much space exactly?</p>

<p>Well, 32 bit integers which can store up to 2^32 each take up 4 bytes because each byte is 8 bits, and 32 / 8 = 4</p>

<p>If we wanted to store each prime under 2,000,000,000, we would have to store over 98,000,000,000. This takes up more space, and is slower at runtime than a bitset, which is explained below.</p>

<p>This approach will take 98,000,000 integers of space (each is 32 bits, which is 4 bytes), and when we check at runtime, we will need to check every integer in the array until we find it, or we find a number that is greater than it.</p>

<p>For example, say I give you a small list of primes: [2, 3, 5, 7, 11, 13, 17, 19]. I ask you if 15 is prime. How do you tell me?</p>

<p>Well, you would go through the list and compare each to 15.</p>

<p>Is 2 = 15?</p>

<p>Is 3 = 15?</p>

<p>. . .</p>

<p>Is 17 = 15?</p>

<p>At this point, you can stop because you have passed where 15 would be, so you know it isnât prime.</p>

<p>Now then, letâs say we use a list of bits to tell you if the number is prime. The list above would look like:</p>

<p>001101010001010001010</p>

<p>This starts at 0, and goes to 19</p>

<p>The 1s mean that the index is prime, so count from the left: 0, 1, 2</p>

<p><strong>001</strong>101010001010001010</p>

<p>The last number in bold is 1, which indicates that 2 is prime.</p>

<p>In this case, if I asked you to check if 15 is prime, you donât need to go through all the numbers in the list; All you need to do is skip to 0 . . . 15, and check that single bit.</p>

<p>And for memory usage, the first approach uses 98000000 integers, whereas this one can store 32 numbers in a single integer (using the list of 1s and 0s), so we would need
2000000000/32=62500000 integers.</p>

<p>So it uses about 60% as much memory as the first approach, and is much faster to use.</p>

<p>We store the array of integers from the second approach in a file, then read it back when you run.</p>

<p>This uses 250MB of ram to store data on the first 2000000000 primes.</p>

<p>You can further reduce this with wheel sieving (like what you did storing (prime-1)/2)</p>

<p>I'll go a little bit more into wheel sieve.</p>

<p>You got it right by storing (prime - 1)/2, and 2 being a special case.</p>

<p>You can extend this to <strong>p#</strong> (the product of the first <strong>p</strong> primes)</p>

<p>For example, you use <strong>(1#)*k+1</strong> for numbers <strong>k</strong></p>

<p>You can also use the set of linear equations <strong>(n#)*k+L</strong>, where <strong>L</strong> is the set of primes less than <strong>n#</strong> and 1 excluding the first <strong>n</strong> primes.</p>

<p>So, you can also just store info for <strong>6*k+1</strong> and <strong>6*k+5</strong>, and even more than that, because <strong>L={1, 2, 3, 5}{2, 3}</strong></p>

<p>These methods should give you an understanding of some the methods behind it.</p>

<p>You will need someway to implement this bitset, such as a list of 32 bit integers, or a string.</p>

<p>Look at: <a href="https://pypi.python.org/pypi/bitarray" rel="nofollow">https://pypi.python.org/pypi/bitarray</a> for a possible abstraction</p>

python c++ prime-factoring factoring