Qbasicnews.com
January 22, 2022, 08:34:23 AM *
Welcome, Guest. Please login or register.

Login with username, password and session length
News: Back to Qbasicnews.com | QB Online Help | FAQ | Chat | All Basic Code | QB Knowledge Base
 
   Home   Help Search Login Register  
Pages: [1] 2
  Print  
Author Topic: Prime Contest  (Read 15071 times)
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« on: October 12, 2006, 08:56:19 AM »

Make a program that calculates the value of the 1,000th, the 10,000th and the 100,000th prime. If FB is used the 1,000,000th prime must be calculated too.
The first prime is 2, the second is 3 and so on...
The number of primes up to the integer x can be aproximated by x/(log(x)-1) (QB's log)
The winner is who has the three primes right and faster.

References:
All you wanted to know about primes and never dared to ask http://primes.utm.edu/

The idea is to generate primes and count them up to the required count. There are a lot of optimizations possible in this generation and counting so I hope we will get some interesing sources.

We can open 3 categories according to the speed and memory limitations: Qbasic, Qb4.5 and FreeBASIC

EDITED: Added an additional request for FB, without it even an unoptimized souce takes less than 0.5 second.

EDIT2: For reference, the results are
Code:

the    1000 th prime is      7919
the   10000 th prime is    104729
the  100000 th prime is   1299709
the 1000000 th prime is  15485863

of course your program must FIND these results!
Logged

Antoni
nkk_kan
Member
*
Posts: 73


« Reply #1 on: October 12, 2006, 09:47:36 AM »

err winer? :lol:
Logged

\__/)
(='.'=) Copy bunny into your signature to
(")_(") help him gain world domination.
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #2 on: October 12, 2006, 09:51:30 AM »

Thanks!
Logged

Antoni
yetifoot
Ancient Guru
****
Posts: 575



« Reply #3 on: October 13, 2006, 07:53:32 AM »

Heres my simple FB prime program, using trial divison.  It basically combines the only two rules i know, that theres no need to look at even numbers (except 2), and that you only need to check if n is divisible by any prime numbers smaller than sqr(n)

It takes about 0.5 seconds to get the 100,000 prime, and about 12.5 to get the 1,000,000 on my P4 1.8Ghz

I've been reading the link Antoni posted and may try again once i learn some different methods.

Code:
Dim Shared prime_list(1 To 1000000) As uInteger
Dim Shared prime_list_count As uInteger

Dim As uInteger prime_val, i, j, n
Dim As Integer is_prime

prime_list(1) = 2    ' Setup the first two values
prime_list(2) = 3
prime_list_count = 2
prime_val = 3

Do
  prime_val += 2     ' Add 2, we don't need to look at even numbers
  is_prime = -1
  n = prime_val
  i = 1
  j = Int(sqr(n)) + 1 ' Only need to check if n is divisible by any prime smaller than sqr(n)
  While prime_list(i) < j
    If (n mod prime_list(i)) = 0 Then
      is_prime = 0
      Exit While
    End If
    i += 1
  Wend
  If is_prime Then
    prime_list_count += 1
    prime_list(prime_list_count) = prime_val
  End If
Loop Until prime_list_count = 1000000

Print prime_list(1000)
Print prime_list(10000)
Print prime_list(100000)
Print prime_list(1000000)
Logged

EVEN MEN OF STEEL RUST.
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #4 on: October 13, 2006, 10:41:55 AM »

Great, we have one entry!
But sieves are faster...
Logged

Antoni
yetifoot
Ancient Guru
****
Posts: 575



« Reply #5 on: October 13, 2006, 05:35:41 PM »

Heres another one, a fairly standard sieve, using a bit array.  I tried it just using a regular array which will take much more memory, and i though would be faster, but this bit based method won in the timings.

I also tried to change IsGood to a macro, but that actually slowed it down, I think thats P4 weirdness.

The only real optimization i did for this was the multiples of two, which because i used a bit array, i just iterated through, and masked out all multiple of two numbers.

This one runs in about 3.5 to 4 seconds, so 1/3 of the time of my first attempt, and i'm sure there are better methods yet...

Code:
Const MaxPrime = 1000000             ' The max value for P
Const MaxVal   = MaxPrime * 16       ' The max value for N (a bit of a cheat using * 16...)
Const Max32    = (MaxVal \ 32) + 1   ' The number of 32-bit vars needed to store bit array

Dim Shared BitsArray(0 To Max32) As uInteger

#macro MarkGood(n)
  Scope
    Dim As uInteger p = (n) \ 32, o = (n) mod 32
      BitsArray(p) = BITRESET(BitsArray(p), o)
  End Scope
#endmacro

