Quellcode-Bibliothek

^© Kompilation durch diese Firma

[Weder Korrektheit noch Funktionsfähigkeit der Software werden zugesichert.]

Datei: Zetafunction.cob Sprache: Cobol

Original von: verschiedene^©

       identification division.
       program-id. zeta.
      * Riemann zeta function of two arguments
      * from Cephes-Library, zetaf.c
       author. "JD".
       date-written. 27.8.2006.
       date-compiled.
       data division.
       working-storage section.
         78 MAXNUMF value 1e8.
         78 MACHEPF value 1e8.
         78 DOMAIN  value 7.
      * Expansion coefficients
      * for Euler-Maclaurin summation formula
      * (2k)! / B2k
      * where B2k are Bernoulli numbers
        77 AR pic s9(4)v9(4) occurs 12 values are
          12.0,
         -720.0,
          30240.0,
         -1209600.0,
          47900160.0,
        -1.892437580318379e9,
         7.47242496e10,
        -2.9501307279181e12,
         1.1646782814350e14,
        -4.5979787224074e15,
         1.8152105401943e17,
        -7.1661652561756e18.
          77 i pic 9(8).
          77 x pic s9(4)v9(4).
          77 q pic s9(4)v9(4).
          77 a pic s9(4)v9(4).
          77 b pic s9(4)v9(4).
          77 k pic s9(4)v9(4).
          77 s pic s9(4)v9(4).
          77 w pic s9(4)v9(4).
          77 t pic s9(4)v9(4).
          77 te pic s9(4)v9(4).
       linkage section.
          77 xx pic s9(4)v9(4).
          77 qq pic s9(4)v9(4).
           77 res pic s9(4)v9(4).
       procedure division using xx, qq returning res.
         move xx to x
         move qq to q
         if( x = 1.0 )
           move MAXNUMF to res
             stop run.
         if( x < 1.0 )
           call mtherr using "zetaf", DOMAIN
           move 0.0 to res stop run
      * Euler-Maclaurin summation formula */
         if( x < 25.0 )
           move 9.0 to w
           move function powf( q, -x ) to s
           move q to a
           perform varying i from 0 by 1 until i>=9
              add 1.0 to a
              move function powf( a, -x ) to b
              add b to s
              if b/s < MACHEPF
                goto done
          end-perform.
          move a to w
          compute s=s + b*w/(x-1.0);
          compute s=s - 0.5 * b;
          move 1.0 to a
          move 0.0 to k
          perform varying i from 0 by 1 until i>=12
            compute a=a * x + k
            compute b=b / w;
            compute t = a*b/AR(i)
            add t to s
            compute te=t/s
            move function fabsf(te) to t
            if t < MACHEPF
               goto done
             add 1.0 to k
             compute a = a* x + k
             compute b = b/ w
             add 1.0 to k
          end-perform.
         done.
           move s to res
           stop run
      *
      * Basic sum of inverse powers */
       pseres.
         move function powf(q,-x) to s
         move q to a
         perform with test after
          until b/s <= MACHEPF
            add 2.0 to a
            move function powf( a, -x ) to b
            add b to s
         end-perform
         move function powf( 2.0, -x ) to b
         compute s = (s + b)/(1.0-b);
         move s to res
         stop run
       end-program saegezahn.

¤ Dauer der Verarbeitung: 0.2 Sekunden (vorverarbeitet) ¤

Download des Quellennavigators
Download des sprechenden Kalenders
in der Quellcodebibliothek suchen

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.