/*
 * Copyright 1/5/95 Zhishun Alex Liu.  All Rights Reserved.
 *
 ***************************
 * WORK IN PROGRESS
 ***************************
 *
 * Written by Zhishun Alex Liu under direction of Professor William Kahan.
 * The author's current electronic mail address as of December 1994 is
 * Alex.Liu@Eng.Sun.COM.
 *
 * Use of this software is granted with the understanding that all recipients
 * should regard themselves as participants in an ongoing research project and
 * hence should feel obligated to report their experiences (good or bad) with
 * this elementary functions test suite to the author.
 */
#ifndef lint
static char sccsid[] = "@(#)findpi.c	1.1 (BeEF) 1/5/95";
#endif
#include <stdio.h>
#include <assert.h>
#include "sys.h"
#ifdef __STDC__
extern void exit(int);
#endif

extern int sig_bits;

static LONG piov2_p[] = {	/* no leading sign + power-of-2 info */
	0x1921fb, 0x54442, 0xd1846, 0x9898c, 0xc5170,
	0x1b839, 0xa2520, 0x49c11, 0x14cf9, 0x8e804,
	0x177d4, 0xc7627, 0x3644a, 0x29410, 0xf31c6,
	0x809bb, 0xdf2a3, 0x3679a, 0x74863, 0x66056,
	0x14dbe, 0x4be28, 0x6e9fc, 0x26ada, 0xdaa38,
	0x48bc9, 0x0b6ae, 0xcc4bc, 0xfd8de, 0x89885,
	0xd34c6, 0xfdad6,
};
#define	PCS_PI2	(sizeof(piov2_p) / sizeof(piov2_p[0]))
#define	HML	3
#define	MAX_PI2	4
#define	PI2x(i)	(npiov2 + HML * (i - 1))
static GENERIC npiov2[PCS_PI2 + HML * MAX_PI2 + 1];	/* i * %pi / 2, i = 1..MAX_PI2 */

static int
#ifdef __STDC__
mkpi2(GENERIC *_hi, GENERIC *_lo, char *sp)
#else
mkpi2(_hi, _lo, sp)
GENERIC *_hi, *_lo;			/* input arguments */
char *sp;
#endif
{
	register int i, digit, bit, nbits = 0;
	GENERIC hi = *_hi, lo = *_lo;

	hi -= ONE;		/* leading bit */
	assert(hi >= ZERO && hi < ONE);
	*sp++ = '0'; *sp++ = 'x'; *sp++ = '1';
	for (i = 1, digit = 0; i <= ((2 * sig_bits + 19) / 20) * 20; i++) {
		hi += hi; lo += lo;
		bit = hi ? (int) hi : (int) lo;
		digit <<= 1; digit += bit;
		if (bit) {
			if (hi)
				hi -= (GENERIC) bit;
			else
				lo -= (GENERIC) bit;
		}
		if (!nbits && hi == ZERO && lo == ZERO)
			nbits = i + 1;
		if (i % 4 == 0) {
			*sp++ = "0123456789abcdef"[digit];
			digit = 0;
			if (i % 20 == 0) {
				if (nbits)
					break;
				else {
					*sp++ = ',';
					*sp++ = ' '; *sp++ = '0'; *sp++ = 'x';
				}
			}
		}
	}
	*sp = '\0';
	return nbits ? nbits : 2 * sig_bits;
}

