Page 2 of 3 FirstFirst 123 LastLast
Results 11 to 20 of 23

Thread: Factoring 1

  1. #11
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    I've got it down to 5m49s now (for the primes up to 10 million) and I'm guessing that's it for straight-up thinBasic.

    I removed the need for an SQR in the main calculating loop, replacing it by an addition, a multiplication and a compare (the bit with srt and prd). This is also very useful if (when?) I convert the main loop to assembler. After having done some 330 million numbers the checking runs at some 6,000 numbers per second. Almost another 4 billion numbers to do...

    The program is more complicated than it needs to be because it works with a variable number of assumed primes (variable nap). Be careful changing that number of assumed primes because if you do it will delete any work done with another number of assumed primes. Decoding the primes file with a higher number of assumed primes is slower, but you win a lot of memory.

    I'm working on an include script that can use/process the file generated with this program. Once that's finished I'll put it in a topic in the math subforum. I will definitely complete the entire 32-bits primes file, possibly with the main bit in assembler. I'll make the primes file (102,976,561 bytes, 98.2 MB) available through P2P (.torrent) then.
    '----------------------------------------------------------------------------
    ' !All32BitPrimes
    '
    ' Create a file with all 32-bit prime numbers.
    '
    ' ©2011 Johannes
    '----------------------------------------------------------------------------
     
    Uses "Console","File"
     
    ' Number of assumed primes. Must be between 3 and 6.
    Word nap = 6
     
    ' File name for prime number file.
    String fn_primes = APP_SourcePath+"Primes.primes"
     
    ' Auto-save interval in seconds.
    Quad svi = 60
     
    ' Some initialisations.
    DoEvents(Off)
    PrintL "DETERMINE ALL 32-BIT PRIMES"
    PrintL "==========================="
    DWord i,j,k,n,nbr
    Word srt,prd
    Quad che,mxi = 4294967296
    HiResTimer_Init
    nap=Min(Max(nap,3),6)
     
    ' Check work done so far.
    If Not FILE_Exists(fn_primes) Then FILE_Save(fn_primes,Chr$(nap))
    Integer hnd=FILE_Open(fn_primes,"BINARY")
    String chk=FILE_Get(hnd,1)
    FILE_Close(hnd)
    If chk<>Chr$(nap) Then FILE_Save(fn_primes,Chr$(nap))
     
    ' Determine all 16-bit primes using Euler's Sieve.
    Print "Determining all 16-bit prime numbers: "
    Word sieve(65535)
    For i=2 To 65535
      sieve(i)=i
    Next i
    n=65534
    For i=2 To 255
      If sieve(i) Then
        For j=Int(65535/i) To i Step -1
          If sieve(j) Then sieve(i*j)=0:n-=1
        Next j
      EndIf
    Next i
    Word p16b(n)
    j=0
    For i=2 To 65535
      If sieve(i) Then
        j+=1
        p16b(j)=sieve(i)
      EndIf
    Next i
    ReDim sieve()
    PrintL TStr$(n)+" primes found."
     
    ' Print assumed primes and prime storage information.
    Print "Working with"+Str$(nap)+" assumed primes:"
    Ext fract = 100.0
    For i=1 To nap
      j=p16b(i)
      fract*=(j-1)/j
      Print Str$(j)
    Next i
    PrintL "."
    PrintL "Fraction of stored numbers is "+Format$(fract,"0.000")+"%."
     
    ' Determine sizes of number block and bit mask.
    Word blk,bts,msk = 1
    If nap Then
      For i=1 To nap
        blk*=p16b(i)
        bts*=p16b(i)-1
      Next i
    EndIf
    msk=bts/8
     
    ' Check number of already considered numbers.
    Quad hgh = (FILE_Size(fn_primes)-1)*blk/msk
    fract=Min(100,hgh/42949672.96)
    Print "Every":If msk>1 Then Print Str$(msk)
    Print " byte":If msk>1 Then Print "s"
    Print " define":If msk=1 Then Print "s"
    PrintL Str$(blk)+" numbers."
    PrintL "Total file size for all 32-bit primes is"+Str$(msk*Ceil(mxi/blk)+1)+" bytes."
    PrintL "Currently"+Str$(hgh)+" numbers have been considered."
    PrintL "Overall progress is "+Format$(fract,"0.000")+"%."
     
    ' Check if all work has been done.
    If hgh>=mxi Then
      PrintL "All 32-bit primes have been determined."
      WaitKey
      Stop
    EndIf
     
    ' Construct offset table.
    Long ofs(bts)
    Boolean prm
    k=0
    PrintL "Creating offset table ("+TStr$(bts)+" entries)."
    For i=1 To blk
      prm=TRUE
      For j=1 To nap
        If Mod(i,p16b(j))=0 Then prm=FALSE:Exit For
      Next j
      If prm Then k+=1:ofs(k)=i
    Next i
     
    ' Reserve (a whole lotta) memory for primes list and load existing data.
    Byte dta(msk*Ceil(mxi/blk))
    If FILE_Size(fn_primes)>1 Then Poke$ VarPtr(dta(1)),Mid$(FILE_Load(fn_primes),2)
     
    ' Append 16-bit primes if necessary.
    Byte byt
    If hgh<65536 Then
      PrintL "Appending 16-bit primes."
      For i=1 To msk*Ceil(65536/blk)
        dta(i)=0
      Next i
      For i=1 To n
        j=Array Scan ofs, =Mod(p16b(i),blk)
        If j Then
          j-=1
          k=msk*Int(p16b(i)/blk)+Int(j/8)+1
          j=Mod(j,8)
          byt=1
          SHIFT LEFT byt,j
          dta(k)=(dta(k) Or byt)
        EndIf
      Next i 
      hgh=blk*Int(65536/blk)
    EndIf
     
    ' Main loop for checking successive numbers.
    PrintL "Autosaving every"+Str$(svi)+" seconds."
    svi*=1000000
    Print "Determining further primes"
    k=msk*hgh/blk+1
    byt=1
    Quad ctr
    srt=Int(Sqr(hgh))
    HiResTimer_Get
    While hgh<mxi
      For i=1 To bts
        che=hgh+ofs(i)
        If che>=mxi Then Exit For
        nbr=che
        If srt<65535 Then
          prd=srt+1
          If prd*prd<=nbr Then srt=prd
        EndIf
        prm=TRUE
        j=1
        While p16b(j)<=srt And j<=n
          If Mod(nbr,p16b(j))=0 Then prm=FALSE:Exit While
          j+=1 
        Wend
        If prm Then dta(k)=(dta(k) Or byt)
        SHIFT LEFT byt,1
        If byt=0 Then byt=1:k+=1
      Next i
      hgh+=blk
      ctr+=HiResTimer_Delta
      If ctr>=svi Then ctr=0:SavePrimes():Print "."
    Wend
    PrintL
     
    ' Save primes list at program end.
    PrintL "Saving primes list."
    SavePrimes()
    PrintL "All 32-bits primes have been determined."
     
    DoEvents(On)
    WaitKey
     
    Sub SavePrimes()
      FILE_Save(fn_primes,Chr$(nap)+Peek$(VarPtr(dta(1)),msk*Int(hgh/blk)))
    End Sub
    
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  2. #12
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    I worked on something, just to see how fast I could find all of the primes through 10,000,019.

    I used a vector, called, "primes-vector", with 10,000,020 elements (the extra element is because Racket indexes vectors from 0).

    The procedure, "make-primes-vector", sets the value of an index to #t if the index is prime, and to #f if the index is composite.

    Below, in the REPL interactions, I timed the procedure (in seconds), and then counted the number of primes, as a check.

    (Now am too tired for anything else.)



    ; code --------------------------------------------------------------------------------------
    
    #lang racket
    
    (define primes-vector 0)
    (define max-val 0)
    
    (define (make-primes-vector n)
      (set! max-val n)
      (set! primes-vector (make-vector (+ max-val 1) #t))
      (vector-set! primes-vector 0 #f)
      (vector-set! primes-vector 1 #f)
      (define max-prime (floor (sqrt n)))  
      (define (do-prime a)
        (cond
          ((> a max-prime) void)
          (else
           (remove-multiples a 2)
           (do-prime (next-prime (+ a 1))))))  
      (define (remove-multiples b c)
        (define index (* b c))
        (cond
          ((> index max-val) void)
          (else
           (vector-set! primes-vector index #f)
           (remove-multiples b (+ c 1)))))  
      (define (next-prime d)
        (if (vector-ref primes-vector d) d
            (next-prime (+ d 1))))
      (do-prime 2))
    
    (define (count-primes)
      (define count 0)
      (define (iter-count i)
        (cond
          ((> i max-val) count)
          (else
           (when (vector-ref primes-vector i) (set! count (+ count 1)))
           (iter-count (+ i 1)))))
      (iter-count 2))
    
    (define (time-make-primes-vector n)
      (let ((t1 (current-milliseconds)))
        (make-primes-vector n)
        (let ((t2 (current-milliseconds)))
          (let ((tt (/ (- t2 t1) 1000.0)))
            tt))))
    
    ; REPL interactions -------------------------------------------------------------------------
    
    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (time-make-primes-vector 10000019)
    8.861
    > (count-primes)
    664580
    >
    
    Last edited by danbaron; 05-05-2011 at 11:05.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  3. #13
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by danbaron View Post
    I worked on something, just to see how fast I could find all of the primes through 10,000,019.
    Ah! That's a different kettle of fish.

    You're doing a sieve, presupposing 10,000,019 as a limit. In thinBasic it took me 11.88 seconds and that includes keeping a tally of the number of primes (also 664,580 for me). Please wait while I do this in PowerBASIC. (Racket also compiles, so it's only fair.)
    '----------------------------------------------------------------------------
    ' !Primes10M
    '
    ' Calculate primes up to 10 million.
    '
    ' ©2011 Johannes
    '----------------------------------------------------------------------------
     
    DoEvents(Off)
    HiResTimer_Init
    HiResTimer_Get
     
    ' Determine primes using Euler's Sieve.
    Long i,j,n,x=10000019
    Byte s(n)=1
    s(1)=0:x-=1
    For i=2 To Int(Sqr(n))
      If s(i) Then
        For j=Int(n/i) To i Step -1
          If s(j) Then s(i*j)=0:x-=1
        Next j            
      EndIf
    Next i
     
    Quad t=HiResTimer_Delta
    DoEvents(On)
     
    MsgBox 0,"Time:"+Str$(t/1000000)+$CRLF+"Primes:"+Str$(x)
    
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  4. #14
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by Johannes View Post
    Please wait while I do this in PowerBASIC.
    This is ridiculous.
    '----------------------------------------------------------------------------
    ' !Primes10M
    '
    ' Calculate primes up to 10 million.
    '
    ' ©2011 Johannes
    '----------------------------------------------------------------------------
    #COMPILE EXE
    #DIM NONE
     
    FUNCTION PBMAIN()
    DIM i AS LONG
    DIM j AS LONG
    DIM n AS LONG
    DIM t1 AS LONG
    DIM t2 AS LONG
    DIM x AS LONG
     
    t1=TIMER
     
    ' Determine primes using Euler's Sieve.
    n=10000019
    x=10000019
    DIM s(n) AS BYTE
    FOR i=0 TO n
      s(i)=1
    NEXT i
    s(0)=0
    s(1)=0:x=x-1
    FOR i=2 TO INT(SQR(n))
      IF s(i) THEN
        FOR j=INT(n/i) TO i STEP -1
          IF s(j) THEN s(i*j)=0:x=x-1
        NEXT j
      END IF
    NEXT i
     
    t2=TIMER
     
    MSGBOX "Time:"+STR$(t2-t1)+$CRLF+"Primes:"+STR$(x)
     
    END FUNCTION
    
    Please note that I do everything between the timing: creating the array, filling it with ones, keeping track of the number of primes.

    PowerBASIC times in whole seconds. For those 664,580 primes the time is 0 (zero) seconds. I upped the limit to 100 million and the time was 1 second (5,761,455 primes). Then I upped the limit to 1 billion and I got a significant time: 11 seconds (50,847,534 primes).
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  5. #15
    thinBasic author ErosOlmi's Avatar
    Join Date
    Sep 2004
    Location
    Milan - Italy
    Age
    57
    Posts
    8,777
    Rep Power
    10
    Just in case, for those lazy like me, thinBasic has native IsPrime function

    HiResTimer_Init
    HiResTimer_Get
    '---Determine primes
    Long x      = 10000019 
    Long Count      
    Long nPrimes
    For Count = 1 To x
      If IsPrime(Count) Then Incr nPrimes
    Next
    Quad t=HiResTimer_Delta
    MsgBox 0, "Time:" + Str$(t/1000000) + $CRLF + "Primes:" + Str$(nPrimes)
    
    www.thinbasic.com | www.thinbasic.com/community/ | help.thinbasic.com
    Windows 10 Pro for Workstations 64bit - 32 GB - Intel(R) Xeon(R) W-10855M CPU @ 2.80GHz - NVIDIA Quadro RTX 3000

  6. #16
    thinBasic author ErosOlmi's Avatar
    Join Date
    Sep 2004
    Location
    Milan - Italy
    Age
    57
    Posts
    8,777
    Rep Power
    10
    To me, PowerBasic as always been the best native compiler, ever.
    Not just because it is fast, creates optimized code and is reliable, but mainly because they introduced, during the years, many new commands always keeping very high compatibility with previous code.
    www.thinbasic.com | www.thinbasic.com/community/ | help.thinbasic.com
    Windows 10 Pro for Workstations 64bit - 32 GB - Intel(R) Xeon(R) W-10855M CPU @ 2.80GHz - NVIDIA Quadro RTX 3000

  7. #17
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by Johannes View Post
    After having done some 330 million numbers the checking runs at some 6,000 numbers per second. Almost another 4 billion numbers to do...
    By now my program has determined the primes up to 700 million, one-sixth of the entire work load. Total running time up to now somewhere around 20 or 30 hours. The checking slows down proportional to SQR(n) of course.

    And then I wrote the sieve program in thinBasic and then PowerBASIC (see above). And I realised that if I used bits instead of bytes I could fit those 4 billion numbers in memory. Which makes a sieve possible.
    '----------------------------------------------------------------------------
    ' !All32BitPrimes
    '
    ' Create a file containing all 32-bit primes using Euler's Sieve.
    '
    ' ©2011 Johannes
    '----------------------------------------------------------------------------
     
    #COMPILE EXE
    #DIM NONE
     
    FUNCTION PBMAIN()
      DIM i AS DWORD
      DIM j AS DWORD
      DIM k AS DWORD
      DIM l AS DWORD
      DIM maxb AS DWORD
      DIM maxn AS DWORD
      DIM q AS QUAD
      DIM t1 AS DWORD
      DIM t2 AS DWORD
      DIM t3 AS DWORD
      DIM x AS DWORD
     
      ' Upper limit for prime generation.
      maxn = 4294967295
    '  maxn = 1073741823
    '  maxn = 268435455
    '  maxn = 67108863
    '  maxn = 16777215
    '  maxn = 4194303
    '  maxn = 1048575
    '  maxn = 262143
    '  maxn = 65535
     
      ' Define assumed primes.
      DIM ap(5) AS DWORD
      ARRAY ASSIGN ap() = 2,3,5,7,11,13
     
      ' Determine storage parameters.
      DIM blk AS DWORD
      DIM bts AS DWORD
      DIM msk AS DWORD
      blk = 1
      bts = 1
      FOR i=0 TO 5
        blk = blk*ap(i)
        bts = bts*(ap(i)-1)
      NEXT i
      msk = bts/8
     
      ' Construct offset table.
      DIM ofs(bts-1) AS DWORD
      k = 0
      FOR i=1 TO blk
        l = 1
        FOR j=0 TO 5
          IF (i MOD ap(j))=0 THEN l = 0 : EXIT FOR
        NEXT j
        IF l<>0 THEN ofs(k) = i : k=k+1
      NEXT i
     
      ' Create array for sieve.
      q = maxn
      q = blk*INT((q+blk)/blk)
      q = INT((q+7)/8)
      maxb = q
      DIM sieve(maxb-1) AS BYTE
      q = maxn
      INCR q
      SHIFT RIGHT q,3
      FOR i=0 TO q-1
        sieve(i) = &hFF
      NEXT i
     
      t1 = TIMER
     
      ' Zero and one are not prime numbers.
      sieve(0)=sieve(0) AND &hFC
      x = maxn-1
      ' Determine primes using Euler's Sieve.
      FOR i=0 TO INT(SQR(maxn))
        k = i
        SHIFT RIGHT k,3
        l = 1
        SHIFT LEFT l,(i AND 7)
        IF (sieve(k) AND l)<>0 THEN
          j = INT(maxn/i)
          WHILE j>=i
            k = j
            SHIFT RIGHT k,3
            l = 1
            SHIFT LEFT l,(j AND 7)
            IF (sieve(k) AND l)<>0 THEN
              k = i*j
              SHIFT RIGHT k,3
              l = 1
              SHIFT LEFT l,(i*j AND 7)
              sieve(k) = sieve(k) XOR l
              x = x-1
            END IF
            j = j-1
          WEND
        END IF
      NEXT i
     
      t2 = TIMER
     
      ' Extract relevant bits for output file.
      DIM i1 AS QUAD
      DIM k1 AS QUAD
      DIM m AS DWORD
      DIM src AS QUAD
      DIM tgt AS DWORD
      q = maxn
      q = msk*INT((q+blk)/blk)
      DIM dta(q-1) AS BYTE
      i1 = 0
      tgt = 0
      m = 1
      WHILE i1<maxn
        FOR j=0 TO bts-1
          src = i1+ofs(j)
          k1 = src
          SHIFT RIGHT k1,3
          l = 1
          SHIFT LEFT l,(src AND 7)
          IF (sieve(k1) AND l)<>0 THEN
            dta(tgt) = dta(tgt) OR m
          END IF
          SHIFT LEFT m,1
          IF m=256 THEN m = 1 : tgt = tgt+1
        NEXT j
        i1 = i1+blk
      WEND
     
      t3 = TIMER
     
      ' Display statistics.
      MSGBOX "Maximum:"+STR$(maxn)+$CRLF+"Sieve time:"+STR$(t2-t1)+$CRLF+"Primes:"+STR$(x)+$CRLF+$CRLF+"Save time:"+STR$(t3-t2)+$CRLF
     
      ' Save file.
      OPEN CURDIR$+"/PB_Primes.primes" FOR BINARY ACCESS WRITE AS #13
      PUT$ #13,CHR$(6)
      PUT$ #13,PEEK$(VARPTR(dta(0)),q)
      SETEOF #13
      CLOSE #13
     
    END FUNCTION
    
    A little over 5 million numbers per second, and the program is linear. So where the old program takes twice as long to proces a number four times as big, this new program wades through 5 million checks per second if the number is 1 or if it's 1 billion.

    Generating the file for 268 million numbers took 52 seconds. The conversion between the two massive arrays took another 3 seconds. That transfer is bit by bit and I didn't really optimise it. Although only set bits (primes) are transferred, we're still talking about 14,630,843 primes (around 4.8 million bit-by-bit transfers per second).

    Projected time to determine all 4 billion primes: 14 minutes.

    Just in time: my program finished for primes up to one billion-and-a-bit. A total of 54,400,028 primes, determined in 3m34s. Another 12 seconds for convsersion to the output format, resulting file is 24.5 MB. And my (slow) thinBasic program that checks the compressed prime files says it's all in order.
    Last edited by Johannes; 06-05-2011 at 08:49. Reason: Bugs in source
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  8. #18
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    all_of_them.jpg

    And the primes check out. All 203,280,221 of them. Times are in seconds.

    I'll whip up a thinBasic include script that can use that 98.2-MB file. IsPrime, NextPrime, PrevPrime, NthPrime, perhaps some other stuff I can think of. When that's finished I'll put that script along with the compiled PowerBASIC program (less than 10 kB in an archive) in the math subforum.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  9. #19
    thinBasic author ErosOlmi's Avatar
    Join Date
    Sep 2004
    Location
    Milan - Italy
    Age
    57
    Posts
    8,777
    Rep Power
    10
    Consider that, if you like, I can add those keywords as thinBasic native commands in Core engine.

    Thanks a lot
    Eros
    www.thinbasic.com | www.thinbasic.com/community/ | help.thinbasic.com
    Windows 10 Pro for Workstations 64bit - 32 GB - Intel(R) Xeon(R) W-10855M CPU @ 2.80GHz - NVIDIA Quadro RTX 3000

  10. #20
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    Quote Originally Posted by Johannes View Post
    You're doing a sieve, presupposing 10,000,019 as a limit. In thinBasic it took me 11.88 seconds and that includes keeping a tally of the number of primes (also 664,580 for me). Please wait while I do this in PowerBASIC. (Racket also compiles, so it's only fair.)
    If you mean "fair" in the sense of a competition, then, I think none of this is fair. Thinbasic is an interpreter, Racket is compiled into byte code, and, as far as I know, Powerbasic is compiled into machine code; and, we didn't create any of them. Additionally, the timings are for different machines. Additionally, we didn't create the algorithms we used. Etc. (My guess is that, "keeping a tally of the number of primes", does not slow down the procedure very much.)

    And, I question whether speed itself, is self-evidently an admirable and worthy goal. I can make an analogy to driving a car. What good is it being the first one to reach a destination, if you have nothing to do when you get there? Or, I guess, equivalently, what good is it being the first one, to reach a random destination?

    Last edited by danbaron; 06-05-2011 at 08:19.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

Page 2 of 3 FirstFirst 123 LastLast

Similar Threads

  1. Haskell: factoring into primes
    By danbaron in forum Science
    Replies: 0
    Last Post: 12-06-2011, 08:47
  2. Factoring 2
    By danbaron in forum Science
    Replies: 0
    Last Post: 09-05-2011, 06:05

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •