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

[Gerald H]

(02-16-2019 07:57 PM)Albert Chan Wrote:  

(02-16-2019 07:09 PM)Gerald H Wrote:  "Set at random

N := 803189 * 485909

the product of two primes, how many bases are actually witnesses? Try running PRIME?, I will be very surprised if the 35s is fast enough for you to repeat the test until a 1 appears."

Re-reading your quote, I was mistaken.
I read a bit too fast, and think you are saying until 1 [witness] appears.

Changing the word witness to non-witness made the statement un-ambiguous.
"Until 1 pseduoprime appears" also work.

(02-16-2019 06:18 PM)Gerald H Wrote:  Should you inspect the programme you'll find the largest small factor that is tested for is

999999

If the code checked this far, why the need for SPRP test ?
Any 12 digits integer, with small factor upto 999999 checked, already proved primality.

Although the programme may check for

999999

as a factor it does not check for all factors upto 999999.

[Gerald H]

(02-17-2019 04:15 AM)Thomas Klemm Wrote:  Meanwhile I extended the search a bit:

Code:

       r     n       k
       12.02 1563151 130048
       12.02 1844267 153456
       16.02 1911397 119284
       16.03 1600117 99844
       30.05 1968557 65520
       32.05 1663213 51892
       32.06 1389581 43348
       35.77 1922801 53748
       40.06 2078737 51892
       40.07 1333603 33280
       43.71 1777793 40676
       53.41 2102437 39364
       60.09 1999759 33280
       60.10 1589531 26448
       60.11 1275347 21216
       70.10 2121013 30256
       70.12 1546021 22048
       84.15 1455451 17296
       93.50 1459009 15604
       …

For \(1563151=1021×1531\) or \(1844267=1109×1663\) it's realistic to hit a strong liar after a few attempts.

Cheers
Thomas

Thank you for the very interesting statistics, Thomas, you have done more work on that than I have for such small numbers.

Indeed, as the number tested becomes smaller a single test becomes less reliable.

On the 35s my interest is for numbers of the form

(5 or 6 digit prime)*(5 or 6 digit prime)

& for such the programme as is has been satisfactory.

It would be very nice of you if you could produce statistics for numbers in that range.

The programme at

http://www.hpmuseum.org/forum/thread-4236.html

contains some useful calculations & sadly no users, if indeed there are any, have suggested improvements - perhaps you could assist?

[Gerald H]

"Although the programme may check for

999999

as a factor it does not check for all factors upto 999999."

Not correct, sorry - To some degree the programme does check for all factors up to 999999.

eg If you enter

55021677489

the programme very quickly finds its largest factor & returns

0

[Albert Chan]

Trivia: searching backwards, found a particular bad composite, with many SPRP non-witnesses.

N = 999999 512881 = 881917 * 1133893 = A*B

For PRP test:
Total PRP non-witnesses = gcd(A-1,N-1) gcd(B-1,N-1) = 125988 * 125988 = 15872 976144
Percent of hitting PRP non-witness = (125988² - 2) / (N-3) ≈ 1.59%

For SPRP test:
Tried first 10 million bases, SPRP non-witnesses = 59708
Percent of hitting SPRP non-witness ≈ 59708/1e7 ≈ 0.60%

Note: 0.60% assumed statstical trend continued all the way to N-2

[Gerald H]

(02-17-2019 10:58 PM)Albert Chan Wrote:  Trivia: searching backwards, found a particular bad composite, with many SPRP non-witnesses.

N = 999999 512881 = 881917 * 1133893 = A*B

For PRP test:
Total PRP non-witnesses = gcd(A-1,N-1) gcd(B-1,N-1) = 125988 * 125988 = 15872 976144
Percent of hitting PRP non-witness = (125988² - 2) / (N-3) ≈ 1.59%

For SPRP test:
Tried first 10 million bases, SPRP non-witnesses = 59708
Percent of hitting SPRP non-witness ≈ 59708/1e7 ≈ 0.60%

Note: 0.60% assumed statstical trend continued all the way to N-2

A very nice example of a difficult case - But the chance of testing this number, or any other constructed to be difficult cases, when feeding the programme a random number is very small.

Similarly the likelihood of finding a number as in posting #43 is very small.

[fred_76]

Hi !

After a long period without programming I was finally able to complete the Miller-Rabin primality test on the HP35s. It runs very fast, even if the code is not optimized, but from 9 digits numbers and above. For numbers less than 9 digits, the optimized brut force is faster. For numbers above 11 digits, I still have problems to calculate the exponentiation without overflow therefore I can't conclude.

Compared to the "brut force V3" :

