Prime Number Search

[SlideRule]

A excerpt from Interface Age, 1982-03, pg.12

… Instead of finding the PRIME numbers, it is quicker to search for the numbers that are NOT PRIME. The program to do this is listed here. …
      10 PRINT "Starting."
      20 DEFINT A-Z
      30 DIM A(1000)
      40 FOR I=2 TO 500
      50 FOR J=2 TO 32
      60 IF I*J>1000 GOTO 90
      70 A(I*J)=1
      80 NEXT J
      90 NEXT I
     100 FOR I=1 TO 1000
     110 IF A(I)=0 THEN PRINT I;
     120 NEXT I
     130 PRINT "Finished."
  Although the problem is trivial, the exercise emphasizes the importance of the approach adopted - a fairly simple revision of basic ideas …

BEST!
SlideRule

error corrected

[Maximilian Hohmann]

Hello,

back in the old days of antiquity this method used to be called "the sieve of Erathostenes". Since then (Erathosthenes lived between 275 and 194 BC) this has been part of every programming course for every programming language ever written...

Regards
Max

[Allen]

(07-19-2024 01:40 PM)SlideRule Wrote:  ...it is quicker to search for the numbers that are NOT PRIME....

for \( n > 2^{32} \) this is less true.

Far more interesting rabbit hole than S.o.E on a calculator: https://github.com/bbuhrow/yafu

I think

Code:


100 FOR 1=1 TO 1000
110 IF A(I)=0 THEN PRINT I;
120 NEXT I

should be

Code:


100 FOR I=1 TO 1000         # <----   change int 1 to variable "I"
110 IF A(I)=0 THEN PRINT I;
120 NEXT I

[SlideRule]

Thanks, correction made.

BEST!
SlideRule

[ttw]

I have been using a modified "Wheel" method. It's an extension of Eratsthonese Sieve using irregular skipping. There's a trade-off in the complexity of search and speed. For a binary computer, I use the 2-3-5 (or 30) sieve; All primes are of the form {1, 7, 11, 13, 17, 19, 23, 29)+30*K for all K starting at 0. (The 2, 3, and 5 are done separately). I like to use this on a binary computer there are 8 primes in the wheel so each byte represents 30 numbers. It's not too efficient; 8 primes out of thirty so a 15 to 4 compression and speed. This can be expanded, by including 7 to give 48 out of 210. On the HP50g (EMU48) I usually go one more step to 480 out of 2310 numbers. One doesn't actually scan 2310 numbers but instead generates blocks of 480 which are possible primes, then tests. Most of the long-range searches I have read about went up to 30030 size. (HP50g is fast enough but doesn't have enough memory.)

As I have usually needed primes of restricted classes I have made a few modifications. For Sophie Germain primes (or "safe" primes), or for primes of the form 4k+1, there's a simple trick. Take P=(1,7,11,13,17,19,23,29)*30k and create 2*p+1 or 60k+(2,14,22,26,34,38,46,58)+1 or 60ik+(3,15,23,27,35,39,47,59). Primes of form 30k+(1,7,13,17,19) cannot generate Sophie Germain Primes. (I used 7 and 11 and did a search starting around 2^60 or 2^120.)

I have another, unpublished, method that I'll introduce here. It's fairly efficient but there are some big improvements that I haven't implemented yet.

One generates Stern's Diatomic Series using the two-term recursion: S(n+1)=S(n)+S(n-1)-2*(S(n-1) mod S(n)). (S. Northshield, Stern’s diatomic sequence 0,1,1,2,1,3,2,3,1,4,.... Amer. Math. Monthly 117 (2010, no. 7,581-598.) Starting with S0=1 and S1=1 gives the first terms. However, one can start with any 2 numbers that are relatively prime.

The terms of the Stern Diatomic Sequence taken in pairs (S(n) and S(n+1)) or (S(n-1) and (S(n)) generate every fraction A/B in lowest terms exactly once. The order is according to the sum of the partial fractions of A/B. This means one can quickly generate some prime-based linear-congruential generators with good properties. Two thirds of the generated numbers are odd (as each fraction occurs only once, even/even cannot occur. The safe prime associated with a Sophie Germain prime must be congruent to 11 modulo 12 (or maybe it's the Sophie Germain prime itself) so one can check each generated number for being congruent to 11 mod 12. As the parity pattern is odd-odd-even, one can unroll a loop 3 times and only check odd numbers. Unfortunately, the congruence mod 12 doesn't occur every 12 occurrences as it would in normal ordering of integers. I generated about 10 130 bit Sophie Germain primes in about 10 minutes on a cell phone using this.

Enhancements are available but I haven't worked things out yet. So far, I only use it to generate lots of pseudo-random number generators. The divisibility of the Stern Diatomic Sequence (not to be confused with a Stern Sequence or a Sturm Sequence or Sturm Sequence...) follows a complicated pattern but one not to hard to program.

[John Keith]

This is all very interesting, it deserves to be expanded into an article if you feel like doing so. I am particularly interested in seeing the HP 50 programs that you mention.

(07-19-2024 11:57 PM)ttw Wrote:  One generates Stern's Diatomic Series using the two-term recursion: S(n+1)=S(n)+S(n-1)-2*(S(n-1) mod S(n)).

The recurrence in the OEIS page, a(2*n) = a(n), a(2*n+1) = a(n) + a(n+1) is faster to compute because it doesn't require MOD. Here is a program implementing the recurrence. It will run on any RPL calculator, returning the first 2n+1 terms of the sequence.

Code:


\<< \-> n
  \<< 0 1 1 2 n
    FOR k k PICK k PICK DUP ROT + SWAP
    NEXT n 2 * 1 + \->LIST
  \>>
\>>

[ttw]

That recursion is much faster, but the memory grows linearly. I run about 10,000 steps/minute in my latest ranges (about 112 to 127 bits). Without storage, the depth recursion is somewhere about 140 steps. About half the mods have zero for the quotient which helps. I did have a recursion formula for the mod function but it got a bit complicated. That's why I used the limited context formula.

The other big advantage is that the sum of partial quotients in S(i)/S(i+1) is constant. This gives good parameters for a linear congruential pseudo-random number generator. One doesn't need to search for a "good" multiplier. I need to look further, but the ordering seems to come from taking the partial quotients and combining two (or maybe more) early quotients to get the next term. There may be a useful procedure from this.

It's a nice sequence. The Stern-Brocot tree and the Calkin-Wilf tree can be derived from the sequence. Also Minkowski's ? function. It provides a 1-1 mapping from rationals to integers.

I'll post some of the code later. I cannot get the computers to speak to the EMU48 program for some reason. I am re-creating a high-speed Fortran code for multiple precision integers that I used to use. Then I can run these on a PC (which I tested as at least 1000x the EMU48.) It's just a hobby.