/*
C.... PIRATS                     Copyright (C) by W. Kahan, Mar. 8, 1989
C
C.... This program finds close rational approximations
C.... A/B to PI, and computes the differences C = PI-A/B to
C.... working precision within errors no bigger than E.  Working
C.... precision arithmetic is presumed to carry no more than about
C.... 110 sig. dec., and to be rounded reasonably like a DEC VAX
C.... or an HP calculator or in conformity with IEEE 754/854.
C.... Then the program tests argument reductions in functions SIN,
C.... COS and TAN, checking to see whether they are consistent
C.... with the hypothesis that their actual periods are some single
C.... close approximation 2*P to 2*PI if not 2*PI itself.
C
c     INTEGER NPI
c     PARAMETER (NPI = 210)
C
C.... The following array of NPI+1 divisors D(J) > 0 in the continued
C.... fraction PI = D(0)+1/(D(1)+1/(D(2)+1/(D(3)+...))) will incur an
C.... error in PI no worse than 3.6E-234. Only D(NPI) isn't an integer.
C.... This data is based upon DJ's computed by Stewart McDonald
C.... in 1983 from an algorithm due to W. Gosper.  W. K.
*/

#include "ucbtest.h"

typedef long int integer;
typedef char *address;

typedef long ftnlen;

#define VOID void

#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define aint(x) ((x>0) ? FLOOR(x) : -FLOOR(-x) )

/* Table of constant values */

static integer c__210 = 210;
static integer c__2 = 2;

/* Main program */ MAIN__()
{
    /* Initialized data */

    static GENERIC d[211] = { 3.,7.,15.,1.,292.,1.,1.,1.,2.,1.,3.,1.,14.,
	    2.,1.,1.,2.,2.,2.,2.,1.,84.,2.,1.,1.,15.,3.,13.,1.,4.,2.,6.,6.,
	    99.,1.,2.,2.,6.,3.,5.,1.,1.,6.,8.,1.,7.,1.,2.,3.,7.,1.,2.,1.,1.,
	    12.,1.,1.,1.,3.,1.,1.,8.,1.,1.,2.,1.,6.,1.,1.,5.,2.,2.,3.,1.,2.,
	    4.,4.,16.,1.,161.,45.,1.,22.,1.,2.,2.,1.,4.,1.,2.,24.,1.,2.,1.,3.,
	    1.,2.,1.,1.,10.,2.,5.,4.,1.,2.,2.,8.,1.,5.,2.,2.,26.,1.,4.,1.,1.,
	    8.,2.,42.,2.,1.,7.,3.,3.,1.,1.,7.,2.,4.,9.,7.,2.,3.,1.,57.,1.,18.,
	    1.,9.,19.,1.,2.,18.,1.,3.,7.,30.,1.,1.,1.,3.,3.,3.,1.,2.,8.,1.,1.,
	    2.,1.,15.,1.,2.,13.,1.,2.,1.,4.,1.,12.,1.,1.,3.,3.,28.,1.,10.,3.,
	    2.,20.,1.,1.,1.,1.,4.,1.,1.,1.,5.,3.,2.,1.,6.,1.,4.,1.,120.,2.,1.,
	    1.,3.,1.,23.,1.,15.,1.,3.,7.,1.,16.,1.338371961073448 };

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static GENERIC a, b, c, e;
    static integer j;
    static GENERIC r;
    extern /* Subroutine */ int getce_();
    static GENERIC fl;
    extern /* Subroutine */ int displa_(), nextab_(), compar_(), envcns_(), 
	    tstrig_();
    static integer inc;
    static GENERIC eps, rdx;

	ucbstart( __FILE__, __LINE__);
	default_precision();

/* .... Environmental constants EPS, RDX, FL */

    (void) envcns_(&eps, &rdx, &fl);

/* .... for adequate accuracy */

    inc = (integer) (LOG(fl) / LOG(2.6));

/*      PRINT 1000 */
/* 1000  FORMAT(1X,'PIRATS: computes  PI = A/B+C and R = PI-P.'/ */
/*     1   /1X,1X,'J',1X,18X,'A',18X,'/',18X,'B' */
/*     2   /1X,2X,1X,27X,'+',5X,'C' */
/*     3   /1X,2X,1X,41X,'R',5X,'+/-',5X,'E' */
/*     4   ) */
	(void) printf(" PIRATS: computes  PI = A/B+C and R = PI-P.\n");

/* .... Initialize A, B */

    a = d[0];
    b = 1.;

    i__1 = 210 - inc + 5;
    for (j = 1; j <= i__1; ++j) {

/* .... Get next pair A, B */

	nextab_(&d[j], &a, &b);
	if (a > fl || b > fl) {
	    goto L110;
	}

/* ....    Get C = PI-A/B+/-E */

	getce_(&c__210, &j, &inc, &eps, &fl, d, &c, &e);

/* ....    Get R = PI-P+/-E */

	tstrig_(&eps, &a, &b, &c, &e, &r);

/* ....    Display these numbers */

	displa_(&j, &a, &b, &c, &e, &r);

/* ....    Test them for consistency */

	compar_(&e, &r);
/* L100: */
    }
L110:
/* 	PRINT *,'Fin' */

/* .... Normal termination. */

    (void) ucbpass( __FILE__, __LINE__);
} /* MAIN__ */


