Page 1 of 2 12 LastLast
Results 1 to 10 of 14

Thread: Riemann zeta function in O2

  1. #1

    Riemann zeta function in O2

    also works for fractional arguments.
    <edit> just timed the zeta function in the gsl library and it's almost three times faster than this implementation, bummer.
    some timings


    O2 858 nSec
    gsl 343 nSec
    FreeBASIC 378 nSec


    my routine compiled with FreeBASIC is only 10% slower than gsl, that's not bad.</edit>

     
    Uses "Console" , "oxygen"
    Dim i , j, t As Long 
    Dim scr As String  
    Dim As Long pZeta
    Dim As Long pFinish 
    Dim x, y As Double
    scr = "function zeta(byval n as double) as double link #pZeta
        'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
        'first we sum the first ten terms of the Zeta series,
        'then use the Euler-Maclaurin summation formula.
        dim as double c(7)
        dim as integer i, k, k1, nc
        dim as double s, dx, n1
        c(1) =  8.3333333333333330e-02 'B(2)/2!
        c(2) = -1.3888888888888890e-03 'B(4)/4!
        c(3) =  3.3068783068783070e-05 'B(6)/6!
        c(4) = -8.2671957671957670e-07 'B(8)/8!
        c(5) =  2.0876756987868100e-08 'B(10)/10!
        c(6) = -5.2841901386874930e-10 'B(12)/12!
        c(7) =  1.3382536530684680e-11 'B(14)/14!
        nc=10
        s=0
        for i=1 to nc
          s=s+1/i^n
        next
                           'f(k)=(k+10)^-n
        s=s+nc^(1-n)/(n-1) 'int(f(k)dk,k=0..infinity)=10^(1-n)/(n-1)
        dx=-n*nc^(-n-1)    'first derivative of f(k)
        s=s-.5*nc^(-n)     'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
        n1=nc^(-2)
        k=1
        k1=2
        for i=1 to 7
          s=s-dx*c(i)       'next term of the series
          dx=dx*n1*(n+k1)*(n+k) 'next derivative of f(k)
          k=k+2
          k1=k1+2
        next
        function = s
    end function
    sub finish() link #pFinish
      terminate
    end sub
    "
    O2_Asm scr
    If O2_Error Then
      MsgBox 0,O2_Error
      Stop
    Else
      O2_Exec
    End If
    Declare Function Zeta(ByVal Double ) As Double At pZeta
    Declare Sub Finish() At pFinish
    For j=2 To 20
      x=Zeta(j)
      Console_WriteLine j, x
    Next 
    t=GetTickCount
    For i=1 To 10000
      For j=2 To 50
        x=Zeta(j)
        'Console_WriteLine x
      Next
    Next
    t=GetTickCount-t
    Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" nSec."
    Console_WriteLine "All done. Press any key to finish"
    Console_WaitKey
    Finish()
    
    Last edited by jack; 28-12-2013 at 00:45.

  2. #2
    thinBasic author ErosOlmi's Avatar
    Join Date
    Sep 2004
    Location
    Milan - Italy
    Age
    57
    Posts
    8,777
    Rep Power
    10
    I think 1/4 of the time is taken by thinBasic just for interpreting and calling O2 function in the loop:
    For i=1 To 10000
      For j=2 To 50
        x=Zeta(j)
      Next
    Next
    

    Below code is 100% O2 code and it takes about 600 ms instead of 800 ms
    I moved constant declarations outside zeta function (they are evaluated just once during O2 just in time compilation.
    Created a test lRun function to see how much time was taken by thinBasic for interpretation in the main loop.

    Uses "Console" , "oxygen"
    Dim i , j, t As Long
    Dim scr As String 
    Dim As Long pZeta
    Dim As Long pRun
    Dim As Long pFinish 
    Dim x, y As Double
    
    
    '---O2 code
    scr = "
    
    
        Dim As Double c(7)
        c(1) =  8.3333333333333330e-02 'B(2)/2!
        c(2) = -1.3888888888888890e-03 'B(4)/4!
        c(3) =  3.3068783068783070e-05 'B(6)/6!
        c(4) = -8.2671957671957670e-07 'B(8)/8!
        c(5) =  2.0876756987868100e-08 'B(10)/10!
        c(6) = -5.2841901386874930e-10 'B(12)/12!
        c(7) =  1.3382536530684680e-11 'B(14)/14!
    
    
    function zeta(byval n as double) as double link #pZeta
        'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
        'first we sum the first ten terms of the Zeta series,
        'then use the Euler-Maclaurin summation formula.
        Dim As Integer i, k, k1, nc
        Dim As Double s, dx, n1
        
        nc=10
        s=0
        For i=1 To nc
          s=s+1/i^n
        Next
                           'f(k)=(k+10)^-n
        s=s+nc^(1-n)/(n-1) 'int(f(k)dk,k=0..infinity)=10^(1-n)/(n-1)
        dx=-n*nc^(-n-1)    'first derivative of f(k)
        s=s-.5*nc^(-n)     'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
        n1=nc^(-2)
        k=1
        k1=2
        For i=1 To 7
          s=s-dx*c(i)       'next term of the series
          dx=dx*n1*(n+k1)*(n+k) 'next derivative of f(k)
          k=k+2
          k1=k1+2
        Next
        Function = s
    End Function
    
    
    function lRun() as long link #pRun
      dim i as long
      dim j as long
      dim x as double
      For i=1 To 10000
        For j=2 To 50
          x=Zeta(j)
        Next
    Next
    end function
    
    
    Sub finish() link #pFinish
      terminate
    End Sub
    "
    
    
    '---JIT compile O2 code
    O2_Asm scr
    If O2_Error Then
      MsgBox 0,O2_Error
      Stop
    Else
      O2_Exec
    End If
    
    
    '---Start thinBasic code
    Declare Function Zeta(ByVal Double ) As Double At pZeta
    Declare Function lRun() As Long At pRun
    Declare Sub Finish() At pFinish
    For j=2 To 20
      x=Zeta(j)
      Console_WriteLine j, x
    Next
    
    
    t=GetTickCount
    lRun()
    'For i=1 To 10000
    '  For j=2 To 50
    '    x=Zeta(j)
    '  Next
    'Next
    t=GetTickCount-t
    
    
    Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
    Console_WriteLine "All done. Press any key to finish"
    Console_WaitKey
    Finish()
    
    Do you know where I can get a pre-compiled version (DLL) of GSL library?

    Thanks a lot
    Eros
    Last edited by ErosOlmi; 30-12-2013 at 11:04.
    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

  3. #3
    1. Smart optimizing static C compilers like GCC used to compile GSL may unroll user loops transparently in full or in batches. This is done in order to reduce the overhead and resources involved in reproducing such high-level constructs in low-level machine code they are resolved to by the compiler. This may apply to both the inner loop of Zeta function proper and the outer loop of test script.

    2. Such compilers may also transparently optimize the low-level code to use a faster SSE floating-point instruction set if such a capability is detected automatically. Modern Pentium and Athlon CPU's do support this instruction set in some form or other.

    These two points may account for why original GSL performs faster than FreeBasic whose compiler is somewhat "dumber". JIT compilers like O2 (and DynC, for that matter ) won't perform any such optimizations at all for the sake of much, much faster compilation speed.

    Quote Originally Posted by ErosOlmi View Post
    Do you know where I can get a pre-compiled version (DLL) of GSL library?
    Here's everything you may need to get a 32-bit Windows binary of GSL, the GNU scientific library for numerical computing, working for you.

    Last edited by mike lobanovsky; 30-12-2013 at 22:29.
    Mike
    (3.6GHz i5 Core Quad w/ 16GB RAM, nVidia GTX 1060Ti w/ 6GB VRAM, x64 Windows 7 Ultimate Sp1)

  4. #4
    thank you Eros and mike for your comments, Eros if you have MinGW installed on your system then you can compile the GSL sources no problem but you may want to run the strip --unneeded command afterwards to reduce the size of the dll.
    <off topic> one of my favorite tricks is to
    convert a static library to a dll by first running ar x libname.a and then gcc -shared -static -o libname.dll *.o </off topic>

    I have attached the latest version of the GSL dll.

    forgot post the gsl test code

    
    
    
    'gsl test     
    Uses "Math", "Console" 
    
    
    Type gsl_sf_result
      Vl As Double
      Er As Double
    End Type
    
    
    Declare Function gamma CDECL Lib "gsl.dll" Alias "gsl_sf_gamma" ( ByVal x As Double) As Double 
    Declare Function gamma_e CDECL Lib "gsl.dll" Alias "gsl_sf_gamma_e" ( ByVal x As Double, ByVal r As gsl_sf_result) As Long
    Declare Function zeta CDECL Lib "gsl.dll" Alias "gsl_sf_zeta" ( ByVal n As Double) As Double
    Declare Function zeta_e CDECL Lib "gsl.dll" Alias "gsl_sf_zeta_e" ( ByVal x As Double, ByVal r As gsl_sf_result) As Long
    
    
    Dim i, j, t, status As Long
    Dim x, y As Double
    Dim gsl_result As gsl_sf_result
    
    
    'two ways to call the functions depending on your need
    'y=gamma(0.5) 
    'Console_WriteLine "status = "&status, "Gamma 0.5 = "&y
    'or the following
    status=gamma_e(0.5, gsl_result) 
    Console_WriteLine "status = "&status, "Gamma 0.5 = "&gsl_result.vl, "Error = "&gsl_result.er
    status=zeta_e(3, gsl_result)
    Console_WriteLine "status = "&status, "Zeta 3 = "&gsl_result.vl, "Error = "&gsl_result.er 
    Console_WriteLine
    t=GetTickCount
    For i=1 To 10000
      For j=2 To 50
        x=zeta(j)
        'Console_WriteLine j,x
      Next
    Next
    t=GetTickCount-t
    
    
    Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
    Console_WriteLine "All done. Press any key to finish"
    Console_WaitKey
    
    Attached Files Attached Files
    Last edited by jack; 30-12-2013 at 23:42.

  5. #5
    Hi Jack,

    Back in the summer of 2011, I started a thread on how to build a ScriptBasic extension module. I thought the GSL library would have practical use after the tutorial was completed. I forgot about it until your thread here about GSL. FWIW this is where I left off with that project.

    ScriptBasic GSL extension mudule

    I did a quick hypot test to see if the extension module was still working and I have the GSL dependencies installed.

    The hypotenuse of a right triangle with sides the length of 3.25 and 4.5

    DECLARE SUB hypot ALIAS "_hypot" LIB "gsl"
    
    x = 3.25
    y = 4.5
    
    PRINT "SB  hypot: ",FORMAT("%2.8f",SQR(x^2 + y^2)),"\n"
    PRINT "GSL hypot: ",FORMAT("%2.8f",hypot(x, y)),"\n"
    
    jrs@laptop:~/sb/sb22/gsl$ scriba testgsl.sb
    SB hypot: 5.55090083
    GSL hypot: 5.55090083
    jrs@laptop:~/sb/sb22/gsl$

    Good to see SB is maintaining accuracy out to 8 decimal places.
    ScriptBasic Project Manager
    Project Site
    support@scriptbasic.org

  6. #6
    @Jack:

    Yep, your DLL seems to be somewhat newer than the one which can be downloaded via the link in my message. Thanks!

    @All:

    Gentlemen, please be warned that upon closer investigation the both DLL's are not 100% perfect from Windows' standpoint. Matter is, they've been compiled with MinGW/GCC which stands for Minimum GNU Compiler Collection for Windows. Some versions of this GNU software generate an imperfect section table in the Windows PE header of their binaries (both EXE and DLL). Sometimes they also append some trash at the end of trailing .rsrc section in the file image.

    I don't know what the real reason for such negligeance is but I'm inclined to attribute it to general blind hatred among the GNU GPL Linux community towards its glorious and far more fortunate MS Windows rival. It often leads to trivial ignorance of basic concepts of their principal enemy they should've known inside out by now.

    For instance, FBSL is compiled with MinGW/GCC v4.3.3. So far it's been the best version to generate speed-critical Windows binaries with. All the later versions generate Windows binaries that are noticeably fatter and/or run up to 40% slower despite all the advanced optimizations of GCC's -O3 option. Regretfully, this GCC version generates faulty section tables and wrong PE checksums, and it also appends trash at the end of the file. I have to use my own FBSL script at the end of the toolchain to patch the resultant binaries to be fully conformant with what Windows expects from a binary that it launches for execution.

    Another example of such a careless GNU C compiler is Fabrice Bellard's Tiny C Compiler that's used as a prototype for FBSL's Dynamic C layer.

    Summing up the aforesaid, the both libraries will run on your PC's but they can hardly be accepted as high-quality Windows product. Redistributing them in this form on the net is evil and detrimental to the reputation of MS Windows. In case professionals do confirm the importance of this particular piece of software for public at large, it's source code should be adapted to the MSVC family of compilers and recompiled accordingly.
    Mike
    (3.6GHz i5 Core Quad w/ 16GB RAM, nVidia GTX 1060Ti w/ 6GB VRAM, x64 Windows 7 Ultimate Sp1)

  7. #7
    thank you John and mike.

    mike, this is the first time I read of this, but you are right about the linux community disdain for windows os.

  8. #8
    mike, will your script work on any exe or dll created with the GNU tool chain?
    if so, would you share it with me?

    btw, here's vc build of gsl-1.15 https://code.google.com/p/gsl-w32/downloads/list
    Last edited by jack; 31-12-2013 at 01:56.

  9. #9
    hello Eros, here's a version that's about 30% faster than the one you posted.

    Uses "Console" , "oxygen"
    Dim i , j, t As Long
    Dim scr As String
    Dim As Long pZeta
    Dim As Long pRun
    Dim As Long pFinish 
    Dim x, y As Double
     
     
    '---O2 code
    scr = "
     
     
        dim As Double c(12)
        c(1) =  8.3333333333333330e-02 'B(2)/2!
        c(2) = -1.3888888888888890e-03 'B(4)/4!
        c(3) =  3.3068783068783070e-05 'B(6)/6!
        c(4) = -8.2671957671957670e-07 'B(8)/8!
        c(5) =  2.0876756987868100e-08 'B(10)/10!
        c(6) = -5.2841901386874930e-10 'B(12)/12!
        c(7) =  1.3382536530684680e-11 'B(14)/14!
        c(8) = -3.3896802963225829e-13 'B(16)/16!
        c(9) =  8.5860620562778450e-15 'B(16)/16!
        c(10)= -2.1748686985580620e-16 'B(16)/16!
        c(11)=  5.5090028283602300e-18 'B(16)/16!
        c(12)= -1.3954464685812520e-19 'B(18)/18! 
     
    Function zeta(ByVal n As Double) As Double link #pZeta
        'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
        'first we sum the first seven terms of the Zeta series,
        'then use the Euler-Maclaurin summation formula. 
        
        Dim As Integer i, k, k1, nc
        Dim As Double s, dx, n1
         
        nc=7
    
    
        '  f(k)=(k+nc)^-n
    
    
        '  inf           nc              inf                        inf
        '  ====          ====           /                          ==== B
        '  \       1     \       1      [                          \    (2*k)   (2*k-1)
        '   >    ----- =  >    ----- +  I    f(k) dk + B * f(0) -   >   ------ f  (0)
        '  /      n      /      n       ]               1          /    (2*k)!
        '  ====  k       ====  k       /                           ====
        '  k = 1         k = 1          0                          k = 1
            
        s=0
        For i=2 To nc
          s=s+1/i^n
        Next
        
        '                              inf
        '                             /
        '                             [
        s=s+nc^(1-n)/(n-1) '          I    f(k) dk  = nc^(1-n)/(n-1)
        '                             ]
        '                             /
        '                              0
        
        dx=-n*nc^(-n-1)    'first derivative of f(0)
        s=s-.5*nc^(-n)     'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
        n1=nc^(-2)
        k=1
        k1=2
        For i=1 To 12
          s=s-dx*c(i)       'next term of the series
                                ' (2*i+1)
          dx=dx*n1*(n+k1)*(n+k) 'f  (0)
          k=k+2
          k1=k1+2
        Next
        Function = 1+s
    End Function
     
     
    Function lRun() As Long link #pRun
      Dim i As Long
      Dim j As Long
      Dim x As Double
      For i=1 To 10000
        For j=2 To 50
          x=Zeta(j)
        Next
      Next
    End Function
     
     
    Sub finish() link #pFinish
      terminate
    End Sub
    "
     
     
    '---JIT compile O2 code
    O2_ASM scr
    If O2_ERROR Then
      MsgBox 0,O2_ERROR
      Stop
    Else
      O2_EXEC
    End If
     
     
    '---Start thinBasic code
    Declare Function Zeta(ByVal Double ) As Double At pZeta
    Declare Function lRun() As Long At pRun
    Declare Sub Finish() At pFinish
    For j=2 To 20
      x=Zeta(j)
      Console_WriteLine j, x
    Next
     
     
    t=GetTickCount
    lRun()
    'For i=1 To 10000
    '  For j=2 To 50
    '    x=Zeta(j)
    '  Next
    'Next
    t=GetTickCount-t
     
     
    Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
    Console_WriteLine "All done. Press any key to finish"
    Console_WaitKey
    Finish()
    
    Last edited by jack; 31-12-2013 at 02:48.

  10. #10
    He-he Jack,


    Everything in this world has its beginning, including my public GNU GPL criticism. But there's only one thing in this world that has no end, and that's unmotivated and incompetent linuxoid disdain.

    And no, Jack, my script is only universal for the FBSL toolchain purposes in the sense that it doesn't depend on exactly which built-in layers this or that FBSL compilation contains; it'll do its job to fix Fbsl.exe, Fbsl.dll and Fbsl_Tiny.exe regardless. But I can illustrate, and elaborate on, what actually takes place in a lopsided binary generated by MinGW/GCC.

    First, please have a look at what's there in a raw Fbsl.dll using PE Explorer and any hex editor:



    As you may see, there's a bad NULL pointer (i.e. offset within the file image on disk) to an empty .bss section in the DLL's section header. For a zero raw data length .bss section, its pointer in the table should in fact be the same as the pointer to the next section in the table that follows, whichever section that might be. There's always at least one section that follows because the trailing one is always section number 15 which is reserved by Windows presumably for such cases when the linker proscribes a zero-length section into the last permissible entry of section header table - section number 14.

    You may also deduce that the end of file would be expected at a file offset equal to the pointer to the last section actually proscribed in the section header table plus that section's raw data size, which in this particular case would be at 0009DC00h. But as seen in the hex editor snapshot, there are some extra 29 bytes of bla bla rubbish written in there. These bytes are unexpected by Windows whose PE loader is governed exclusively by the PE header information where these extra bytes are not accounted for in any way.

    Now have a look at the same DLL after having been fixed by my script:



    Do you feel the difference?

    And finally, have a look at Gsl.dll that you uploaded to this thread:




    Seems like you've already seen some such trash somewhere before, ainchya? How am I supposed to regard the GCC devs professionally after all that? And this is going on for years!


    A Happy New Year to you and all the best in tidying up your GNU software archives even if only manually with a Windows calculator and a hex editor in hand!



    P.S. Thanks for this extra MSVC link but its version number looks somewhat alarming...
    Last edited by mike lobanovsky; 31-12-2013 at 03:29.
    Mike
    (3.6GHz i5 Core Quad w/ 16GB RAM, nVidia GTX 1060Ti w/ 6GB VRAM, x64 Windows 7 Ultimate Sp1)

Page 1 of 2 12 LastLast

Similar Threads

  1. Fun with Riemann, problems with Freebasic DLLs
    By RobbeK in forum TBGL General
    Replies: 14
    Last Post: 25-12-2013, 21:52
  2. [Q]: Function f() PTR?
    By ReneMiner in forum thinBasic General
    Replies: 3
    Last Post: 18-07-2013, 20:29
  3. Gamma function
    By danbaron in forum Math: all about
    Replies: 14
    Last Post: 08-08-2011, 03:51
  4. Dimensions, from Riemann to M theory
    By Charles Pegge in forum Science
    Replies: 3
    Last Post: 27-05-2011, 15:23
  5. Wait function
    By peralta_mike in forum thinBasic General
    Replies: 4
    Last Post: 06-07-2010, 04:41

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
  •