(11C) Random Prime Number Generator

[Gamo]

Program to Generate Random Prime Number.

Procedure: Fix 0

[E] display Prime Number

Example:
This program use Random Number Generator range approximately
from 1000 to 9999 and some time the result will get below 1000.


Below result is randomly difference on each run

[E] display 7247
[E] display 7013
[E] display 5351
[E] display 3209

Program:

Code:


LBL E
RAN#
EEX 4
x
INT
1001
X<>Y
X>Y (Test 7)
GTO 2
GTO E
----------------------------
LBL 2
2 ÷
INT
2 x 1 +
STO 0
3
STO 1
-----------------------------
LBL 4
RCL 0
RCL 1
÷
LSTx
X>Y (Test 7)
GTO 0
X<>Y
INT
LSTx
X=Y (Test 5)
GTO E
2
STO+1
GTO 4
-----------------------------
LBL 0
RCL 0
RTN

Remark: This program will run much faster on the emulator.

Gamo

[Dieter]

(10-20-2018 01:08 PM)Gamo Wrote:  This program use Random Number Generator range approximately
from 1000 to 9999 and some time the result will get below 1000.

The program can and will generate numbers as low as 1 (!which is not prime BTW).

But it can be changed to work correctly. You have to adjust the random number generator. Let me know if you're interested in an updated version.

Dieter

[Gamo]

Thanks Dieter

Sure do, I like to see the update version.

Thank You

Gamo

[Thomas Klemm]

We can restrict the search to numbers that aren't divisible by 2, 3 or 5.
Thus we only choose numbers of the sequence: 7, 11, 13, 17, 19, 23, 29, 31, …
On the other hand we only use these numbers to check if the chosen random number is divisible.

Initialisation
Store the following numbers into these registers:

1: 31
2: 29
3: 23
4: 19
5: 17
6: 13
7: 11
8: 7

Code:

▸LBL A
 RAN#  8  ×  1  +  STO I
 RAN#  333  ×  INT  30  ×
 RCL (i)  +  STO 0
 CLx  STO 9
▸LBL 0
 8  STO I
▸LBL 1
 RCL 0
 RCL 9  RCL (i)  +
 ÷
 LSTx  x>y  GTO 2
 R↓  FRAC  x=0  GTO A
 DSE  GTO 1
 30  STO+ 9
 GTO 0
▸LBL 2
 RCL 0  RTN

Drawback: 2, 3 and 5 are never chosen.
Though I haven't tested I assume it runs faster.

Kind regards
Thomas

[Albert Chan]

Hi, Gamo

Trivia: if I am betting for a number (< 10000) with this random prime generator, I'd pick 9587

Closest prime below 9587 is 9551, with a big gap of 36.
All 36 random integers (9552 to 9587) will go for 9587.

With prime table, removing prime gap bias is trivial: pick a random index, then lookup the prime.

(* Mathematica: primes below 10000 (1229 in total), all equally likely *)

Prime[Random[Integer, {1, 1229}]]

Update: without prime table, a simple way is keep picking random number until it is prime.
I think Klemm prime-wheel code do this ...

[John Keith]

(10-21-2018 01:31 PM)Albert Chan Wrote:  (* Mathematica: primes below 10000 (1229 in total), all equally likely *)

Prime[Random[Integer, {1, 1229}]]

Or on the HP Prime (pun intended):

ithprime(rand(1229)+1)

[Dieter]

(10-21-2018 08:45 AM)Thomas Klemm Wrote:  We can restrict the search to numbers that aren't divisible by 2, 3 or 5.
Thus we only choose numbers of the sequence: 7, 11, 13, 17, 19, 23, 29, 31, …
On the other hand we only use these numbers to check if the chosen random number is divisible.

Good idea. But the program does not do what it is supposed to. ;-)

(10-21-2018 08:45 AM)Thomas Klemm Wrote:  Drawback: 2, 3 and 5 are never chosen.

That's not a problem: The goal is generating prime numbers with four digits, i.e. somewhere between 1000 and 9999.

This means that the progam has to be modified. First a random integer between 33 and 332 has to be be chosen, this is multiplied by 30 and finally one of the primes in R1...R8 is added. This returns a "prime candidate" between 997 and 9991. Since 997 is prime this is a possible result outside of the desired range. The only way to correct this seems to be a simple test whether R0 is less than 1000. Or even less than 1009 which is the first prime in this range. ;-)

This means the program would start like this:

Code:

▸LBL A
 RAN#  8  ×  1  +  STO I
 RAN#  300  ×  INT  33 + 30  ×
 RCL (i)  +  STO 0
...

There is another point where the program may not work as intended. The random number can be as large as 0,99999 99999. Multiplying this by 300 with 10-digit precision returns 300 exactly. Thus the largest prime candidate is 10021. Replacing the constant 300 with something like 299,999 sbould fix this.

Dieter

[Dieter]

(10-21-2018 03:22 AM)Gamo Wrote:  Sure do, I like to see the update version.

OK - but this is not as sophisticated as the program Thomas posted. ;-)

For four-digit prime numbers your program requires two changes.
First change the program start this way:

Code:

LBL E
RAN#
4500
x 
INT
2
x
1001
+
STO 0

This generates an odd random number between 1001 and 9999. Take a closer look at how it works: first a random integer between 0 and 4499 is generated. Then this is doubled, which yields an even integer between 0 and 8998. Finally 1001 is added, resulting in an odd random integer between 1001 and 9999.

Note: due to roundoff issues the largest prime candidate can be 10001. But since 10001 is not a prime number it will not be returned.

Finally remove the complete LBL 3 part at the end. This seems to be left from another prime number program. Replace the "GTO 3" line with "GTO E".

After these changes the program will generate four-digit random integers.

Dieter

[Dieter]

(10-21-2018 01:31 PM)Albert Chan Wrote:  Trivia: if I am betting for a number (< 10000) with this random prime generator, I'd pick 9587

Closest prime below 9587 is 9551, with a big gap of 36.
All 36 random integers (9552 to 9587) will go for 9587.

With prime table, removing prime gap bias is trivial: pick a random index, then lookup the prime.

To be honest, I have no idea what you want to say here – ?!?
How is this related to this thread's subject?

(10-21-2018 01:31 PM)Albert Chan Wrote:  Update: without prime table, a simple way is keep picking random number until it is prime.
I think Klemm prime-wheel code do this ...

Sure, this also is the idea behind Gamo's original program. Generate an odd random integer (here within the wrong range, though) and check if it's prime. If not, try another one.

Dieter

[Albert Chan]

(10-21-2018 05:20 PM)Dieter Wrote:  

(10-21-2018 01:31 PM)Albert Chan Wrote:  Update: without prime table, a simple way is keep picking random number until it is prime.
I think Klemm prime-wheel code do this ...

Sure, this also is the idea behind Gamo's original program. Generate an odd random integer (here within the wrong range, though) and check if it's prime. If not, try another one.

My mistake. I was thrown by Gamo code "LBL3 2 STO+0"
I thought it re-use the random number by searching forward ...

(10-21-2018 04:41 PM)Dieter Wrote:  There is another point where the program may not work as intended. The random number can be as large as 0,99999 99999. Multiplying this by 300 with 10-digit precision returns 300 exactly.

Just curious, is it a guess, or does RAN# really get that high ?

FYI, if above were done with binary math, it will not happen.
Example, for IEEE double, 1 ULP below 1.0 = 1 - 1/2^53

