Re: [bug-gawk] pi in gawk

[Nelson H. F. Beebe]

To answer Francky's original question: no, awk has no built-in named
mathematical constants.

Arnold suggests getting a value of $\pi$ like this:

>> gawk 'BEGIN { printf ("%.17g\n", atan2(0, -1)) }'

While that is correctly mathematically, it relies on unknown accuracy
of underlying math library functions, and thus, is a bad idea in
portable code.

The best way to represent constants, like $\pi$, $\sqrt{2}$, and so
on, is as a sum of exact and approximate parts, so that compile-time,
or run-time, evaluation of the expression gets the result correct to
the last bit.

Here are two snippets of C code, one for decimal arithmetic, and one
for binary (both assuming IEEE 754 64-bit formats):

    static const fp_t PI_HI = FP(3141592653589793.0) / FP(1.e+15);
    static const fp_t PI_LO = FP(2.384626433832795e-16);

    static const fp_t PI_HI = FP(7074237752028440.0) / FP(2251799813685248.0);
    static const fp_t PI_LO = FP(1.2246467991473532e-16);

The FP() macro wrapper supplies a suitable type suffix.

For awk, which has a default of double, and that is IEEE 754 64-bit
with a 53-bit significand on almost all systems today, one could write

        PI_HI = 7074237752028440.0 / 2251799813685248.0
        PI_LO = 1.2246467991473532e-16
        PI = PI_HI + PI_LO

The high part is defined as a rational number that is exactly
representable in the base and precision of the host arithmetic.  Both
numerator and denominator are also exactly representable, and not
subject to base-conversion errors.

The low part is an accurate correction that can be specified to as
many digits as desired, although there are historical machines that
choked, or worse, got the wrong answer, if you specified `too many'
digits in a numeric constant.

How did I get the numeric values in the first place?  I used a
symbolic algebra system:

% maple
> evalf(Pi, 35);
                      3.1415926535897932384626433832795029
% math
In[1]:= N[Pi, 35]
Out[1]= 3.1415926535897932384626433832795029

% maxima
(%i1) fpprec : 35$
(%i2) bfloat(atan2(0,-1));
(%o2)               3.1415926535897932384626433832795029b0

% gp
? \p 35
   realprecision = 38 significant digits (35 digits displayed)
? Pi
%1 = 3.1415926535897932384626433832795029

The first two are commercial, and widely used and documented in
hundreds of books.  The second two are opensource or freeware; maxima
is available as an optional package in most GNU/Linux distributions.

It doesn't take much more work in those systems to find rational
approximations to constants, along with a correcting term, but you do
have to get somewhat familiar with their idiosyncratic syntax.

Here is how to do it in Maple:

> Digits := 35:

> printf("%a.0 / %a.0\n", trunc(Pi * 2**51), 2**51):
7074237752028440.0 / 2251799813685248.0

> printf("%.35g\n", Pi - trunc(Pi * 2**51) / 2**51):
1.224646799147353177e-16

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: address@hidden  -
- 155 S 1400 E RM 233                       address@hidden  address@hidden -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------