#macro MarkBad(n)
  Scope
    Dim As uInteger p = (n) \ 32, o = (n) mod 32
      BitsArray(p) = BITSET(BitsArray(p), o)
  End Scope
#endmacro

Function IsGood(ByVal n As uInteger) As Integer
  Dim As uInteger p = n \ 32, o = n mod 32
    Return NOT BIT(BitsArray(p), o)
End Function

Dim As uInteger i, n1, n2, count, mask = &H55555555

For i = 0 To Max32 - 1 ' Mark off all multiples of two quickly using a mask of 10...
  BitsArray(i) = mask
Next i

MarkGood(2) ' Restore 2 as a prime
MarkBad(1)  ' Make 1 not a prime

count = 1 ' start count offset at 1 to account for 2 being a prime

For n1 = 3 To MaxVal
  If IsGood(n1) Then
    count += 1
    If count = 1000 Then Print n1
    If count = 10000 Then Print n1
    If count = 100000 Then Print n1
    If count = 1000000 Then
      Print n1
      Exit For ' We've found the 1 millionth prime, we can quit
    End If
    For n2 = (n1 + n1) To MaxVal Step n1  ' work from n+n to max marking off multiples
      MarkBad(n2)
    Next n2
  End If
Next n1
Logged

EVEN MEN OF STEEL RUST.
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #6 on: October 13, 2006, 05:39:23 PM »

This is an interesting source by Rich Geldreich anyone can find at ABC packets. It uses an original idea and it's probably the only way to do it in QBasic 1.1, because of the 160K memory limits.

It ran in 14 secondes in QB1.1 and in 2 seconds compiled in QB4.5


Code:

'Prime tally using a moving window version of the Erathostenes' Sieve
'Antoni Gual 10/2006 for the comtest at QBN. Qbasic1.1 version  
'----------------------------------------------------------------
'A true bit sieve would be faster, but the memory sizes in QB1.1
'require bold ideas. This one was created for QB by Rich Geldreich in
'1992 from an idea in Donald Knuth's TAOCP.
'
'In a normal sieve each prime found is used in turn to mark all its
'composites thus the complete sieve must be hold in memory tor the final
'tally. In Rich's version all primes found so far are used at the same
'time to mark composites in the same moving slice of ths sieve, the
'numbers left unmarked are primes,and they can be counted as the slice
'progresses.
'In fact there is no data representing the sieve slice...only a priority
'queue that keeps the primes and it's factors used in the present sieve
'slice. This queue has to be dimensioned to hold all primes up to the
'square root of the maximum prime, the present size of 4096 would allow
'for primes up to 2^31.

'Additional optimizations;
'  Multiples of 2 and 3 are skipped
'  A prime p starts to sieve at p*p, because p*a  for a<p will be found
'  by a.
'  The heap is an udt but is kept in separate arrays for speed.


DEFINT A-Z

DECLARE SUB PutPrime (a&)
DECLARE FUNCTION GetPrime& ()

CONST heapsize = 4096

'Priority queue
DIM heapq(1 TO heapsize) AS LONG
DIM HeapQ1(1 TO heapsize) AS LONG
DIM HeapQ2(1 TO heapsize) AS LONG

DIM SHARED n AS LONG
DIM t AS LONG
DIM Q AS LONG, Q1 AS LONG, Q2 AS LONG
DIM TQ AS LONG, TQ1 AS LONG
DIM u AS LONG, primepos AS LONG, cnt AS LONG

primepos = 1000

n = 5
d = 2
r = 1
t = 25
heapq(1) = 25
HeapQ1(1) = 10
HeapQ2(1) = 30

cnt = 2

DO
  DO
    Q = heapq(1)
    Q1 = HeapQ1(1)
    Q2 = HeapQ2(1)

    TQ = Q + Q1
    TQ1 = Q2 - Q1

    '***Insert Heap(1) into priority queue
    i = 1
    DO
        j = i * 2
        IF j <= r THEN
           
            IF j < r THEN
                IF heapq(j) > heapq(j + 1) THEN
                  j = j + 1
                END IF
            END IF

            IF TQ > heapq(j) THEN
                heapq(i) = heapq(j)
                HeapQ1(i) = HeapQ1(j)
                HeapQ2(i) = HeapQ2(j)
                i = j
            ELSE
                EXIT DO
            END IF
        ELSE
            EXIT DO
        END IF
    LOOP
    heapq(i) = TQ
    HeapQ1(i) = TQ1
    HeapQ2(i) = Q2
    '***

  LOOP UNTIL n <= Q

  DO WHILE n < Q
    cnt = cnt + 1
    IF cnt < heapsize THEN heapq(cnt - 2) = n
    IF cnt = primepos THEN
       PRINT USING "The  ####### th prime is ######### "; primepos; n
       IF primepos = 100000 THEN PRINT "Ended": SYSTEM
       primepos = primepos * 10
    END IF
    n = n + d
    d = 6 - d
  LOOP

  IF n = t THEN
    u = heapq(r + 1)
    t = u * u

    '***Find location for new entry
    j = r + 1
    DO
      i = j \ 2
      IF i = 0 THEN
        EXIT DO
      END IF
      IF heapq(i) <= t THEN
        EXIT DO
      END IF
      heapq(j) = heapq(i)
      HeapQ1(j) = HeapQ1(i)
      HeapQ2(j) = HeapQ2(i)
      j = i
    LOOP
    '***
    heapq(j) = t
    IF (u MOD 3) = 2 THEN
      HeapQ1(j) = 2 * u
    ELSE
      HeapQ1(j) = 4 * u
    END IF
    HeapQ2(j) = 6 * u

    r = r + 1
   
  END IF
  n = n + d
  d = 6 - d

