Tek-Tips is the largest IT community on the Internet today!
Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!
-
Congratulations IamaSherpa on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!
prime numbers 1
- Thread starter shanley06
- Start date
Well, I'm not quite sure what you mean by "logarthmic integrals and x/ln x", but here's how I would do it:
A number is prime if all the prime numbers less than the square root of it do not divide into that number.
So I would take the square root of the number you want to test, calculate all the prime numbers below the square root, and then see if any of the prime numbers divide into the test number.
To calculuate all the prime numbers below the square root, I would loop from 3 to the square root and test each number with the modulus of all the numbers below it's square root.
That's sort of confusing, hopefully this code will help:
A number is prime if all the prime numbers less than the square root of it do not divide into that number.
So I would take the square root of the number you want to test, calculate all the prime numbers below the square root, and then see if any of the prime numbers divide into the test number.
To calculuate all the prime numbers below the square root, I would loop from 3 to the square root and test each number with the modulus of all the numbers below it's square root.
That's sort of confusing, hopefully this code will help:
Code:
top:
PRINT "number?"
INPUT number
max = SQR(number)
'loop only to the square root
primeslength = -1
prime = 1
FOR i = 2 TO max
isprime = 1
FOR j = 0 TO primeslength
IF i MOD primes(j) = 0 THEN isprime = 0: EXIT FOR
'test against every prime before it
NEXT j
IF isprime = 1 THEN
primeslength = primeslength + 1
primes(primeslength) = i
'store the prime number
IF number MOD i = 0 THEN prime = 0: EXIT FOR
'test the number against the prime number
END IF
NEXT i
IF prime = 0 THEN
PRINT number; "is not prime"
ELSE
PRINT number; "is prime"
END IF
GOTO top- Thread starter
- #3
Yeah I've tried that too...It works good until I get into the ten-thousands then it starts to take about a minute to calculate and I just don't have that kind of time to sit at my computer for hours on end just to find out if a really big number is prime...I know there is some good formula out there for finding it fast because Mathematica find if a ten-thousand digit number is prime in about a second
qbasicking
Programmer
I think that is the best way to do it, the only way to make it faster would be to run it on a newer machine that has a mathematical co-processor.
Sometimes it takes even supercomputers days to figure out if a number is prime or not.
Sometimes it takes even supercomputers days to figure out if a number is prime or not.
shanley06,
I think you could use pre-made list of prime numbers... That is, compute it once and store in array (or even read from disk at startup).
So your program will just search the list.
>I know there is some good formula out there for finding it fast
I really doubt that it is.
Quote (backwards translation from Russian, sorry):
"All known methods of building prime table is nothing more then variations on Eratosthenes sieve". Etudes for programmers, Charles Wetherell, 1978. I know the book is pretty old - I just happened to read it today.
(Good book. Reads like Agatha Christie for me ;-) !)
I think you could use pre-made list of prime numbers... That is, compute it once and store in array (or even read from disk at startup).
So your program will just search the list.
>I know there is some good formula out there for finding it fast
I really doubt that it is.
Quote (backwards translation from Russian, sorry):
"All known methods of building prime table is nothing more then variations on Eratosthenes sieve". Etudes for programmers, Charles Wetherell, 1978. I know the book is pretty old - I just happened to read it today.
(Good book. Reads like Agatha Christie for me ;-) !)
Search for
prime factorization algorithm
Got:
Reading
got:
"The best known classical method for factoring numbers is currently the number field sieve [Lenstra and Lenstra, 1993]"
Now you can search for "number field sieve"...
The only thing that I think it's for really BIG numbers (hundreds of digits)... So it might be too complex for our Basic case ;-)
prime factorization algorithm
Got:
Reading
got:
"The best known classical method for factoring numbers is currently the number field sieve [Lenstra and Lenstra, 1993]"
Now you can search for "number field sieve"...
The only thing that I think it's for really BIG numbers (hundreds of digits)... So it might be too complex for our Basic case ;-)
-
1
- #7
Finding all primes is different from checking a number for primality.
For the first problem, the only solution is Erathostenes Sieve and it's variants.
For primality check there are many methods, some of them are not exact, they give you a probability of primality. The exact methods relay in high level maths, so they are out of reach for me. The more comprehensive place about primes i know is
www.utm.edu/research/primes/
Antoni
For the first problem, the only solution is Erathostenes Sieve and it's variants.
For primality check there are many methods, some of them are not exact, they give you a probability of primality. The exact methods relay in high level maths, so they are out of reach for me. The more comprehensive place about primes i know is
www.utm.edu/research/primes/
Antoni
Actually, the formula I implemented above is quite quick! If you modify the array and the number variable:
dim primes(10000) as long
dim number as long
and then put a loop around the main code like:
for number = 5000000 to 6000000
...
next number
Modify the code so it only prints when the number is prime
it works fast - several numbers printed per second starting at 5000000, well above ten-thousand.
Impressive, considering it starts over every loop calculating the first prime number! That means it's calculating all prime numbers between 0 and the square root of five million in a fraction of a second. It would run even faster if it was optimized (ie didn't start over each time, didn't test even numbers, etc).
I can tell you didn't try the code, because there's an error in it - it doesn't even work up to 2,000 because I forgot to redimension the array
Of course, the modifications I just mentioned fix that.
Try it before telling me it's slow
dim primes(10000) as long
dim number as long
and then put a loop around the main code like:
for number = 5000000 to 6000000
...
next number
Modify the code so it only prints when the number is prime
it works fast - several numbers printed per second starting at 5000000, well above ten-thousand.
Impressive, considering it starts over every loop calculating the first prime number! That means it's calculating all prime numbers between 0 and the square root of five million in a fraction of a second. It would run even faster if it was optimized (ie didn't start over each time, didn't test even numbers, etc).
I can tell you didn't try the code, because there's an error in it - it doesn't even work up to 2,000 because I forgot to redimension the array
Of course, the modifications I just mentioned fix that.
Try it before telling me it's slow
- Thread starter
- #9
- Thread starter
- #10
Reading this thread, following Agual's link (very interesting), I end up writing a program that compared speed of testing prime numbers. Some code based on post of mgh730 and shanley06.
It finds 1000 primes starting from 5000000.
Here's the results on my machine (some 300 or 500 Mhz, Win98, QB1.1):
Here's the program.
It finds 1000 primes starting from 5000000.
Here's the results on my machine (some 300 or 500 Mhz, Win98, QB1.1):
Code:
Primality tests: test for speed
(seconds for finding 1000 primes starting from 5000000)
4 methods. Please wait...
Trial Division Time: 7.089844
Wheel factorisation - small wheel Time: 3.191406
Wheel factorisation - bigger wheel Time: 3.398438
sieve of Eratosthenes Time: 1.699219Here's the program.
Code:
'Primality tests: test speed
'4 methods (well, two weels go for 1.5 so 3.5 methods)
'FUNCTION primeTrialDivision, requires nothing
'Trial division method, based on code by shanley06
'This is the simplest approach.
'Wheel factorisation. Thanks to link by Agual.
'FUNCTION primeSmallWheel, requires nothing
'Simple program
'FUNCTION primeBiggerWheel, uses some global arrays (initialised below)
'to store wheel data. Program becomes more complicated...
'probably this eats up supposed speed gain.
'FUNCTION primeEratosthenes, uses global array to store already found primes
'based on code by mgh730
DIM i AS INTEGER
'initialisation and code ------------
'wheel data for bigger wheel: 2,3,5
'3 base primes
DATA 3
DATA 2,3,5
'8 spokes
DATA 8
DATA 1, 7, 11, 13, 17, 19, 23, 29
'Wheel variables and initialisation ----------
COMMON SHARED NBase AS INTEGER, NSpokes AS INTEGER, modulo AS INTEGER
'--- COMMON statement for Eratosthenes-based algorithm - BASIC have me put it here
COMMON SHARED primeslength AS INTEGER
READ NBase
DIM SHARED basePrimes(NBase) AS LONG
FOR i = 1 TO NBase: READ basePrimes(i): NEXT i
READ NSpokes
DIM SHARED spokes(NSpokes) AS LONG
FOR i = 1 TO NSpokes: READ spokes(i): NEXT i
modulo = 1
FOR i = 1 TO NBase: modulo = modulo * basePrimes(i): NEXT i
'--- end of wheel initialisation segment
'--- array for primes for Eratosthenes-based algorithm
DIM SHARED primes(2000) AS LONG 'It shows that it used 332 of them
primeslength = 0
primes(0) = 2
'test code starts here ----------------------------
CLS
PRINT "Primality tests: test for speed"
PRINT "(seconds for finding 1000 primes starting from 5000000)"
PRINT "4 methods. Please wait..."
PRINT
DIM number AS LONG
'if you want to see primes list, set wantToSeePrimes to 1
wantToSeePrimes = 0
method$ = "Trial Division"
PRINT method$; TAB(40);
i = 0
t = TIMER
FOR number = 5000000 TO 6000000
IF primeTrialDivision(number) THEN
IF wantToSeePrimes THEN PRINT number
i = i + 1
END IF
IF i = 1000 THEN EXIT FOR
NEXT number
PRINT "Time: "; TIMER - t
method$ = "Wheel factorisation - small wheel"
PRINT method$; TAB(40);
i = 0
t = TIMER
FOR number = 5000000 TO 6000000
IF primeSmallWheel(number) THEN
IF wantToSeePrimes THEN PRINT number
i = i + 1
END IF
IF i = 1000 THEN EXIT FOR
NEXT number
PRINT "Time: "; TIMER - t
method$ = "Wheel factorisation - bigger wheel"
PRINT method$; TAB(40);
i = 0
t = TIMER
FOR number = 5000000 TO 6000000
IF primeBiggerWheel(number) THEN
IF wantToSeePrimes THEN PRINT number
i = i + 1
END IF
IF i = 1000 THEN EXIT FOR
NEXT number
PRINT "Time: "; TIMER - t
method$ = "sieve of Eratosthenes"
PRINT method$; TAB(40);
i = 0
t = TIMER
FOR number = 5000000 TO 6000000
IF primeEratosthenes(number) THEN
IF wantToSeePrimes THEN PRINT number
i = i + 1
END IF
IF i = 1000 THEN EXIT FOR
NEXT number
PRINT "Time: "; TIMER - t
FUNCTION primeBiggerWheel (n AS LONG)
'wheel factorisation (as I got it)
'[URL unfurl="true"]http://primes.utm.edu/glossary/page.php?sort=WheelFactorization[/URL]
'bigger weel
'primes 2, 3, 5
DIM x AS LONG, xx AS LONG, sq AS LONG, i AS INTEGER
sq = SQR(n)
primeBiggerWheel = 1
'first check base primes
FOR i = 1 TO NBase
IF (n MOD basePrimes(i)) = 0 THEN primeBiggerWheel = 0: EXIT FUNCTION
NEXT i
'then go by spoke
'first loop - special (don't want miss 1 or check for it every time)
FOR i = 2 TO NSpokes 'so skip spokes(1), that is 1
IF (n MOD spokes(i)) = 0 THEN primeBiggerWheel = 0: EXIT FUNCTION
NEXT i
'main loop
x = modulo
DO
FOR i = 1 TO NSpokes
xx = x + spokes(i)
IF xx > sq THEN EXIT FUNCTION
IF (n MOD xx) = 0 THEN primeBiggerWheel = 0: EXIT FUNCTION
NEXT i
x = x + modulo
LOOP
END FUNCTION
FUNCTION primeEratosthenes (number AS LONG)
'sieve of Eratosthenes
'based by code by mgh730
DIM max AS LONG, i AS LONG, j AS INTEGER
max = SQR(number)
'loop only to the square root
primeEratosthenes = 1
'first, check by stored numbers
FOR j = 0 TO primeslength
IF number MOD primes(j) = 0 THEN primeEratosthenes = 0: EXIT FUNCTION
NEXT j
'next check the rest
FOR i = primes(primeslength) + 1 TO max
isprime = 1
FOR j = 0 TO primeslength
IF i MOD primes(j) = 0 THEN isprime = 0: EXIT FOR
'test against every prime before it
NEXT j
IF isprime = 1 THEN
primeslength = primeslength + 1
primes(primeslength) = i
'store the prime number
IF number MOD i = 0 THEN primeEratosthenes = 0: EXIT FUNCTION
'test the number against the prime number
END IF
NEXT i
END FUNCTION
FUNCTION primeSmallWheel (n AS LONG)
'wheel factorisation (as I got it)
'[URL unfurl="true"]http://primes.utm.edu/glossary/page.php?sort=WheelFactorization[/URL]
'small weel
'primes 2 and 3. This would give us two spokes: 1 and 5
DIM x AS LONG, sq AS LONG
sq = SQR(n)
primeSmallWheel = 1
'first check base primes
IF (n MOD 2) = 0 THEN primeSmallWheel = 0: EXIT FUNCTION
IF (n MOD 3) = 0 THEN primeSmallWheel = 0: EXIT FUNCTION
'then go by spoke (1st and 5th in a table width 6)
x = 5
DO
IF (n MOD x) = 0 THEN primeSmallWheel = 0: EXIT FUNCTION
x = x + 2
IF (n MOD x) = 0 THEN primeSmallWheel = 0: EXIT FUNCTION
x = x + 4
LOOP UNTIL x > sq
END FUNCTION
FUNCTION primeTrialDivision (n AS LONG)
'based on code by shanley06
DIM x AS LONG
primeTrialDivision = 1
FOR x = 2 TO SQR(n)
IF (n MOD x) = 0 THEN primeTrialDivision = 0: EXIT FUNCTION
NEXT x
END FUNCTIONtsh: A primality test is a program to test if 23453453312452345241123 is prime. Yours are prime listings.
Add this one to your collection...
Antoni
Add this one to your collection...
Code:
WIDTH , 50: CLS
CONST maxlng = &H7FFFFFFF
PRINT "PRIME CALCULATOR BY ANTONI GUAL agual@eic.ictnet.es"
DO
PRINT "Maximum prime(5 to "; maxlng; "): "; : INPUT mp&
LOOP UNTIL mp& >= 5 AND mp& <= maxlng
CONST stor = 4750
DIM pr(stor + 1) AS LONG
DIM pr1(stor + 1) AS LONG
pr(1) = 2 * 3
pr1(1) = 3 * 3
ind& = 1
ind2% = 1
pr1(2) = maxlng
sqre& = 3 * 3
test& = 5
OPEN "nul" FOR OUTPUT AS #1
t! = TIMER
PRINT 2, 3,
CNT% = 2
DO
DO
a& = pr1(1)
IF a& >= test& THEN EXIT DO
b& = pr(1)
a& = a& + b&
i% = 1
DO
j% = i% + i%
IF j% > ind2% THEN EXIT DO
IF pr1(j%) > pr1(j% + 1) THEN j% = j% + 1
IF a& <= pr1(j%) THEN EXIT DO
pr(i%) = pr(j%)
pr1(i%) = pr1(j%)
i% = j%
LOOP
pr1(i%) = a&
pr(i%) = b&
LOOP
IF a& > test& THEN
ind& = ind& + 1
PRINT test&,
'the print #1 "." is just for W2000, you can erase it if you want
CNT% = CNT% + 1: IF CNT% = 250 THEN CLS : CNT% = 0: PRINT #1, ".";
IF ind& <= stor THEN pr(ind&) = test&: IF ind2% = stor THEN EXIT DO
END IF
test& = test& + 2
IF test& > sqre& THEN
GOSUB addnew
IF LEN(INKEY$) THEN EXIT DO
END IF
LOOP UNTIL test& > mp&
PRINT : PRINT : PRINT "time for "; ind& + 1; "primes: "; TIMER - t!; "seconds"
a$ = INPUT$(1)
END
addnew:
ind2% = ind2% + 1
a& = pr(ind2%)
sqre& = a& * a&
j% = ind2%
DO
i% = j% \ 2
IF i% = 0 THEN EXIT DO
IF pr1(i%) <= sqre& THEN EXIT DO
pr(j%) = pr1(i%)
pr1(j%) = pr1(i%)
j% = i%
LOOP
pr1(j%) = sqre&
pr(j%) = a& + a&
pr1(j% + 1) = maxlng
RETURN
viewstack:
CLS : FOR ii% = 0 TO ind2% + 1: PRINT pr(ii%), pr1(ii%): NEXT
PRINT test&
a$ = INPUT$(1)
RETURNAntoni
Antoni,
you said:
>A primality test is a program to test if 23453453312452345241123 is prime. Yours are prime listings.
Nope, they are not.
Prime listing is used just to test speed.
In each loop there is a call to function "prime#method#", that takes one number and returns 0 or 1.
Thanks for the program.
I tried your prigram; I saw it run.
But I have no idea how it works - I haven't found MOD or "/" - just one "\", and it was "\2".
Indeed!
Is there some complicated math behind it?
Wonderful!
you said:
>A primality test is a program to test if 23453453312452345241123 is prime. Yours are prime listings.
Nope, they are not.
Prime listing is used just to test speed.
In each loop there is a call to function "prime#method#", that takes one number and returns 0 or 1.
Thanks for the program.
I tried your prigram; I saw it run.
But I have no idea how it works - I haven't found MOD or "/" - just one "\", and it was "\2".
Indeed!
Is there some complicated math behind it?
Wonderful!
No complicated math, just a way of building Erathostenes sieve without actually saving it. I keep a priority queue of prime multiples, which I increment by just adding the prime. All gaps the sequence leaves are primes. I't an old method, I think it is in Knuth's "Art of Computer Programming" .
Antoni
Antoni
- Thread starter
- #15
- Status