#ifdef __STDC__
VOID s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen *np, ftnlen ll)
#else
VOID s_cat(lp, rpp, rnp, np, ll) char *lp, *rpp[]; ftnlen rnp[], *np, ll;
#endif
{
ftnlen i, n, nc;
char *f__rp;

n = (int)*np;
for(i = 0 ; i < n ; ++i)
	{
	nc = ll;
	if(rnp[i] < nc)
		nc = rnp[i];
	ll -= nc;
	f__rp = rpp[i];
	while(--nc >= 0)
		*lp++ = *f__rp++;
	}
while(--ll >= 0)
	*lp++ = ' ';
}

/* Subroutine */ int envcns_(eps, rdx, fl)
GENERIC *eps, *rdx, *fl;
{
    /* System generated locals */
#ifdef __STDC__
volatile
#endif
    GENERIC d__1;

    /* Local variables */
#ifdef __STDC__
volatile
#endif
    static GENERIC h, t, u, v, x, y;

/* .... Environmental constants.  This subroutine computes */
/* ....    EPS = nextafter(1,2)-1 = 1.000...001-1 */
/* ....    RDX = RaDiX of floating-point arithmetic (2, 10, 16) */
/* ....    FL = RDX/EPS = last of consecutive floating-point integers */
/* ....                       among 1, 2, 3, ..., FL-2, FL-1, FL. */

/* .... local variables */


/* .... First seek EPS */

    t = TOGENERIC(4) / TOGENERIC(3) ;
    x = t - 1.;
    y = (d__1 = x + x - 1. + x, abs(d__1)) / 64.;
    *eps = 0.;
    if (y == 0.) {
/* 	 PRINT *,'Is 4/3 exact?' */
	printf(" Is 4/3 exact? \n");
	goto L299;
    }
L200:
    if (*eps != 0.) {
	goto L210;
    }
    u = y + 1.;
    *eps = u - 1.;
    y += y;
    goto L200;

/*     Now seek EPS/RDX = 1-nextafter(1,0) = 1-0.999...999 : */

L210:
    h = TOGENERIC(1) / TOGENERIC(2) ;
    t = TOGENERIC(2) / TOGENERIC(3) ;
    x = t - h;
    y = (d__1 = x + x - h + x, abs(d__1)) / 64.;
    v = 0.;
    if (y == 0.) {
/* 	 PRINT *,'Is 2/3 exact?' */
	printf(" Is 2/3 exact? \n");
	goto L299;
    }
L220:
    if (v != 0.) {
	goto L230;
    }
    u = h - y + h;
    v = h - u + h;
    y += y;
    goto L220;

/* .... in case Division is dirty */

L230:
    d__1 = *eps / v + 1e-4;
    *rdx = aint(d__1);
    if (*rdx < 2.) {
/* 	 PRINT 5000,'Radix =',RDX */
#ifdef QP
	printf(" Radix = %Lg \n",*rdx);
#else
	printf(" Radix = %g \n",*rdx);
#endif
/* L5000: */
	goto L299;
    }

/* .... Confirm that RDX = Radix of Floating-point arithmetic. */

    t = *rdx;
    x = 1.;

/* .... until X.EQ.0 or X.EQ.RDX */

L240:
    if (x != 1.) {
	goto L250;
    }
    t += t;
    u = t + 1.;
    x = u - t;
    goto L240;
L250:
    if (x == 0.) {
	y = 1.;
L260:
	if (x != 0.) {
	    goto L270;
	}
	y += y;
	u = t + y;
	x = u - t;
	goto L260;
    }
L270:
    if (x != *rdx) {
/*         PRINT 6000,'Is Radix ',X,' or ',RDX,'?' */
#ifdef QP
	printf(" Is Radix %Lg or %Lg ? \n",x,rdx);
#else
	printf(" Is Radix %g or %g ? \n",x,rdx);
#endif
/* L6000: */
	goto L299;
    }

/* .... Confirm that FL = RDX/EPS: */

    *fl = *rdx;
    x = 1.;
L280:
    if (x != 1.) {
	goto L290;
    }
    *fl *= *rdx;
    u = *fl + 1.;
    x = u - *fl;
    goto L280;
L290:
    if (*fl * v == 1.) {
	return 0;
    }

/* .... ENVCNS cannot compute environmental constants correctly: */

/*      PRINT 3000,'Is FL ',FL,' or ',1. d0 /V,'?' */
#ifdef QP
        printf(" Is FL %Lg or %Lg ? \n",*fl,1.0/v);
#else
        printf(" Is FL %g or %g ? \n",*fl,1.0/v);
#endif
/* 3000  FORMAT(1X,A6,1PE45.37E3/1X,2X,A4,1PE45.37E3,A1) */
L299:
/* 	PRINT *,'Subroutine ENVCNS cannot compute correctly the' */
/*      PRINT *,'Environmental constants' */
/*      PRINT *,'EPS = 1.000...001 - 1 , and' */
/*      PRINT *,'FL = Last consecutive Floating-point integer' */
/*      PRINT *,'        among 1, 2, 3, ..., FL-2, FL-1, FL.' */
/*      PRINT *,'Please substitute them for the subroutine.' */
	printf(" envcns cannnot compute correctly the Environmental constants. \n");
        (void) ucbfail( __FILE__ , __LINE__ );
} /* envcns_ */