LOOP
Logged

Antoni
yetifoot
Ancient Guru
****
Posts: 575



« Reply #7 on: October 13, 2006, 06:18:50 PM »

I've just been looking at some of the previous posts about primes, seeing what methods other people used.  I found a lot by you Antoni!, and some other interesting things I may try and add to my program.

I came across this too which made me laugh

http://members.surfeu.fi/kklaine/primebear.html
Logged

EVEN MEN OF STEEL RUST.
yetifoot
Ancient Guru
****
Posts: 575



« Reply #8 on: October 14, 2006, 12:12:18 AM »

Thats a nice one Antoni, works fast, i'm still trying to understand how it works.

I improved my second one, its a bit faster now, but I still need to learn more to make it go even faster.  Some of the code wasn't necessary, and i even forgot to ignore multiples of 2.

Code:
Const MaxPrime = 1000000             ' The max value for P
Const MaxVal   = MaxPrime * 16       ' The max value for N (a bit of a cheat using * 16...)
Const Max32    = (MaxVal \ 32) + 1   ' The number of 32-bit vars needed to store bit array

Dim Shared BitsArray(0 To Max32) As uInteger
Dim As uInteger n1, n2, count, p, o, n1x2, n1x3
Dim As uInteger steps(0 To 47) = { _
 2,  4,  2,  4,  6,  2,  6,  4,  2,  4,  6,  6,  2,  6,  4,  2, _
 6,  4,  6,  8,  4,  2,  4,  2,  4,  8,  6,  4,  6,  2,  4,  6, _
 2,  6,  6,  4,  2,  4,  6,  2,  6,  4,  2,  4,  2, 10,  2, 10  _
} ' This lookup table is used to calculate the step, to avoid multiples of 2, 3, 5, and 7
  ' any more that that and the table becomes very large (the one inc. 11 is 480 entrys)
Dim As Integer curr_step
Dim As Integer prime_to_find = 1000

count = 4 ' start count offset to account for 2, 3, 5 and 7 being prime

n1 = 11   ' start at 11 due to count starting at 4
While n1 <= MaxVal
  p = n1 shr 5  ' \ 32      'p is integer postition in bitarray
  o = n1 and 31 ' mod 32    'o is bit offset
  If NOT BIT(BitsArray(p), o) Then ' If the bit isn't set then it hasn't been struck out
    count += 1
    If count = prime_to_find Then
      Print Using "###,###,###th prime - ###,###,###"; prime_to_find; n1
      If prime_to_find = 1000000 Then Exit While
      prime_to_find *= 10
    End If
    If n1 <= (sqr(MaxVal) + 1) Then ' Only strike out multiples of primes <= sqr(MaxVal)
      n1x2 = n1 + n1                ' (the +1 is just to account for any rounding, may not
      n1x3 = n1x2 + n1              '  be needed?)
      ' we don't need to step by n1, as that will wastefully look at even numbers (odd+odd=even)
      ' same goes for start pos, n1 is odd, so 2*n1 not needed, start at 3*n1
      For n2 = n1x3 To MaxVal Step n1x2  ' work from 3n to max marking off multiples
        p = n2 shr 5  ' \ 32      'p is integer postition in bitarray
        o = n2 and 31 ' mod 32    'o is bit offset
        BitsArray(p) = BITSET(BitsArray(p), o) ' Set the bit to show its bad
      Next n2
    End If
  End If
  ' we can step by set amounts, to avoid multiples of 2, 3 and 5
  n1 += steps(curr_step)
  curr_step += 1
  If curr_step = 48 Then curr_step = 0
Wend


EDIT: added a check for n1 <= sqr(maxval)
Logged

EVEN MEN OF STEEL RUST.
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #9 on: October 14, 2006, 06:13:17 AM »

