Page 2 of 2 FirstFirst 12
Results 11 to 14 of 14

Thread: Riemann zeta function in O2

  1. #11
    thank you mike for the info.

    Happy new year

  2. #12
    thinBasic author ErosOlmi's Avatar
    Join Date
    Sep 2004
    Location
    Milan - Italy
    Age
    57
    Posts
    8,777
    Rep Power
    10
    Jack,

    I've developed Zeta function inside thinBasic Math module using you latest O2 code.
    It seems it is a little bit faster even if interpreted: 430 ms here on my PC and results seem OK.

    To test the script, leave thinBasic_Math.dll in the same directory of the script: thinBasic will search for needed modules in script directory first.

    If OK for you, I will release it in next thinBasic update.

    Ciao
    Eros
    Attached Files Attached Files
    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. #13
    Eros
    yes of course it's OK to use my code, but perhaps some boundary checks should be included, maybe you did that already, for double x>52 the result is simply 1 also my function does not handle negative arguments, perhaps it should.
    on second thought, if it's OK with you I will write a Zeta function using Chebyshev approximation, it should be about 8 times faster and I will make it more bullet proof.
    Last edited by jack; 31-12-2013 at 13:24.

  4. #14
    hello Eros.

    here's the zeta function using Chebyshev approximations, it's 3 times faster for positive arguments than the one you used in the Math module, also works for negative arguments.
    code is in PowerBasic Console compiler.

    #COMPILE EXE
    #DIM ALL
    
    
    FUNCTION gamma(BYVAL y AS DOUBLE) AS DOUBLE
        DIM Pi AS STATIC DOUBLE
        DIM sq2pi AS STATIC DOUBLE
        DIM g AS STATIC DOUBLE
        DIM t AS DOUBLE, w AS DOUBLE, gam AS DOUBLE, x AS DOUBLE
        DIM i AS LONG, cg AS LONG
        DIM c(15) AS STATIC DOUBLE
    
    
        Pi    = 3.1415926535897932385
        sq2pi = 2.50662827463100050241577    'sqrt(2Pi)
        g     = 607/128 ' best resu'ts when 4<=g<=5
        x = y-1
        cg = 14
        'Lanczos approximation for the complex plane
        'calculated using vpa digits(256)
        'the best set of coeffs was selected from
        'a solution space of g=0 to 32 with 1 to 32 terms
        'these coeffs really give superb performance
        'of 15 sig. digits for 0<=real(z)<=171
        'coeffs should sum to about g*g/2+23/24
    
    
        'http://www.numericana.com/answer/info/godfrey.htm
    
    
          c( 1) =        0.99999999999999709182
          c( 2) =       57.156235665862923517
          c( 3) =      -59.597960355475491248
          c( 4) =       14.136097974741747174
          c( 5) =       -0.49191381609762019978
          c( 6) =        0.000033994649984811888699
          c( 7) =        0.000046523628927048575665
          c( 8) =       -0.000098374475304879564677
          c( 9) =        0.00015808870322491248884
          c(10) =       -0.00021026444172410488319
          c(11) =        0.00021743961811521264320
          c(12) =       -0.00016431810653676389022
          c(13) =        0.000084418223983852743293
          c(14) =       -0.000026190838401581408670
          c(15) =        0.0000036899182659531622704
    
    
        IF ( x < 0 ) THEN
            x=-x
            IF FRAC(x)=0 THEN
                gam=10^4932
            ELSE
                t = c(1)
                FOR i=1 TO cg
                    t = t + c(i+1)/(x+i)
                NEXT
                w = x + g + 0.5
                gam=w^(x+0.5)*EXP(-w)*sq2pi*t
                gam = Pi*x/(gam*SIN(Pi*x))
            END IF
        ELSE
            t = c(1)
            FOR i=1 TO cg
                t = t + c(i+1)/(x+i)
            NEXT
            w = x + g + 0.5
            gam=w^(x+0.5)*EXP(-w)*sq2pi*t
        END IF
        FUNCTION = gam
    END FUNCTION
    
    FUNCTION zeta(BYVAL n AS DOUBLE) AS DOUBLE
        'Riemann zeta function evaluated using tables from
        'Chebyshev Approximations for the Riemann Zeta Function
        'By W. J. Cody,* K. E. Hillstrom,* and Henry C. Thacher, Jr.**
    
    
        'except for the second set which was computed using Maple.
    
    
        DIM p(8) AS STATIC DOUBLE
        p(0) = 12871681214.82446392809
        p(1) = 13753969320.37025111825
        p(2) = 5106655918.364406103683
        p(3) = 856147100.2433314862469
        p(4) = 74836181.24380232984824
        p(5) = 4860106.585461882511535
        p(6) = 273957.4990221406087728
        p(7) = 4631.710843183427123061
        p(8) = 57.87581004096660659109
    
    
        DIM q(8) AS STATIC DOUBLE
        q(0) = 25743362429.64846244667
        q(1) = 5938165648.679590160003
        q(2) = 900633037.3261233439089
        q(3) = 80425366.34283289888587
        q(4) = 5609711.759541920062814
        q(5) = 224743.1202899137523543
        q(6) = 7574.578909341537560115
        q(7) =-23.73835781373772623086
        q(8) = 1.000000000000000000000
    
    
        DIM p1(8) AS STATIC DOUBLE
        p1(0) = 0.1811021757722157193018e-1
        p1(1) =-0.1781077292672229847925e-2
        p1(2) = 0.3271270368846556844317e-2
        p1(3) =-0.1183785392295491855720e-2
        p1(4) = 0.2811461476445431284369e-3
        p1(5) =-0.4540699985455522705413e-4
        p1(6) = 0.4921410048947204581619e-5
        p1(7) =-3.3060968398926623965732e-7
        p1(8) = 1.0649769599519517944407e-8
    
    
        DIM q1(8) AS STATIC DOUBLE
        q1(0) =-.95084345744468539815114
        q1(1) = .23724094500966523346054
        q1(2) = .12719149298560634692396
        q1(3) =-0.3742718269184690837860e-1
        q1(4) =-0.2745688015667984911198e-1
        q1(5) =-0.7479966479870980001960e-2
        q1(6) =-0.1221861376076105986461e-2
        q1(7) =-0.1189441500804843551022e-3
        q1(8) =-0.6171014074208758609628e-5
    
    
        DIM p2(7) AS STATIC DOUBLE
        p2(0) = 1.66156480515774675916e-11
        p2(1) =-4.68068827660654526862e-09
        p2(2) = 5.83519727319147047318e-07
        p2(3) =-4.17644012643145602124e-05
        p2(4) = 1.85468422843597959483e-03
        p2(5) =-5.11288800220490240591e-02
        p2(6) = 8.10450231751100353193e-01
        p2(7) =-5.69951948768478922618e+00
    
    
        DIM q2(9) AS STATIC DOUBLE
        q2(0) =-6.99562633519191654964e-10
        q2(1) =-1.77757961895149256941e-08
        q2(2) =-9.82231825734078036442e-07
        q2(3) =-2.84927282759096487594e-05
        q2(4) =-5.81727909388048093531e-04
        q2(5) =-1.15848749169766585807e-02
        q2(6) =-1.28149124051978195742e-01
        q2(7) =-1.11913057349097709324e+00
        q2(8) =-7.67928761604628812537e-01
        q2(9) = 1.00000000000000000000e+00
    
    
        DIM p3(8) AS STATIC DOUBLE
        p3(0) = 1.0314487718885971168e-15
        p3(1) =-5.1258461396468824062e-13
        p3(2) = 1.1294879419487354786e-10
        p3(3) =-1.4423466537313095228e-08
        p3(4) = 1.1682467698445809766e-06
        p3(5) =-6.1497516799031480614e-05
        p3(6) = 2.0559467798883032750e-03
        p3(7) =-3.9933942939466886853e-02
        p3(8) = 3.4523497673617845708e-01
    
    
        DIM q3(8) AS STATIC DOUBLE
        q3(0) = 5.9395941728841905020e-11
        q3(1) =-6.0475535907999180572e-09
        q3(2) = 3.6468020866838856275e-07
        q3(3) =-1.2945690556801181241e-05
        q3(4) = 3.2018949847022925001e-04
        q3(5) =-5.0780155709999407748e-03
        q3(6) = 5.4962890788158726560e-02
        q3(7) =-3.2451761115597241852e-01
        q3(8) = 1.0000000000000000000e+00
        DIM pi AS STATIC DOUBLE
        DIM pn AS DOUBLE, qn AS DOUBLE, n2 AS DOUBLE, y AS DOUBLE
        DIM j AS LONG
        pi=3.141592653589793
        IF n<0 THEN
            n2=ABS(n)
            IF FRAC(n2)=0 THEN
                j=n2
                IF j MOD 2=0 THEN
                    FUNCTION=0
                    EXIT FUNCTION
                END IF
            END IF
            n2=1+n2
            FUNCTION=2^n*pi^(-n2)*COS(0.5*pi*n2)*gamma(n2)*zeta(n2)
        ELSEIF n>55 THEN
            FUNCTION=1
        ELSE
        SELECT CASE n
            CASE n=0.5 TO 5
                pn=p(0)
                qn=q(0)
                n2=n
    
    
                FOR j=1 TO 8
                    pn+=n2*p(j)
                    qn+=n2*q(j)
                    n2*=n
                NEXT
                qn*=(n-1)
                FUNCTION=pn/qn
            CASE n=5 TO 11
                y=n/3-8/3
                n2=y*y
                pn=p1(0)+p1(1)*n
                qn=q1(0)+q1(1)*n
                FOR j=2 TO 8
                    pn+=n2*p1(j)
                    qn+=n2*q1(j)
                    n2*=y
                NEXT
                FUNCTION=pn/qn+1
            CASE n=11 TO 25
                y=1/n
                n2=y
                pn=p2(0)+p2(1)*n2
                qn=q2(0)+q2(1)*n2
                FOR j=2 TO 7
                    n2*=y
                    pn+=n2*p2(j)
                    qn+=n2*q2(j)
                NEXT
                n2*=y
                qn+=n2*q2(8)
                n2*=y
                qn+=n2*q2(9)
                FUNCTION=1+2^(y*pn/qn-n)
            CASE n=25 TO 55
                y=1/n
                n2=y
                pn=p3(0)+p3(1)*n2
                qn=q3(0)+q3(1)*n2
                FOR j=2 TO 8
                    n2*=y
                    pn+=n2*p3(j)
                    qn+=n2*q3(j)
                NEXT
                FUNCTION=1+2^(y*pn/qn-n)
        END SELECT
        END IF
    END FUNCTION
    
    
    FUNCTION PBMAIN () AS LONG
        DIM s AS DOUBLE, t AS DOUBLE
        DIM i AS LONG, j AS LONG
    
    
        t=TIMER
        FOR i=1 TO 10000
            FOR j=2 TO 50
                s=zeta(j)
    '            ?j,s
            NEXT
        NEXT
        t=TIMER-t
    
    
        PRINT t*1000;" nSec"
    
    
        INPUT t
    END FUNCTION
    
    fixed copy/paste mistake.
    Last edited by jack; 01-01-2014 at 15:21.

Page 2 of 2 FirstFirst 12

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
  •