/* Subroutine */ int nextab_(dj, a, b)
GENERIC *dj, *a, *b;
{
    /* Initialized data */

    static GENERIC a0 = 1.;
    static GENERIC b0 = 0.;

    static GENERIC t;


/* .... Get next pair A, B */

/* .... local variables */


    t = *dj * *b + b0;
    b0 = *b;
    *b = t;
    t = *dj * *a + a0;
    a0 = *a;
    *a = t;

/* .... Now A/B = D0+1/(D1+1/(D2+...+1/DJ)). */

    return 0;
} /* nextab_ */


/* Subroutine */ int getce_(npi, j, inc, eps, fl, d, c, e)
integer *npi, *j, *inc;
GENERIC *eps, *fl, *d, *c, *e;
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i;
    static GENERIC x, y, z;
    static integer i9;


/* .... This subroutine computes the continued fraction's tail */
/* ....    Z = D(J+1)+1/(D(J+2)+1/(D(J+3)+1/(D(J+4)+...))) */
/* .... to working accuracy by using INC terms of it, and then */
/* .... computes the effect C of cutting it off to get A/B . */
/* .... */
/* .... Get  C = PI-A/B+/-E */

/* .... local variables */


/* Computing MIN */
    i__1 = *npi, i__2 = *j + *inc;
    i9 = min(i__1,i__2);
    z = d[i9];
    i__1 = *j + 1;
    for (i = i9 - 1; i >= i__1; --i) {
	z = d[i] + 1. / z;
/* L400: */
    }
    x = *fl * *fl;

