[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)

[Albert Chan]

(04-11-2019 03:18 PM)fred_76 Wrote:  Modular multiplication

Based on the fact that if you write :
A = xk+y
B = uk+u
where k is a base number, I took k = (999 999 999 999)^0.5 = 1 000 000

Then

A*B = (xk+y)*(uk+v) = xuk²+xvk+yuk+yv or better = k(xuk+xv+yu)+yv

MOD function is computational expensive.
You can rearrange the formula to reduce MOD count per multiply from 4 to 3

A*B = (xk + y)*(uk − v) = (xk uk) − (y v) + (y uk − v xk)

Base number must be powers of 10, otherwise MOD won't work.
Similarly, IEEE double must pick 2^26 as base, sqrt(2^53-1) ~ 94906266 won't work.
Since (2^26)^2 < 2^53-1, some adjustment is needed for huge IEEE double.

This is my mulmod(a, b, n) Lua code:

PHP Code:

local function mulmod(a, b, n)          -- a*b (mod n)
    local m = a*b
    if m <= 2^53 then return m%n end    -- a*b is exact
    if m > 4.e31 then a=n-a; b=n-b end  -- assume n >= max(a,b)
    local a2 = 2^26 - a % 2^26
    local b2 = b % 2^26
    a = a + a2
    b = b - b2
    m = a*b2 - a2*b
    a = (a*b)%n - (a2*b2)%n         -- a = [-n+1,n-1]
    if m >= 0 then m=m%n else m = n - (-m)%n end
    b = a + m                       -- b = [-n+1,2n-1]
    if b >= n then return a-n+m end
    if b >= 0 then return b else return b+n end
end 


Quote:Prime routine

I used the bases shown here (thank you Albert Chan).

To cover *full* 53 bits IEEE double range, I use above link with 5 best bases.
However, that does not quite reach 53 bits. One simple check pushed it over 53 bits

Snippet from my isprime(n)

PHP Code:

-- Below good for full 53 bits (2^53-1 = 9007199254740991). 
-- Other "holes": 9049203443108251, 9325578861942661, 9718260168425581, 10064706355310041 ...

return issprp(2, n)
   and issprp(4130806001517, n)
   and issprp((149795463772692064 % n) - 4, n)
   and issprp((186635894390467040 % n) - 3, n)
   and issprp((3967304179347716096 % n) - 291, n)
   and n ~= 7999252175582851 


Note that the huge numbers are represented *exactly* in binary.
Also, adjustment (-4,-3,-291) must be negative, to avoid 53-bits overflow.

lua> p = require 'nextprime'
lua> x = 2^53-1
lua> for i=1,10 do x=p.prevprime(x); print(x) end
9007199254740881
9007199254740847
9007199254740761
9007199254740727
9007199254740677
9007199254740653
9007199254740649
9007199254740623
9007199254740613
9007199254740571

Edit: nextprime.lua in https://github.com/achan001/PrimePi