[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
02-17-2019, 05:22 AM
Post: #41
 Gerald H Senior Member Posts: 1,438 Joined: May 2014
RE: [HP35s] Program for prime number (brut force)
(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."

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.
02-17-2019, 05:46 AM
Post: #42
 Gerald H Senior Member Posts: 1,438 Joined: May 2014
RE: [HP35s] Program for prime number (brut force)
(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

contains some useful calculations & sadly no users, if indeed there are any, have suggested improvements - perhaps you could assist?
02-17-2019, 04:33 PM (This post was last modified: 02-17-2019 04:37 PM by Gerald H.)
Post: #43
 Gerald H Senior Member Posts: 1,438 Joined: May 2014
RE: [HP35s] Program for prime number (brut force)
"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
02-17-2019, 10:58 PM (This post was last modified: 02-17-2019 11:47 PM by Albert Chan.)
Post: #44
 Albert Chan Senior Member Posts: 904 Joined: Jul 2018
RE: [HP35s] Program for prime number (brut force)
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
02-18-2019, 06:09 AM
Post: #45
 Gerald H Senior Member Posts: 1,438 Joined: May 2014
RE: [HP35s] Program for prime number (brut force)
(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.
04-09-2019, 03:52 PM (This post was last modified: 04-10-2019 07:31 AM by fred_76.)
Post: #46
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
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)
04-09-2019, 05:22 PM
Post: #47
 ttw Member Posts: 192 Joined: Jun 2014
RE: [HP35s] Program for prime number (brut force)
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.
04-10-2019, 07:26 AM (This post was last modified: 04-10-2019 07:32 AM by fred_76.)
Post: #48
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
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.
04-10-2019, 03:56 PM (This post was last modified: 04-10-2019 06:07 PM by fred_76.)
Post: #49
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
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...
04-10-2019, 09:24 PM
Post: #50
 ttw Member Posts: 192 Joined: Jun 2014
RE: [HP35s] Program for prime number (brut force)
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.
04-11-2019, 02:52 AM (This post was last modified: 04-13-2019 01:54 AM by Albert Chan.)
Post: #51
 Albert Chan Senior Member Posts: 904 Joined: Jul 2018
RE: [HP35s] Program for prime number (brut force)
(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
04-11-2019, 07:42 AM
Post: #52
 ttw Member Posts: 192 Joined: Jun 2014
RE: [HP35s] Program for prime number (brut force)
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
04-11-2019, 03:18 PM (This post was last modified: 04-11-2019 05:25 PM by fred_76.)
Post: #53
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
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.

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
04-13-2019, 06:24 PM (This post was last modified: 10-09-2019 11:08 PM by Albert Chan.)
Post: #54
 Albert Chan Senior Member Posts: 904 Joined: Jul 2018
RE: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
(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 endend

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
04-15-2019, 03:55 PM (This post was last modified: 04-19-2019 05:46 PM by Albert Chan.)
Post: #55
 Albert Chan Senior Member Posts: 904 Joined: Jul 2018
RE: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 2 Guest(s)