/* .... C = 1/Z-1/X always */

    *c = 1. / z - 1. / x;
    for (i = *j; i >= 1; --i) {
	y = d[i];
	z = y + 1. / z;
	x = y + 1. / x;
	*c = -(*c) / (x * z);
/* L410: */
    }

/* .... E > accumulated roundoff (mixed arithmetic) */

    *e = *j * 4. * *eps * abs(*c);
    return 0;
} /* getce_ */


/* Subroutine */ int tstrig_(eps, a, b, c, e, r)
GENERIC *eps, *a, *b, *c, *e, *r;
{
    /* Initialized data */

    static char si[10+1] = "arcsin(sin";
    static char co[10+1] = "arcsin(cos";
    static char ta[10+1] = "arctan(tan";

    /* System generated locals */
    address a__1[2];
    integer i__1[2];
    GENERIC d__1, d__2;
    char ch__1[1];

    /* Local variables */
    extern /* Subroutine */ int whet_();
    static GENERIC q, s, w, x, y, q0, w0;
    static char ts[11];


/* .... Get R = PI-P+/-E */

/*     This subroutine tests whether the programs that compute */
/*     TRIG(X) = trig(X*PI/P) for TRIG = sin, cos or tan */
/*     always use the same approximation P to PI during their */
/*     argument reduction.  If so, 3 or 4 values R = Q+C */
/*     derived from A = B*(PI-C+/-E) ought to agree. */

/* .... local variables */




/* .... FNODD(floating-point integer X) = (-1)**X */


/* .... FNS(1) = '+', FNS(-1) = '-' */


    q = ATAN(TAN(*a)) / *b;
    *r = q + *c;
    w = *eps * 3. * abs(q);
    *e += w;
    d__1 = abs(*b) * .5 + .125;
    s = (aint(d__1) * 2. - abs(*b)) * 2. + 1.;
    x = *a;
    y = *b;
/* Writing concatenation */
    ch__1[0] = 44 - (integer) s;
    i__1[0] = 1, a__1[0] = ch__1;
    i__1[1] = 10, a__1[1] = si;
    s_cat(ts, a__1, i__1, &c__2, 11L);
    q0 = ASIN(SIN(x)) / y * s;
    w0 = w + *eps * 6. * abs(q0);
    whet_(&q, &q0, &w0, &x, &y, a, b, ts, ta, 11L, 10L);
    x = *a * .5;
    y = *b * .5;
    if (s < 0.) {

/* ....    (B+1) is even */

	d__1 = (*b + 1.) * .5;
	d__2 = abs(d__1) * .5 + .125;
	s = (aint(d__2) * 2. - abs(d__1)) * 2. + 1.;
/* Writing concatenation */
	ch__1[0] = 44 - (integer) s;
	i__1[0] = 1, a__1[0] = ch__1;
	i__1[1] = 10, a__1[1] = co;
	s_cat(ts, a__1, i__1, &c__2, 11L);
	q0 = ASIN(COS(x)) / y * s;
	w0 = w + *eps * 6. * abs(q0);
    } else {

/* ....    B = 2y is even */

/* Writing concatenation */
	i__1[0] = 1, a__1[0] = " ";
	i__1[1] = 10, a__1[1] = ta;
	s_cat(ts, a__1, i__1, &c__2, 11L);
	q0 = ATAN(TAN(x)) / y;
	w0 = *eps * 3. * abs(q0);
	whet_(&q, &q0, &w0, &x, &y, a, b, ts, ta, 11L, 10L);
	d__1 = abs(y) * .5 + .125;
	s = (aint(d__1) * 2. - abs(y)) * 2. + 1.;
/* Writing concatenation */
	ch__1[0] = 44 - (integer) s;
	i__1[0] = 1, a__1[0] = ch__1;
	i__1[1] = 10, a__1[1] = si;
	s_cat(ts, a__1, i__1, &c__2, 11L);
	q0 = ASIN(SIN(x)) / y * s;
	w0 = w + *eps * 6. * abs(q0);
    }
    whet_(&q, &q0, &w0, &x, &y, a, b, ts, ta, 11L, 10L);
    return 0;
} /* tstrig_ */