Quote:Run time
8 digits 99 999 989 1m 28s => 2m 09s (1.5x slower)
9 digits 999 999 937 4m 36s => 2m 14s (2.1x faster)
10 digits 9 999 999 967 14m 37s => 2m 39s (5.5x faster)
11 digits 99 999 999 977 46m 26s => 2m 57s (15.8x faster)
11.7 digits 500 000 000 023 => 3m 02s
12 digits 999 999 999 989 2h 32m 18s => not working
(2^249999999997 mod 999999999989 does not work well because numbers are too big, I shall change the algorithm)

[ttw]

One can cast out small numbers (perhaps by "wheels" rather than trial division; wheels do several divisions at once) then switch to SQUFOF or Pollard's Rho or Pollard's P-1 method.

However, I digress as I only read part of an answer. The above is for factoring. Prime testing is also interesting. There is a polynomial-time algorithm but it runs in m^6 power or so where m is the length of the candidate.

There are several good tests. The Selfridge-Pomerance-Wagstaff test is rather nice (I'll try it later this year sometime). The Wiki has good article on prime testing.

[fred_76]

TTW, the "optimized brut force" methods is a wheel (2 3 5 7), optimized to be as fast as possible on the HP35s. It will factorise the number when you repeat the program.

The last version is a Miller-Rabin method I tested last week end. The problem is the calculation of the power A^B mod N which fails on the HP35s when the prime number to test is slightly bigger than 500 000 000 000. I'm using the right-to-left algorithm.

[fred_76]

Miller Rabin optimized on my HP35s gives the following results. In brackets the speed compared to the brut force v3).

Code:

6 digits (999 983) in 0' 37" slower than brut force
7 digits (9 999 989) in 1' 9" slower than brut force
8 digits (99 999 989) in 1' 22" (x1.1)
9 digits (999 999 937) in 1' 27" (x3.2)
10 digits (9 999 999 967) in 2' 13" (x6.6)
11 digits (99 999 999 977) in 2' 29" (x19)
12 digits (999 999 999 989) in 2' 37" (x58)

Note that the HP48GX does 999,999,999,989 in ~ 31 seconds... it is 5x faster than the HP35S.

I use the Miller Rabin algorithm already detailed in some page before. I test only a small set of bases as it was proven that such small set is enough :

* if n < 2,047, it is enough to test a = 2;
* if n < 9,080,191, it is enough to test a = 31 and 73;
* if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
* if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;

The MR test is however to be used for numbers greater than about 90 000 000 as the brut force is faster for smaller number.

Well, that's quite fast for such a calculator...

[ttw]

One of the nice things about SQUFOF is that the intermediate results are of the order of SQRT(N) rather than N. It helps; one only needs really long stuff N itself.

[Albert Chan]

(04-10-2019 03:56 PM)fred_76 Wrote:  I use the Miller Rabin algorithm already detailed in some page before. I test only a small set of bases as it was proven that such small set is enough :

* if n < 2,047, it is enough to test a = 2;
* if n < 9,080,191, it is enough to test a = 31 and 73;
* if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
* if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;

Or, you can pick the bases from the records: http://miller-rabin.appspot.com

Assuming 12 digits range:

if n < 132,239, test a = 814494960528
if n < 360,018,361, test a = 1143370 and 2350307676
if n < 154,639,673,381, test a = 15 and 176006322 and 4221622697
if n < 47,636,622,961,201, test a = 2 and 2570940 and 211991001 and 3749873356

Do a = a (mod n). If a = 0 or 1 or n-1, skip the test and return as probable prime.

Edit: with little effort, you can use the bigger base for bigger range.

Example if n < 341,531, test a = 9,345,883,071,009,581,737 → reduced a = ((9345883e12 % n) + 71009581737) % n

[ttw]

Amusing divisibility tricks. Nearly all divisibility tests are the same. I haven't checked to see if a similar test would work for divisibility by various composites (which would allow quick trial divisibility test as part of a primality test or a factoring program.)[/align]


https://arxiv.org/pdf/1903.04903.pdf

[fred_76]

Here is the code of the MR primality check.

I've been inspired by the post of Thomas Klemm. I haven't used the Gerald H arithmetic library (far too complex to me to understand with all the imbricated XEQs... sorry Gerald). That leaves place for further optimisation.

Modular addition

