identification division.
       program-id. piarea.
      * Berechnung von PI nach der Formel
      * F=PI*r**2
      * am Viertelkreis
       author. "JD".
       date-written. 15.8.2007.
       date-compiled.
       data division.
       working-storage section.
       77 n   pic 9(4).
       77 i   pic 9(4).
       77 c   pic 9(4)V9(4).
       77 s   pic 9(4)V9(4).
       77 xi  pic 9(4)V9(4).
       77 yi  pic 9(4)V9(4).
       77 fi  pic 9(4)V9(4).
       78 eps pic 9(4)V9(4) value 0.001.
       77 x   pic 9(4)V9(4).
       77 a   pic 9(4)V9(4).
       linkage section.
       77 pi   pic 9(4)V9(4).
       procedure division returning pi.
      * computes pi = summ of all rectangles fi
      * below quarter circle, xi = breadth, yi=height
         compute n=1/eps
         move 0 to s
         perform varying i from 1 by 1 until i=n
           compute xi=i*eps
           compute x=1-xi*xi
           perform root
           move a to yi
           compute fi=eps*yi
           add fi to s
         end-perform
         add 1 to n
         compute pi=4*s
         stop run.
        root.
         move x to a
         subtract x from 1 giving c
         perform until c<=eps
      *    a*a=x*(1-c), c>=0
           compute a = a*(1+0.5*c)
           compute c = c*c*(0.75+0.25*c)
         end-perform.
      *  x*(1-eps)/x<=a*a<x.
        end-program piareae.