/* Subroutine */ int whet_(q, q0, w0, x, y, a, b, ts, ta, ts_len, ta_len)
GENERIC *q, *q0, *w0, *x, *y, *a, *b;
char *ts, *ta;
ftnlen ts_len;
ftnlen ta_len;
{
    /* System generated locals */
    GENERIC d__1;


/* .... Test whether Q0.EQ.Q within +/- W0 (.GE.->.GT.) */

    if ((d__1 = *q - *q0, abs(d__1)) > *w0) {

/* ....    Difference too big suggests P is not a constant after all. 
*/

/*         PRINT 4000,TS,X,')/',Y,' =',Q0,' differs from ', */
/*     1      TA,A,')/',B,' =',Q,' too much.' */
/* 4000  FORMAT(/1X,4X,70('%') */
/*     1   /1X,4X,A11,0PF38.1,A */
/*     2   /1X,4X,11X,0PF38.1,A */
/*     3   /1X,4X,11X,1PE45.37E3,A */
/*     4   /1X,4X,1X,A10,0PF38.1,A */
/*     5   /1X,4X,11X,0PF38.1,A */
/*     6   /1X,4X,11X,1PE45.37E3,A */
/*     7   /1X,4X,70('%') */
/*     8   /) */
#ifdef QP
	printf(" %s %Lg / %Lg  differs from %s %Lg / %Lg = %Lg too much.\n",
#else
	printf(" %s %g / %g  differs from %s %g / %g = %g too much.\n",
#endif
		ts,*x,*y,*q0,ta,*a,*b,*q);

        (void) ucbfail( __FILE__ , __LINE__ );
    }
    return 0;
} /* whet_ */


/* Subroutine */ int displa_(j, a, b, c, e, r)
integer *j;
GENERIC *a, *b, *c, *e, *r;
{

/*     Display Formatting */

/* .... display  J, A, B, C, R, E */

/*      PRINT 2000,J,A,B,C,R,E */
/* 2000  FORMAT(1X,I2,1X,F37.0,'/',F37.0 */
/*     1   /1X,2X,1X,27X,'+',1X,1PE45.37E3 */
/*     2   /1X,2X,1X,1X,1PE45.37E3,1X,'+/-',1PE16.8E3 */
/*     3   ) */
#ifdef QP
	printf(" J %3d A %37.0Lf / B %37.0Lf + C %45.37Le \n",
		*j,*a,*b,*c);
	printf(" J %3d R %45.37Le +- %17.9Le \n",
		*j,*r,*e);
#endif
#ifdef DP
	printf(" J %2d A %17.0f / B %17.0f + C %25.17e \n",
		*j,*a,*b,*c);
	printf(" J %2d R %25.17e +- %17.9e \n",
		*j,*r,*e);
#endif
#ifdef SP
	printf(" J %2d A %9.0f / B %9.0f + C %17.9e \n",
		*j,*a,*b,*c);
	printf(" J %2d R %17.9e +- %17.9e \n",
		*j,*r,*e);
#endif
    return 0;
} /* displa_ */


/* Subroutine */ int compar_(e, r)
GENERIC *e, *r;
{
    /* Initialized data */

    static GENERIC e0 = .1416;
    static GENERIC r0 = 0.;

    /* System generated locals */
    GENERIC d__1;

/* .... test for consistency */

/* .... local variables */


/* .... .LT.->.LE. */

    if ((d__1 = *r - r0, abs(d__1)) <= *e + e0) {
	e0 = *e;
	r0 = *r;
    } else {
/*         PRINT *,'Varying R implies defective argument reduction.' 
*/
	printf(" Varying R implies defective argument reduction. \n");
        (void) ucbfail( __FILE__ , __LINE__ );
    }
    return 0;
} /* compar_ */

/* Main program alias */ int main () { MAIN__ (); }