Quote from: "yetifoot"

I came across this too which made me laugh

http://members.surfeu.fi/kklaine/primebear.html


Ok,Arktinen Krokotiili Projekti is the winner! Let's close the contest, nothing more more can be done.. :rotfl:  :rotfl:  :rotfl:

EDITED:
Their javascript prime finding algorithm is a little slow..
Code:

function is_x_prime_number(x)
{
var limit=0;
var div=3;
var x_limit = Math.sqrt(x);
while (x%div!=0 && div<x_limit)div+=2;
is_prime = (x%div==0 && x!=div)*1
return is_prime;
}

A simple trial division...
Logged

Antoni
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #10 on: October 14, 2006, 06:50:38 AM »

Now seriously:
Provisional scores for FB:

my (Rich's) entry         Execution time: 2.585 s
yetifoot's second entry Execution time: 4.377 s
yetifoot's second entry Execution time: 2.131 s

Rich Geldreich's method is not that slow after all....
Logged

Antoni
yetifoot
Ancient Guru
****
Posts: 575



« Reply #11 on: October 14, 2006, 09:18:08 AM »

Indeed, I found rich geldreichs method fast, even considering the fact its aimed at qb 1.1 and has to account for that.

I spotted another optimization for my second one, i've just edited the post, as it wasn't enough to warrent posting the whole code again.  I wasn't checking n1 <= sqr(maxval), even though i knew about that trick, i hadn't been able to find where to add it until know.

I think my entry is now as fast as i can get it, but i already thought that about 5 times already, so who knows? Smiley

Is anyone else thinking of entering some code?

Maybe its just me and you Antoni, i fear i'll be in last place if that is the case Smiley
Logged

EVEN MEN OF STEEL RUST.
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #12 on: October 14, 2006, 10:33:49 AM »

Your modified last entry scores 1.27 seconds, it nearly doubles its previous speed...

Im' surprised other people has'nt showed up. Is someone coding something?
Logged

Antoni
Antoni Gual
Na_th_an
*****
Posts: 1434



WWW
« Reply #13 on: October 14, 2006, 04:18:01 PM »

I think i have killed the competition by myself...

First I forgot QB1.1 can't do huge arrays so a sieve up to 1.6M to find the 100,000th prime is out of limits.

Then I posted Rich Geldeich's solution for QB4.5, thinking a sieve was faster. In fact it is..in FreeBASIC. With all that huge arrays, index offseting  and bit twiddling a sieve in QB4.5 is actually slower than Rich's solution.

So perhaps this is the reason because no one else enters.

Here is a sieve for QB 4.5


Code:

'Prime sieve for QB4.5  by Antoni Gual 10/2006
'Run QB4.5 with /AH
 
DEFLNG A-H
DEFINT I-M
amaxp = 1300000
'index offset of -4000 so it does not reach the limit of 16737...
REDIM p(-4000 TO amaxp \ 64 - 4000) AS LONG

'to avoid bit rotations
DIM powers2(31) AS LONG
b = 1
FOR i = 0 TO 31
  powers2(i) = b
  IF i = 30 THEN EXIT FOR
  b = b + b
NEXT
powers2(31) = &H80000000

ctarget = 1000
cnt = 1
asqrt = INT(SQR(amaxp))

FOR b = 3 TO amaxp STEP 2
 IF (p((b \ 64) - 4000) AND powers2((b AND 63) \ 2)) = 0 THEN
   cnt = cnt + 1
   IF cnt = ctarget THEN
     PRINT ctarget; b
     IF ctarget = 100000 THEN SYSTEM
     ctarget = ctarget * 10
   END IF
   IF b <= asqrt THEN
    FOR bb = b * b TO amaxp STEP 2 * b
      p((bb \ 64) - 4000) = p((bb \ 64) - 4000) OR powers2((bb AND 63) \ 2)
    NEXT
   END IF
 END IF
NEXT
Logged

Antoni
yetifoot
Ancient Guru
****
Posts: 575



« Reply #14 on: October 15, 2006, 01:13:51 PM »

There are a few optimizations in there that I hadn't thought of.

1. Using a lookup table for the powers, to make the bit-twiddling faster.

2. Starting the inner loop at n * n, i had discovered that i only needed to start an 3n, but n ^ 2 is even better!

3. Makeing the bitarray \ 64, because theres no need to store data for even numbers, I had thought of this, but I never got it working in my program, that halves the memory overheads.

I added these to my program, but it was only a very very small speed increase (less than 5%), so I haven't updated my code.
Logged

EVEN MEN OF STEEL RUST.
Pages: [1] 2
  Print  
 
Jump to:  

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines Valid XHTML 1.0! Valid CSS!