/* From Rosetta Fortran */
test: procedure options (main);

  declare i fixed binary;

  on underflow ;

  put skip list ('Lanczos', 'Builtin' );
  do i = 1 to 10;
     put skip list (lanczos_gamma(i/3.0q0), gamma(i/3.0q0) );
  end;


lanczos_gamma: procedure (a) returns (float (18)) recursive;
    declare a float (18);
    declare pi float (18) value (3.14159265358979324E0);
    declare cg fixed binary initial ( 7 );

    /* these precomputed values are taken by the sample code in Wikipedia, */
    /* and the sample itself takes them from the GNU Scientific Library */
    declare p(0:8) float (18) static initial
         ( 0.99999999999980993e0, 676.5203681218851e0, -1259.1392167224028e0,
         771.32342877765313e0, -176.61502916214059e0, 12.507343278686905e0,
         -0.13857109526572012e0, 9.9843695780195716e-6, 1.5056327351493116e-7 );

    declare ( t, w, x ) float (18);
    declare i fixed binary;

    x = a;

    if x < 0.5 then
       return ( pi / ( sin(pi*x) * lanczos_gamma(1.0-x) ) );
    else
       do;
          x = x - 1.0;
          t = p(0);
          do i = 1 to cg+2;
             t = t + p(i)/(x+i);
          end;
          w = x + float(cg) + 0.5;
          return ( sqrt(2*pi) * w**(x+0.5) * exp(-w) * t );
       end;
  end lanczos_gamma;

end test;