Quelle  Zetafunction.cob   Sprache: Cobol

       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.