/*
 *  test_cos.c ... check bugs in cosine function
 *  
 *
 *  Created by Tomonori KOUYA 08/02/18.
 *  Copyright 2008 Tomonori KOUYA All rights reserved.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <float.h>
#include <fenv.h>

// get the value of cos by specific rounding mode
// rmode ... FE_TONEREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD
double cos_rmode(double x, int rmode)
{
	int current_rmode;
	double ret;

	// get current rounding mode
	current_rmode = fegetround();

	// change rounding mode 
	fesetround(rmode);

	// cos
	ret = cos(x);

	// restore rounding mode
	fesetround(current_rmode);

	return ret;
}

// maximum value of array arg[] 
double dmax(double arg[], int num)
{
	int i;
	double ret;

	ret = arg[0];
	for(i = 1; i < num; i++)
	{
		if(ret < arg[i])
			ret = arg[i];
	}	

	return ret;
}	

// main function
int main(void)
{
	int times;
	double ret[4], aerr[3], rerr[3], dx = 0.0;

	printf("cos(            x            ) =                           (aerr   , rerr   )\n");
	printf("-----------------------------------------------------------------------------\n");
	for(times = 0; times < 100; times++)
	{
		// in case of dx >> 2.8, big differences are found
		dx = 0.1 * times;

		// calc in 4 rounding modes
		ret[0] = cos_rmode(dx, FE_TONEAREST); // RN mode
		ret[1] = cos_rmode(dx, FE_TOWARDZERO);// RZ mode
		ret[2] = cos_rmode(dx, FE_UPWARD);    // RP mode
		ret[3] = cos_rmode(dx, FE_DOWNWARD);  // RM mode

		// aerr ... absolute differences
		aerr[0] = fabs(ret[0] - ret[1]);
		aerr[1] = fabs(ret[0] - ret[2]);
		aerr[2] = fabs(ret[0] - ret[3]);

		// rerr ... relative differences
		rerr[0] = aerr[0];
		rerr[1] = aerr[1];
		rerr[2] = aerr[2];
		if(fabs(ret[0]) > DBL_EPSILON)
		{
			rerr[0] /= fabs(ret[0]);
			rerr[1] /= fabs(ret[0]);
			rerr[2] /= fabs(ret[0]);
		}

		// print values and differences
		printf("cos(%25.17e) = %25.17e (%7.1e, %7.1e)\n", dx, ret[0], dmax(aerr, 3), dmax(rerr, 3));

		// print only if very very large relative differences are found
		if(dmax(rerr, 3) >= 1.0e-2)
		{
			printf("RN, RZ: %25.17e, %25.17e\n", ret[0], ret[1]);
			printf("RP, RM: %25.17e, %25.17e\n", ret[2], ret[3]);
		}
	}

	return 0;
}