GSL sample: View and Control IEEE754 floating-point arithmetic

Last Update: Oct. 4, 2007

Sample code

ieee754fp.c

/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*        View and control                     */
/*           IEEE754 floating-point arithmetic */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-04                    */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include <fenv.h> // since C99
#include <gsl/gsl_ieee_utils.h>

int main(void)
{
	double a, b = cos(1.0 / 10.0), c = sin(1.0 / 11.0);

	printf("b -> Decimal = %25.17e\n", b);
	printf("     Binary  =  ");
	gsl_ieee_printf_double(&b);
	printf("\n");
	printf("c -> Decimal = %25.17e\n", c);
	printf("     Binary  =  ");
	gsl_ieee_printf_double(&c);
	printf("\n\n");

	/* setup ronding mode */
	/* set IEEE754 double precition floating-point number */
	printf("-- Round to -Infinity\n");
	fesetround(FE_DOWNWARD); // Round to -Infinity
	a = b + c;
	printf("b + c -> Decimal = %25.17e\n", a);
	printf("	 Binary  =  ");
	gsl_ieee_printf_double(&a);
	printf("\n\n");

	printf("-- Round to +Infinity\n");
	fesetround(FE_UPWARD); // Round to +Infinity
	a = b + c;
	printf("b + c -> Decimal = %25.17e\n", a);
	printf("	 Binary  =  ");
	gsl_ieee_printf_double(&a);
	printf("\n\n");

	printf("-- Round to Zero\n");
	fesetround(FE_TOWARDZERO); // Round to Zero
	a = b + c;
	printf("b + c -> Decimal = %25.17e\n", a);
	printf("	 Binary  =  ");
	gsl_ieee_printf_double(&a);
	printf("\n\n");

	printf("-- Round to Nearest (default)\n");
	fesetround(FE_TONEAREST); // Round to Nearest (default)
	a = b + c;
	printf("b + c -> Decimal = %25.17e\n", a);
	printf("	 Binary  =  ");
	gsl_ieee_printf_double(&a);
	printf("\n\n");

	return 0;
}

Example execution

$ ./ieee754fp
b -> Decimal =   9.95004165278025821e-01
     Binary  =   1.1111110101110001001011111001101010000001011111000001*2^-1
c -> Decimal =   9.07839235088703650e-02
     Binary  =   1.0111001111011001110101111110011110010010000101000011*2^-4

-- Round to -Infinity
b + c -> Decimal =   1.08578808878689603e+00
         Binary  =   1.0001010111110110001101010100101110111001110111110100*2^0

-- Round to +Infinity
b + c -> Decimal =   1.08578808878689625e+00
         Binary  =   1.0001010111110110001101010100101110111001110111110101*2^0

-- Round to Zero
b + c -> Decimal =   1.08578808878689603e+00
         Binary  =   1.0001010111110110001101010100101110111001110111110100*2^0

-- Round to Nearest (default)
b + c -> Decimal =   1.08578808878689625e+00
         Binary  =   1.0001010111110110001101010100101110111001110111110101*2^0

<- Go back