Here we reduce first A to A≡N and B to B≡N then calculate
A-N+B (which can't overflow) ≡ N
we could save a calculation of A-N if A+B does not overflow, but the test takes longer than the substraction A-N.

Code:

Line    Code    Comment
A001    LBL A    X=N, Y=B, Z=A returns (A+B)≡N
A002    RMDR    
A003    x<>y    
A004    Last X    
A005    RMDR    
A006    Last X    
A007    -    
A008    Last X    
A009    Rv    
A010    +    
A011    R^    
A012    RMDR    
A013    RTN

Exemple of use :
RCL A
RCL B
RCL N
XEQ A
returns (A+B)≡N in X and N in Last X (y z t contain rubbish)

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
A*A = x²k²+xyk+xyk+y² = k(x²k+xy+xy)+y² Note that you can't use directly 2xyk as it sometimes overflows.

The code is slightly more complex as it treats differently the case (AxB)≡N and (AxA)≡N for speed.

Code:

Line    Code Comment
B001    LBL B    X=N, Y=B, Z=A
B002    CF 0    
B003    STO N    N=module
B004    Rv    
B005    x=y?    
B006    SF 0    
B007    ENTER    
B008    ENTER    
B009    1000000    
B010    STO K    k=10…
B011    INT/    
B012    STO U    
B013    CLx    
B014    Last X    
B015    RMDR    
B016    STO V    
B017    FS? 0    
B018    GTO a=b    
B019    CLx    
B020    Last X    
B021    INT/    
B022    STO X    
B023    CLx    
B024    Last X    
B025    RMDR    
B026    STO Y    
B027    RCL X    A<>B
B028    RCL*U    
B029    RCL N    
B030    RMDR    
B031    RCL*K    xuk
B032    RCL X    
B033    RCL*V    xv
B034    RCL N    
B035    XEQ A    xv+xuk
B036    RCL Y    
B037    RCL*U    yu
B038    RCL N    
B039    XEQ A    yu+xv+xuk
B040    RCL*K    k(yu+xv+xuk)
B041    RCL Y    
B042    RCL*V    yv
B043    RCL N    
B044    XEQ A    k(yu+xv+xuk)+yv
B045    RTN    
B046    RCL U    A=B
B047    RCL*V    uv
B048    ENTER    
B049    ENTER    
B050    RCL N    
B051    XEQ A    2uv
B052    RCL U    
B053    x²        u²
B054    RCL N    
B055    RMDR    
B056    RCL*K    ku²
B057    RCL N    
B058    XEQ A    ku²+2uv
B059    RCL*K    k(ku²+2uv)
B060    RCL V    
B061    x²        v²
B062    RCL N    
B063    XEQ A    k(ku²+2uv)+v²
B064    CF 0    
B065    RTN

Exemple of use :
RCL A
RCL B
RCL N
XEQ B
returns (A*B)≡N in X and N in Last X (y z t contain rubbish)

Modular exponentiation

This one is the regular "right-to-left" modular exponentiation.

Code:

Line    Code    Comment
C001    LBL C    N in X, B in Y, A in Z
C002    STO M    
C003    Rv    
C004    STO B    
C005    Rv    
C006    STO A    
C007    1    
C008    STO R    
C009    RCL B    
C010    x<=0?    
C011    GTO C032    
C012    2    
C013    RMDR    =1 if odd, 0 if even
C014    X>0?    if odd
C015    XEQ C026    
C016    RCL A    
C017    RCL A    
C018    RCL M    
C019    XEQ B    call modular *
C020    STO A    A = (A*A) ≡ M
C021    RCL B    
C022    2    
C023    INT/    
C024    STO B    
C025    GTO C010    
C026    RCL R    
C027    RCL A    
C028    RCL M    
C029    XEQ B    call modular *
C030    STO R    R = (A*R) ≡ M
C031    RTN    
C032    RCL R    
C033    RTN

Exemple of use :
RCL A
RCL B
RCL N
XEQ C
returns (A^B)≡N in X (yzt and last x contain rubbish)

Check if composite

As explained here.

Code:

Line    Code    Comment
D001    LBL D    X=N, Y=D, Z=S, T=A
D002    STO P    stores N in P
D003    Rv    
D004    STO D    stores D in D
D005    Rv    
D006    STO S    stores S in S
D007    x<>y    
D008    RCL P    
D009    RMDR    A≡ N
D010    RCL D    
D011    Last X    
D012    XEQ C    (A^D) ≡ N
D013    STO Z    stores in Z
D014    1    
D015    x=y?    
D016    GTO D031    
D017    1    
D018    RCL-P    
D019    RCL+Z    
D020    x=0?    
D021    GTO D031    
D022    RCL Z    
D023    RCL Z    
D024    RCL P    
D025    XEQ B    X²=N
D026    STO Z    
D027    DSE S    
D028    GTO D017    
D029    1        if true
D030    RTN    
D031    0        if false
D032    RTN

Prime routine

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

Code:

Line    Code    Comment
E001    LBL E    
E002    STO L    
E003    1    
E004    -    
E005    STO E    
E006    0    
E007    STO T    
E008    RCL E    
E009    2    
E010    RMDR    
E011    x>0?    
E012    GTO E019    
E013    RCL E    
E014    2    
E015    INT/    
E016    STO E    
E017    ISG T    increment T
E018    GTO E008    
E019    RCL L    
E020    132239    trigger for 1 base
E021    x>y?    
E022    GTO E088    
E023    Rv    
E024    360018361    trigger for 2 bases
E025    x>y?    
E026    GTO E075    
E027    Rv    
E028    154639673381        trigger for 3 bases
E029    x>y?    
E030    GTO E056    
E031    SF 4        4 bases
E032    SF 3    
E033    SF 2    
E034    SF 1    
E035    2    
E036    XEQ E095    
E037    CF 4    
E038    x>0?    
E039    GTO E100    
E040    2570940    
E041    XEQ E095    
E042    CF 3    
E043    x>0?    
E044    GTO E100    
E045    211991001    
E046    XEQ E095    
E047    CF 2    
E048    x>0?    
E049    GTO E100    
E050    3749873356    
E051    XEQ E095    
E052    CF 1    
E053    x>0?    
E054    GTO E100    
E055    GTO E108    
E056    SF 3        3 bases
E057    SF 2    
E058    SF 1    
E059    15    
E060    XEQ E095    
E061    CF 3    
E062    x>0?    
E063    GTO E100    
E064    176006322    
E065    XEQ E095    
E066    CF 2    
E067    x>0?    
E068    GTO E100    
E069    4221622697    
E070    XEQ E095    
E071    CF 1    
E072    x>0?    
E073    GTO E100    
E074    GTO E108    
E075    SF 2        2 bases
E076    SF 1    
E077    1143370    
E078    XEQ E095    
E079    CF 2    
E080    x>0?    
E081    GTO E100    
E082    2350307676    
E083    XEQ E095    
E084    CF 1    
E085    x>0?    
E086    GTO E100    
E087    GTO E108    
E088    SF 1        1 base
E089    814494960528    
E090    XEQ E095    
E091    CF 1    
E092    x>0?    
E093    GTO E100    
E094    GTO E108    
E095    RCL T    
E096    RCL E    
E097    RCL L    
E098    XEQ D    
E099    RTN    
E100    RCL L    
E101    CF 1    
E102    CF 2    
E103    CF 3    
E104    CF 4    
E105    STO C    
E106    VIEW C    
E107    RTN
E108    RCL L    
E109    STO P    
E110    VIEW P    
E111    RTN

Exemple of use :
RCL N
XEQ E
shows P= N if N is Prime
shows C= N if N is Composite
the flags 1-2-3-4 show the status of the calculations (4 bases max) so you see where it is
the flag 0 is used by the modular multiplication, it blinks so that you know it did not freeze :-)