static void
#ifdef __STDC__
findpi(void)
#else
findpi()
#endif
{
	register int chunks, pieces, i, j;
	GENERIC p2m20, p2m20i, x, y_s2pi, y_s1pi, y_cpi2, pi_del, gtmp;
	char s[PCS_PI2 * (5 + 4) + 1];

	chunks = (sig_bits * HML - 1 + (20 - 1)) / 20;	/* need HML good pieces */
	p2m20 = ld_2_n(-20);

	for (j = 1; j <= MAX_PI2; j++) {
		p2m20i = ONE;
		for (i = 0; i < chunks; i++) {
			PI2x(j)[i] = (GENERIC) ((LONG) j * piov2_p[i]);
			PI2x(j)[i] *= p2m20i;
			PI2x(j)[i] *= p2m20;
			p2m20i *= p2m20;
		}
		pieces = distd(chunks, PI2x(j));
		assert(pieces == HML || pieces == HML + 1);
	}
#if TEST_SIN && TEST_COS
/*
 * We now proceed to discover the PI our user uses.
 *
 * Let	%pi	:= true pi,
 *	Pi	:= %pi rounded to WP and
 *	PI	:= machine pi.
 *
 * Let Sin(x)	:= sin((%pi / PI) * x).
 *
 * Sin(Pi)	= sin(%pi - Pi * (%pi / PI))
 *		= sin((PI - Pi) * (%pi / PI))
 * (approx.)	= (PI - Pi) * (%pi / PI)
 * (approx.)	= PI - Pi				(+)
 *
 * Similarly,
 *
 * Sin(2 * Pi)	= -sin(2 * %pi - 2 * Pi * (%pi / PI))
 *		= -sin(2 * (PI - Pi) * (%pi / PI))
 * (approx.)	= -2 * (PI - Pi) * (%pi / PI)
 * (approx.)	= -2 * (PI - Pi)			(+)
 *
 * and
 *
 * Cos(Pi / 2)	= sin(%pi / 2 - (Pi / 2) * (%pi / PI))
 *		= sin(((PI - Pi) / 2) * (%pi / PI))
 * (approx.)	= ((PI - Pi) / 2) * (%pi / PI)
 * (approx.)	= (PI - Pi) / 2				(+)
 *		= (PI - %pi) / 2 + (%pi - Pi) / 2	(++)
 *
 * Notice that since we surely know Pi, the value of PI can be easily
 * determined by (+).  Since we know (%pi - Pi) to twice the WP accuracy,
 * the difference between PI and true %pi can also be determined by (++).
 *
 * The above 3 instances would give very accurate indication of the value
 * of machine PI our user's implementation of sin() and cos() uses - if
 * there is one.
 */
	x = PI2x(4)[0];
	FUNC_SIN(y_s2pi, x)
	y_s2pi /= -FOUR;
	x = PI2x(2)[0];
	FUNC_SIN(y_s1pi, x)
	y_s1pi /= TWO;
	x = PI2x(1)[0];
	FUNC_COS(y_cpi2, x)
	if (y_s2pi == y_s1pi && y_s1pi == y_cpi2) {
		if (y_cpi2 == PI2x(1)[1]) {
#define	Z "/*\n * The PI sin() & cos() use appears to agree with true PI to\n"
			(void) printf(Z);
#undef	Z
#define	Z " * at least twice the working precision, or %d bits.\n */\n"
			(void) printf(Z, 2 * sig_bits);
#undef	Z
		}
		pi_del = y_cpi2 - PI2x(1)[1];
		pi_del -= PI2x(1)[2];		/* machine PI - true PI rounded */
		if (y_cpi2 < ZERO) {
			GENERIC scale, ulp = getulp(&x, &scale);
			x -= ulp;
			y_cpi2 += ulp;
			assert(y_cpi2 >= ZERO);
		}
		for (gtmp = GABS(pi_del), i = 0; gtmp < ONE; i++)
			gtmp += gtmp;
		i--;
		j = mkpi2(&x, &y_cpi2, s);
		(void) printf("#define\tMPI2_P\t%s\t/* %d(->%d)-bit */\n",
			s, j, i);
	}
	else
		{
		(void) printf("/* Cannot determine the PI sin() & cos() use. */\n");
		ucbfail( __FILE__ , __LINE__ );
	}
#else
	(void) printf(" TEST_SIN %d or TEST_COS %d undefined \n",TEST_SIN,TEST_COS);
	ucbfail( __FILE__ , __LINE__ );
#endif
}

#ifdef SUN_LEGACY_KLUDGE
#undef sparc
#undef mc68000
#define	i386	0	/* kludge to bypass FLOATFUNCTIONTYPE double/int */
#endif
#include <math.h>
void
#ifdef __STDC__
ucbfindpitest_ (int *pargc, char *argv[])
#else
ucbfindpitest_ (pargc, argv)
int *pargc; char *argv[];
#endif
{
#ifdef SUN_LEGACY_KLUDGE	/* i.e. SunOS4.0 and later */
	register int bits;
	int argc;

	argc = *pargc;
	if (argc > 1) {
		bits = atoi(argv[1]);
		if (bits == 53)
			fp_pi = fp_pi_53;
		else if (bits == 66)
			fp_pi = fp_pi_66;		/* SunOS4.0 default */
		else
			fp_pi = fp_pi_infinite;
	}
#endif
	bgn_all();
	findpi();
	ucbpass( __FILE__ , __LINE__ );
	(void) exit(0);
	/* NOTREACHED */
}