Assumed N >= 0, max(N RAN#) = N - N/2^53

If N is pow-of-2, above is exact (no rounding), thus less than N

If N is not pow-of-2, N/2^53 is slightly bigger than 1/2 ULP, thus rounded-up

So, with binary math, max(N RAN#) = N - 1 ULP < N

[Gene]

RAN# is supposed to generate a number 0 < # < 1, but I have never seen 0.000 000 001 or 0.999 999 999. Guess I don't have the patience to do it enough times :-)

[Thomas Klemm]

(10-21-2018 07:02 PM)Albert Chan Wrote:  Just curious, is it a guess, or does RAN# really get that high ?

Does this count?

6.277946006E-1
STO RAN#
RAN#


Cheers
Thomas

[Thomas Klemm]

(10-21-2018 07:47 PM)Gene Wrote:  RAN# is supposed to generate a number 0 < # < 1

Just try this to get 0:

8.603685347E-1
STO RAN#
RAN#


Kind regards
Thomas

[Dieter]

(10-21-2018 07:47 PM)Gene Wrote:  RAN# is supposed to generate a number 0 < # < 1, ...

According to the 15C manual (I assume the 11C behaves the same here) the random number is 0 ≤ r < 1, so it may also be zero. And since the numbers have ten significant digits, the largest possible value is 0,99999 99999.

Dieter

[Gene]

Oops. Missed that 0 <=

:-)

I always think of nine digits for a decimal on these models since they tend to be shown as 0.XXXX etc.

Would a random of 0. 000 000 000 1 be shown in SCI ?

[Albert Chan]

(10-21-2018 04:41 PM)Dieter Wrote:  There is another point where the program may not work as intended. The random number can be as large as 0,99999 99999. Multiplying this by 300 with 10-digit precision returns 300 exactly.

It seems about 50% of the time, N * RAN# = N rounding bug will not occur.

Example: 501 * 0.99999 99999 = 501 - 501/10^10 = 500.9999999 (501 - 1 decimal ULP)

Just like the binary case (post 10), error term slightly bigger than 1/2 decimal ULP, thus rounded-up.

--> If N mantissa fraction M = N/10^ceil(log10(N)), error term = M decimal ULP
--> If M > 1/2, max(N * RAN#) = N - 1 decimal ULP < N

So, a simulated dice = INT(6 * RAN#) + 1 is correct.

But, coin flipping INT(2 * RAN#) may be wrong. (Unless 2 = standing on edge)
Instead, it should remove precision a bit: INT(2 * FRAC(RAN# + 1))

Edit: because of the minus sign in front of error term, rounding rule is half-round-down.
Example: 500 * 0.99999 99999 = 500 - 500/10^10 = 500 - 0 = 500

[Thomas Klemm]

(10-21-2018 10:02 PM)Gene Wrote:  Would a random of 0. 000 000 000 1 be shown in SCI ?

Yes:

0.929424688e-1
STO RAN#
RAN#

[Thomas Klemm]

(10-21-2018 09:09 PM)Dieter Wrote:  According to the 15C manual the random number is 0 ≤ r < 1, so it may also be zero. And since the numbers have ten significant digits, the largest possible value is 0,99999 99999.

The parameters of this PRNG are chosen such that the period is maximal:

(11-06-2014 12:06 AM)Thomas Klemm Wrote:  These numbers fulfill the three requirements of the Hull-Dobell Theorem:

1. c = 1017980433 and m = 10000000000 are relatively prime
2. a - 1 = 1574352260 is divisible by all prime factors of m = 10000000000: 2 and 5
3. a - 1 = 1574352260 = 4*393588065 is a multiple of 4 if m = 10000000000 is a multiple of 4.
   

Thus we can conclude to have maximal period.

Otherwise we couldn't be sure that each possible value is eventually hit.

Cheers
Thomas

[Thomas Klemm]

(10-22-2018 01:02 AM)Albert Chan Wrote:  It seems about 50% of the time, N * RAN# = N rounding bug will not occur.

I wouldn't call this a bug. After all the result is rounded correctly.
It really depends on the mantissa \(m\) of the number.
If \(1 < m \leqslant 5\) then the result is rounded up.
Otherwise it's rounded down.

Obviously I wasn't aware of this fact when I wrote the program.
So, props to Dieter for pointing that out.

A similar thing happens with 0.9999999998.
Here the result is rounded up if \(1 < m \leqslant 2.5\).

And then with 0.9999999997 the result is rounded up if \(1 < m \leqslant 1.666666666\).
Finally with 0.9999999996 the result is rounded up if \(1 < m \leqslant 1.25\).

Cheers
Thomas

[Thomas Klemm]

(10-21-2018 04:41 PM)Dieter Wrote:  That's not a problem: The goal is generating prime numbers with four digits, i.e. somewhere between 1000 and 9999.

Yeah, I should have read the specification more thoroughly.

(10-20-2018 01:08 PM)Gamo Wrote:  This program use Random Number Generator range approximately
from 1000 to 9999 and some time the result will get below 1000.

Well, then I was not so far off.

Cheers
Thomas