One could go slightly faster by storing the constants that are used several times in variables (like 1000000, 1 and 2).

Adding brut force to find factors of for "small numbers"

1) I have also merged the "brut force" with the "MR" so that numbers smaller than 20 359 000 will be analysed by brut force, and greater by MR, because MR is slower than brut force for numbers smaller than 20 359 000.

2) When a number is found to be composite, it goes to brut force to find factors. When a factor is found, the quotient is sent back to MR to check if it is composite or prime, and so on. Brut force is however very slow to find factors if they are big. For example it takes only 34s for MR to say 999 999 998 479 is composite, but 2h30 for Brut Force to find its factors...

Fred

[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

[Albert Chan]

(04-11-2019 03:18 PM)fred_76 Wrote:  Adding brut force to find factors of for "small numbers"

1) I have also merged the "brut force" with the "MR" so that numbers smaller than 20 359 000 will be analysed by brut force, and greater by MR, because MR is slower than brut force for numbers smaller than 20 359 000.

2) When a number is found to be composite, it goes to brut force to find factors. When a factor is found, the quotient is sent back to MR to check if it is composite or prime, and so on. Brut force is however very slow to find factors if they are big. For example it takes only 34s for MR to say 999 999 998 479 is composite, but 2h30 for Brut Force to find its factors...

If you have built "MR", there is no point doing "brute force".
For example, I just created factor.lua, using Pollard's rho Algorithm.

My code do 3 deltas for efficiency, but only 1 is really needed.
Example, this is from my 7 years old laptop.

lua> p = require 'factor'

lua> n = 999999998479
lua> n, '=', table.concat(p.factor(n), ' * ') -- Lua 5.4 time ~ 1.8 msec
999999998479 = 999961 * 1000039

lua> n = 2^53-1
lua> n, '=', table.concat(p.factor(n), ' * ') -- Lua 5.4 time ~ 0.6 msec
9007199254740991 = 6361 * 69431 * 20394401

Edit: updated factor.lua by replacing MOD n for the delta, and use ABS instead. Speed up ~ 5%

Note: Luajit does not work properly with this code. It uses the naive a % b = a - floor(b/a) * a
This is my Floor-Mod patch for Luajit 1.18: https://github.com/LuaJIT/LuaJIT/issues/493

[Gerald H]

My assertion in post #35 is false.

The largest small factor tested for is 199, plus all squares up to 999999^2, so the first input requiring a Rabin test is

211*223,

ie 47,053