Previous Next Contents Index Doc Set Home


Examples

A

This appendix provides examples of how to accomplish some popular tasks. The examples are written either in FORTRAN or ANSI C, and many depend on the current version of libm and libsunmath. These examples were tested with the current C and FORTRAN compilers in the Solaris operating environment.

This appendix is organized into the following sections:

IEEE Arithmetic

page 91

The Math Libraries

page 94

Exceptions and Exception Handling

page 106

Miscellaneous

page 120

IEEE Arithmetic

This section shows how to examine the stored IEEE hexadecimal representation of floating-point numbers.

Examining Stored IEEE hexadecimal Representations

These examples show a way of examining the hexadecimal representation of floating-point numbers. Note that you can also use the debuggers to look at the hexadecimal representation of stored data, as well as write programs based on other ideas.

The following C program prints a double-precision approximation to pi and single-precision infinity: 

#include <math.h>
#include <sunmath.h>
 
main() {
	union {
		float		flt;
		unsigned	un;
	} r;
	union {
		double		dbl;
		unsigned	un[2];
	} d;

	/* double precision */
	d.dbl = M_PI;
	(void) printf("DP Approx pi = %08x %08x = %18.17e \n",
		d.un[0], d.un[1], d.dbl);

	/* single precision */
	r.flt = infinityf();
	(void) printf("Single Precision %8.7e : %08x \n", 
		r.flt, r.un);

	exit(0);
}

The output of this program looks like (SPARC): 

DP Approx pi = 400921fb 54442d18 = 3.14159265358979312e+00 
Single Precision Infinity: 7f800000

And from FORTRAN: 

	program print_ieee_values
c
c the purpose of the implicit statements is to insure
c that the f77_floatingpoint pseudo-intrinsic functions
c are declared with the correct type
c
	implicit real*16 (q)
	implicit double precision (d)
	implicit real (r)
	real*16 z
	double precision x
	real r
c
	z = q_min_normal()
	write(*,7) z, z
 7	format('min normal, quad: ',1pe47.37e4,/,' in hex ',z32.32)
c
	x = d_min_normal()
	write(*,14) x, x
 14	format('min normal, double: ',1pe23.16,' in hex ',z16.16)
c
	r = r_min_normal()
	write(*,27) r, r
 27	format('min normal, single: ',1pe14.7,' in hex ',z8.8)
c
	end

Intel does not support quadruple precision real*16. To run the FORTRAN example on Intel, delete the real*16 declaration and the q_min_normal() function that calculates and prints the quadruple precision value. 

	implicit real*16 (q)
	...
	z = q_min_normal()
	write(*,7) z, z
 7	format('min normal, quad: ',1pe47.37e4,/,' in hex ',z32.32)

The corresponding output looks like (SPARC): 

min normal, quad:   3.3621031431120935062626778173217526026E-4932
   in hex 00010000000000000000000000000000
min normal, double:2.2250738585072014-308 in hex 0010000000000000
min normal, single:1.1754944E-38 in hex 00800000

The Math Libraries

In this section are a few examples that use functions from the math library.

Random Number Generator

The following example calls a random number generator to generate an array of numbers, and uses a timing function to measure the time it takes to compute the EXP of the given numbers: 

#ifdef DP
#define GENERIC double precision
#else
#define GENERIC real
#endif

#define SIZE 400000

	program example
c
	implicit GENERIC (a-h,o-z)
	GENERIC x(SIZE), y, lb, ub
	real time(2), u1, u2
c
c  compute EXP on random numbers in [-ln2/2,ln2/2]
	lb = -0.3465735903
	ub = 0.3465735903
c
c generate array of random numbers
#ifdef DP
	call d_init_addrans()
	call d_addrans(x,SIZE,lb,ub)
#else
	call r_init_addrans()
	call r_addrans(x,SIZE,lb,ub)
#endif
c
c start the clock
	call dtime(time)
	u1 = time(1)
c
c compute exponentials
	do 16 i=1,SIZE
		y = exp(x(i))
 16	continue
c
c get the elapsed time
	call dtime(time)
	u2 = time(1)
	print *,'time used by EXP is ',u2-u1
	print *,'last values for x and exp(x) are ',x(SIZE),y
c
	call flush(6)
	end

Note that compilation must be done using either the -DDP or -DSP flag. The suffix is F instead of f so that the preprocessor is invoked automatically.

This example shows how to use d_addrans to generate blocks of random data uniformly distributed over a user-specified range: 

/*
 * test SIZE*LOOPS random arguments to sin in the range
 * [0, threshold] where
 * threshold = 3E30000000000000 (3.72529029846191406e-09)
 */

#include <math.h>
#include <sunmath.h>
#define SIZE 10000
#define LOOPS 100

main()
{
	double		x[SIZE], y[SIZE];
	int		i, j, n;
	double		lb, ub;
	union {
		unsigned	u[2];
		double		d;
	}  upperbound;

	upperbound.u[0] = 0x3e300000;
	upperbound.u[1] = 0x00000000;

	/* initialize the random number generator */
	d_init_addrans_();

	/* test (SIZE * LOOPS) arguments to sin */
	for (j = 0; j < LOOPS; j++) {

		/*
		* generate a vector, x, of length SIZE, 
		* of random numbers to use as
		* input to the trig functions.
		*/
		n = SIZE;
		ub = upperbound.d;
		lb = 0.0;
		d_addrans_(x, &n, &lb, &ub);

		for (i = 0; i < n; i++)
			y[i] = sin(x[i]);

		/* is sin(x) == x?  It ought to, for tiny x.  */
		for (i = 0; i < n; i++)
			if (x[i] != y[i])
				printf(
					" OOPS: %d sin(%18.17e)=%18.17e \n",
					i, x[i], y[i]);
	}

	ieee_retrospective_();
	exit(0);
}

This FORTRAN example uses some functions recommended by the IEEE standard: 

/*
	Demonstrate how to call 5 of the more interesting IEEE 
	recommended functions from fortran. These are implemented 
	with "bit-twiddling", and so are as efficient as you could 
	hope. The IEEE standard for floating-point arithmetic 
	doesn't require these, but recommends that they be 
	included in any IEEE programming environment.

	For example, to accomplish 
		y = x * 2**n, 
	since the hardware stores numbers in base 2,
	shift the exponent by n places.  
                
	Refer to 

	ieee_functions(3m)
	libm_double(3f)
	libm_single(3f)

	The 5 functions demonstrated here are:

	ilogb(x): returns the base 2 unbiased exponent of x in 
			integer format
	signbit(x): returns the sign bit, 0 or 1
	copysign(x,y): returns x with y's sign bit
	nextafter(x,y): next representable number after x, in 
			the direction y
	scalbn(x,n): x * 2**n

	function        double precision      single precision
	--------------------------------------------------------
	ilogb(x)        i = id_ilogb(x)       i = ir_ilogb(r)
	signbit(x)      i = id_signbit(x)     i = ir_signbit(r)
	copysign(x,y)   x = d_copysign(x,y)   r = r_copysign(r,s)
	nextafter(x,y)  z = d_nextafter(x,y)  r = r_nextafter(r,s)
	scalbn(x,n)     x = d_scalbn(x,n)     r = r_scalbn(r,n)
*/

#include <values.h>

	program ieee_functions_demo
	implicit double precision (d)
	implicit real (r) 
	double precision x, y, z, direction
	real r, s, t, r_direction
	integer i, scale

	print *
	print *, 'DOUBLE PRECISION EXAMPLES:'
	print *

	x = 32.0d0
	i = id_ilogb(x)
	write(*,1) x, i
 1	format(' The base 2 exponent of ', F4.1, ' is ', I2)

	x = -5.5d0
	y = 12.4d0
	z = d_copysign(x,y)
	write(*,2) x, y, z
 2    format(F5.1, ' was given the sign of ', F4.1,
     *   ' and is now ', F4.1)

	x = -5.5d0
	i = id_signbit(x)
	print *, 'The sign bit of ', x, ' is ', i

	x = MINDOUBLE
	direction = -d_infinity()
	y = d_nextafter(x, direction)
	write(*,3) x
 3	format(' Starting from ', 1PE23.16,
     -   ', the next representable number ')
	write(*,4) direction, y
 4	format('    towards ', F4.1, ' is ', 1PE23.16)

	x = MINDOUBLE
	direction = 1.0d0

	y = d_nextafter(x, direction)
	write(*,3) x
	write(*,4) direction, y

	x = 2.0d0
	scale = 3
	y = d_scalbn(x, scale)
	write (*,5) x, scale, y
 5	format(' Scaling ', F4.1, ' by 2**', I1, ' is ', F4.1)

	print *
	print *, 'SINGLE PRECISION EXAMPLES:'
	print *

	r = 32.0
	i = ir_ilogb(r)
	write (*,1) r, i

	r = -5.5
	i = ir_signbit(r)
	print *, 'The sign bit of ', r, ' is ', i

	r = -5.5
	s = 12.4
	t = r_copysign(r,s)
	write (*,2) r, s, t

	r = r_min_subnormal()
	r_direction = -r_infinity()
	s = r_nextafter(r, r_direction)
	write(*,3) r
	write(*,4) r_direction, s

	r = r_min_subnormal()
	r_direction = 1.0e0
	s = r_nextafter(r, r_direction)
	write(*,3) r
	write(*,4) r_direction, s

	r = 2.0
	scale = 3
	s = r_scalbn(r, scale)
	write (*,5) r, scale, y

	print *

	end

The output from this program looks like: 

DOUBLE PRECISION EXAMPLES:

The base 2 exponent of 32.0 is  5
-5.5 was given the sign of 12.4 and is now  5.5
The sign bit of    -5.5000000000000 is   1
Starting from  4.9406564584124654-324, the next representable
   number towards -Inf is  0.0000000000000000E+00
Starting from  4.9406564584124654-324, the next representable
   number towards  1.0 is  9.8813129168249309-324
Scaling  2.0 by 2**3 is 16.0

SINGLE PRECISION EXAMPLES:

The base 2 exponent of 32.0 is  5
The sign bit of    -5.50000 is   1
-5.5 was given the sign of 12.4 and is now  5.5
Starting from  1.4012984643248171E-45, the next representable
   number towards -Inf is  0.0000000000000000E+00
Starting from  1.4012984643248171E-45, the next representable
   number towards  1.0 is  2.8025969286496341E-45
Scaling  2.0 by 2**3 is 16.0

Note:IEEE floating-point exception flags raised: 
   Inexact; Underflow;
See the Numerical Computation Guide, ieee_flags(3M)

ieee_values

This C program calls several of the ieee_values(3m) functions: 

#include <math.h>
#include <sunmath.h>

main() 
{ 
	double x; 
	float r; 

	x = quiet_nan(0); 
	printf("quiet NaN: %.16e = %08x %08x \n", 
		x, ((int *) &x)[0], ((int *) &x)[1]); 

	x = nextafter(max_subnormal(), 0.0);                
	printf("nextafter(max_subnormal,0) = %.16e\n",x);
	printf("                           = %08x %08x\n", 
		((int *) &x)[0], ((int *) &x)[1]); 

	r = min_subnormalf(); 
	printf("single precision min subnormal = %.8e = %08x\n",
		r, ((int *) &r)[0]); 

	exit(0); 
}

When linking, you use both -lm and -lsunmath.

The output is similar to (SPARC): 

quiet NaN: NaN = 7fffffff ffffffff 
nextafter(max_subnormal,0) = 2.2250738585072004e-308 
                           = 000fffff fffffffe 
single precision min subnormal = 1.40129846e-45 = 00000001

(Intel and PowerPC) The output is similar to: 

quiet NaN: NaN = ffffffff 7fffffff 
nextafter(max_subnormal,0) = 2.2250738585072004e-308 
                           = fffffffe 000fffff
single precision min subnormal = 1.40129846e-45 = 00000001

FORTRAN programs that use ieee_values functions should take care to declare the function's type (SPARC): 

	program print_ieee_values
c
c the purpose of the implicit statements is to insure
c that the f77_floatingpoint pseudo-instrinsic
c functions are declared with the correct type
c
	implicit real*16 (q)
	implicit double precision (d)
	implicit real (r)
	real*16 z, zero, one
	double precision x
	real r
c
	zero = 0.0
	one = 1.0
	z = q_nextafter(zero, one)
	x = d_infinity()
	r = r_max_normal()
c
	print *, z
	print *, x
	print *, r
c
	end

(Intel) Intel does not support quadruple precision; real*16. To run the FORTRAN example on Intel, delete the real*16 declarations, the assignment of zero and one, the q_nextafter() function, and the statement to print z. 

	implicit real*16 (q)
    ...
	real*16 z, zero, one
    ...
	zero = 0.0
	one = 1.0
	z = q_nextafter(zero, one)
    ...
	print *, z

The output is similar to (SPARC):

     6.4751751194380251109244389582276466-4966
  Infinity              
     3.40282E+38

ieee_flags -- Rounding Direction

The following example demonstrates how to set the rounding mode to
round towards zero:

#include <math.h>
main() {
	int             i;
	double          x, y;
	char           *out_1, *out_2, *dummy;

	/* get prevailing rounding direction */
	i = ieee_flags("get", "direction", "", &out_1);

	x = sqrt(.5);
	printf("With rounding direction %s, \n", out_1);
	printf("sqrt(.5) = 0x%08x 0x%08x = %16.15e\n",
	       ((int *) &x)[0], ((int *) &x)[1], x);

	/* set rounding direction */
	if (ieee_flags("set", "direction", "tozero", &dummy) != 0)
		printf("Not able to change rounding direction!\n");
	i = ieee_flags("get", "direction", "", &out_2);

	x = sqrt(.5);
	/*
	 * restore original rounding direction before printf, since
	 * printf is also affected by the current rounding direction
	 */
	if (ieee_flags("set", "direction", out_1, &dummy) != 0)
		printf("Not able to change rounding direction!\n");
	printf("\nWith rounding direction %s,\n", out_2);
	printf("sqrt(.5) = 0x%08x 0x%08x = %16.15e\n",
	       ((int *) &x)[0], ((int *) &x)[1], x);

	exit(0); 
}

(SPARC) The output of this short program shows the effects of rounding towards zero: 

demo% cc rounding_direction.c -lm
demo% a.out 
With rounding direction nearest, 
sqrt(.5) = 0x3fe6a09e 0x667f3bcd  = 7.071067811865476e-01 

With rounding direction tozero, 
sqrt(.5) = 0x3fe6a09e 0x667f3bcc  = 7.071067811865475e-01 
demo%

(Intel and PowerPC) The output of this short program shows the effects of rounding towards zero: 

demo% cc rounding_direction.c -lm
demo% a.out 
With rounding direction nearest, 
sqrt(.5) = 0x667f3bcd 0x3fe6a09e  = 7.071067811865476e-01 

With rounding direction tozero, 
sqrt(.5) = 0x667f3bcc 0x3fe6a09e  = 7.071067811865475e-01 
demo%

Set rounding direction towards zero from a FORTRAN program: 

	program ieee_flags_demo
	character*16 out

	i = ieee_flags('set', 'direction', 'tozero', out)
	if (i.ne.0) print *, 'not able to set rounding direction'

	i = ieee_flags('get', 'direction', '', out)
	print *, 'Rounding direction is: ', out

	end

The output is as follows: 

 Rounding direction is: tozero
 Note: Rounding direction toward zero. 
 See the Numerical Computation Guide, ieee_flags(3M)

Exceptions and Exception Handling

ieee_flags -- Accrued Exceptions

Generally, a user program examines or clears the accrued exception bits. Here is a C program that examines the accrued exception flags: 

#include <sunmath.h>
#include <sys/ieeefp.h>

main()
{
	int code, inexact, division, underflow, overflow, invalid;
	double x;
	char *out;

	/* cause an underflow exception */
	x = max_subnormal() / 2.0;

	/* this statement insures that the previous */
	/* statement is not optimized away          */
	printf("x = %g\n",x);

	/* find out which exceptions are raised */
	code = ieee_flags("get", "exception", "", &out);

	/* decode the return value */
	inexact =      (code >> fp_inexact)     & 0x1;
	underflow =    (code >> fp_underflow)   & 0x1;
	division =     (code >> fp_division)    & 0x1;
	overflow =     (code >> fp_overflow)    & 0x1;
	invalid =      (code >> fp_invalid)     & 0x1;

      /* "out" is the raised exception with the highest priority */
	printf(" Highest priority exception is: %s\n", out);

	/* The value 1 means the exception is raised, */
	/* 0 means it isn't.                          */
	printf("%d %d %d %d %d\n", invalid, overflow, division,
		underflow, inexact);

	exit(0);
}

The output from running this program: 

demo% a.out 
x = 1.11254e-308
 Highest priority exception is: underflow
0 0 0 1 1
 Note:IEEE floating-point exception flags raised:   
    Inexact;  Underflow; 
 See the Numerical Computation Guide, ieee_flags(3M)

The same can be done from FORTRAN:

/*
A FORTRAN example that:
     
      *  causes an underflow exception
      *  uses ieee_flags to determine which exceptions are raised
      *  decodes the integer value returned by ieee_flags 
      *  clears all outstanding exceptions

Remember to save this program in a file with the suffix .F, so that
the c preprocessor is invoked to bring in the header file
f77_floatingpoint.h.  
*/

#include <f77_floatingpoint.h>

      program decode_accrued_exceptions
      double precision x
      integer accrued, inx, div, under, over, inv
      character*16 out
      double precision d_max_subnormal

c Cause an underflow exception
      x = d_max_subnormal() / 2.0

c Find out which exceptions are raised
      accrued = ieee_flags('get', 'exception', '', out)

c Decode the value returned by ieee_flags using bit-shift intrinsics
      inx   = and(rshift(accrued, fp_inexact)  , 1)
      under = and(rshift(accrued, fp_underflow), 1)
      div   = and(rshift(accrued, fp_division) , 1)
      over  = and(rshift(accrued, fp_overflow) , 1)
      inv   = and(rshift(accrued, fp_invalid)  , 1)

c The exception with the highest priority is returned in "out"
      print *, "Highest priority exception is ", out

c The value 1 means the exception is raised; 0 means it is not
      print *, inv, over, div, under, inx

c Clear all outstanding exceptions
      i = ieee_flags('clear', 'exception', 'all', out)

      end

The output is as follows: 

 Highest priority exception is underflow       
   0  0  0  1  1

(Intel) While it is unusual for a user program to set exception flags, it can be done. This is demonstrated in the following C example. 

#include <sys/ieeefp.h>
/*
 * from <sys/ieeefp.h>, the exceptions according to bit number are
 *   fp_invalid      = 0,
 *   fp_denormalized = 1,
 *   fp_division     = 2,
 *   fp_overflow     = 3,
 *   fp_underflow    = 4,
 *   fp_inexact      = 5
 */

main()
{
	int             code;
	char           *out;

	if (ieee_flags("clear", "exception", "all", &out) != 0)
	   printf("could not clear exceptions\n");
	if (ieee_flags("set", "exception", "division", &out) != 0)
	   printf("could not set exception\n");
	code = ieee_flags("get", "exception", "", &out);
	printf("out is: %s , fp exception code is: %X \n",
		out, code);

	exit(0);
}

For a SPARC machine, the exception flags have the following definitions: 

/*
 *from <sys/ieeefp.h>, the exceptions according to bit number are
 *   fp_inexact      = 0,
 *   fp_division     = 1,
 *   fp_underflow    = 2,
 *   fp_overflow     = 3,
 *   fp_invalid      = 4
 */

(Intel) The output of the Intel C example is: 

out is: division , fp exception code is: 2

ieee_handler -- Trapping Exceptions

Note - The examples below only apply to the Solaris operating environment.
Here is a FORTRAN program that installs a signal handler to locate an exception (for SPARC and PowerPC systems only): 

      program demo

c declare signal handler function
      external fp_exc_hdl
      double precision d_min_normal
      double precision x

c set up signal handler
      i = ieee_handler('set', 'common', fp_exc_hdl)
      if (i.ne.0) print *, 'ieee trapping not supported here'

c cause an underflow exception (it will not be trapped)
      x = d_min_normal() / 13.0
      print *, 'd_min_normal() / 13.0 = ', x

c cause an overflow exception
c the value printed out is unrelated to the result
      x = 1.0d300*1.0d300
      print *, '1.0d300*1.0d300 = ', x

      end 
c
c the floating-point exception handling function
c
      integer function fp_exc_hdl(sig, sip, uap)

      integer sig, code, addr
      character label*16
c
c The structure /siginfo/ is a translation of siginfo_t from
c 
c
      structure /fault/
         integer address
      end structure
     
      structure /siginfo/
         integer si_signo
         integer si_code
         integer si_errno
         record /fault/ fault
      end structure
  
      record /siginfo/ sip
 
c See  for list of FPE codes
c Figure out the name of the SIGFPE
      code = sip.si_code
      if (code.eq.3) label = 'division'
      if (code.eq.4) label = 'overflow'
      if (code.eq.5) label = 'underflow'
      if (code.eq.6) label = 'inexact'
      if (code.eq.7) label = 'invalid'
      addr = sip.fault.address

c Print information about the signal that happened
      write (*,77) code, label, addr
 77   format ('floating-point exception code ', i2, ',',
     *       a17, ',', ' at address ', z8 )

      end

The output is:

 d_min_normal() / 13.0 =     1.7115952757748-309
floating-point exception code  4, overflow        , at address
1131C
 1.0d300*1.0d300 =     1.0000000000000+300
 Note: IEEE floating-point exception flags raised: 
    Inexact;  Underflow; 
 IEEE floating-point exception traps enabled:
    overflow; division by zero; invalid operation; 
 See the Numerical Computation Guide, ieee_handler(3M),

    ieee_flags(3M)

(SPARC) Here is a somewhat more complex C example: 

/*
 * Generate the 5 IEEE exceptions: invalid, division,
 * overflow, underflow and inexact.
 * Trap on any floating point exception, print a message,
 * and continue.
 * Note that one could also inquire about raised exceptions by
 *    i = ieee_flags("get","exception","",&out);
 * where out will contain the name of the highest exception
 * raised, and i can be decoded to find out about all the
 * exceptions raised.
 */

#include <sunmath.h>
#include <signal.h>
#include <siginfo.h>
#include <ucontext.h>

extern void trap_all_fp_exc(int sig, siginfo_t *sip,
	ucontext_t *uap);

main()
{
	double          x, y, z;
	char           *out;

	/*
	 * Use ieee_handler to establish "trap_all_fp_exc"
	 * as the signal handler to use whenever any floating
	 * point exception occurs.
	 */

	if (ieee_handler("set", "all", trap_all_fp_exc) != 0) {
		printf(" IEEE trapping not supported here.\n");
	}

	/* disable trapping (uninteresting) inexact exceptions */
	if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
		printf("Trap handler for inexact not cleared.\n");

	/* raise invalid */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("1. Invalid: signaling_nan(0) * 2.5\n");
	x = signaling_nan(0);
	y = 2.5;
	z = x * y;

	/* raise division */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("2. Div0: 1.0 / 0.0\n");
	x = 1.0;
	y = 0.0;
	z = x / y;

	/* raise overflow */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("3. Overflow: -max_normal() - 1.0e294\n");
	x = -max_normal();
	y = -1.0e294;
	z = x + y;

	/* raise underflow */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("4. Underflow: min_normal() * min_normal()\n");
	x = min_normal();
	y = x;
	z = x * y;

	/* enable trapping on inexact exception */
	if (ieee_handler("set", "inexact", trap_all_fp_exc) != 0)
		printf("Could not set trap handler for inexact.\n");

	/* raise inexact */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("5. Inexact: 2.0 / 3.0\n");
	x = 2.0;
	y = 3.0;
	z = x / y;

	/* don't trap on inexact */
	if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
		printf(" could not reset inexact trap\n");

	/* check that we're not trapping on inexact anymore */
	if (ieee_flags("clear", "exception", "all", &out) != 0)
		printf(" could not clear exceptions\n");
	printf("6. Inexact trapping disabled; 2.0 / 3.0\n");
	x = 2.0;
	y = 3.0;
	z = x / y;

	/* find out if there are any outstanding exceptions */
	ieee_retrospective_();

	/* exit gracefully */
	exit(0);
}

void trap_all_fp_exc(int sig, siginfo_t *sip, ucontext_t *uap) {
	char           *label = "undefined";

/* see /usr/include/sys/machsig.h for SIGFPE codes */
	switch (sip->si_code) {
	case FPE_FLTRES:
		label = "inexact";
		break;
	case FPE_FLTDIV:
		label = "division";
		break;
	case FPE_FLTUND:
		label = "underflow";
		break;
	case FPE_FLTINV:
		label = "invalid";
		break;
	case FPE_FLTOVF:
		label = "overflow";
		break;
	}

	printf(
	   "   signal %d, sigfpe code %d: %s exception at address %x\n",
	   sig, sip->si_code, label, sip->_data._fault._addr);
}

The output is similar to the following:

1. Invalid: signaling_nan(0) * 2.5
   signal 8, sigfpe code 7: invalid exception at address 10da8
2. Div0: 1.0 / 0.0
   signal 8, sigfpe code 3: division exception at address 10e44
3. Overflow: -max_normal() - 1.0e294
   signal 8, sigfpe code 4: overflow exception at address 10ee8
4. Underflow: min_normal() * min_normal()
   signal 8, sigfpe code 5: underflow exception at address 10f80
5. Inexact: 2.0 / 3.0
   signal 8, sigfpe code 6: inexact exception at address 1106c
6. Inexact trapping disabled; 2.0 / 3.0
Note: IEEE floating-point exception traps enabled:  
   underflow; overflow; division by zero; invalid operation; 
See the Numerical Computation Guide, ieee_handler(3M)

(SPARC) The following program shows how you can use ieee_handler and the include files to modify the default result of certain exceptional situations: 

/*
 * Cause a division by zero exception and use the
 * signal handler to substitute MAXDOUBLE (or MAXFLOAT)
 * as the result.
 * 
 * compile with the flag -Xa
 */

#include <values.h>
#include <siginfo.h>
#include <ucontext.h>

void division_handler(int sig, siginfo_t *sip, ucontext_t *uap);

main() {
	double          x, y, z;
	float           r, s, t;
	char           *out;

	/*
	 * Use ieee_handler to establish division_handler as the
	 * signal * handler to use for the IEEE exception division.
	 */
	if (ieee_handler("set","division",division_handler)!=0) {
		printf(" IEEE trapping not supported here.\n");
	}

	/* Cause a division-by-zero exception */
	x = 1.0;
	y = 0.0;
	z = x / y;

	/*
	 * Check to see that the user-supplied value, MAXDOUBLE,
	 * is indeed substituted in place of the IEEE default
	 * value, infinity.
	 */
	printf("double precision division: %g/%g = %g \n",x,y,z);

	/* Cause a division-by-zero exception */
	r = 1.0;
	s = 0.0;
	t = r / s;

	/*
	 * Check to see that the user-supplied value, MAXFLOAT,
	 * is indeed substituted in place of the IEEE default
	 * value, infinity.
	 */
	printf("single precision division: %g/%g = %g \n",r,s,t);

	ieee_retrospective_();

	exit(0);
}

void division_handler(int sig, siginfo_t *sip, ucontext_t *uap) {
	int             inst; 
	unsigned        rd, mask, single_prec=0;
	float           f_val = MAXFLOAT;
	double          d_val = MAXDOUBLE;
	long           *f_val_p = (long *) &f_val;

	/* Get instruction that caused exception. */
	inst = *(int *) sip->_data._fault._addr;

	/*
	 * Decode the destination register. Bits 29:25 encode the
	 * destination register for any SPARC floating point
	 * instruction.
	 */
	mask = 0x1f;
	rd = (mask & (inst >> 25));

	/*
	 * Is this a single precision or double precision
	 * instruction?  Bits 5:6 encode the precision of the
	 * opcode; if bit 5 is 1, it's sp, else, dp.
	 */
	mask = 0x1;
	single_prec = (mask & (inst >> 5));

	/* put user-defined value into destination register */
	if (single_prec) {
	   uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[rd] =
	      f_val_p[0];

	} else {
	   uap->uc_mcontext.fpregs.fpu_fr.fpu_dregs[rd/2] = d_val;
	}
}

As expected, the output is:

double precision division: 1/0 = 1.79769e+308 
single precision division: 1/0 = 3.40282e+38 
Note: IEEE floating-point exception traps enabled: 
   division by zero; 
See the Numerical Computation Guide, ieee_handler(3M)

ieee_handler -- Abort on Exceptions

You can use ieee_handler to force a program to abort in case of certain floating-point exceptions:

#include <f77_floatingpoint.h>
	program abort
c
	ieeer = ieee_handler('set', 'division', SIGFPE_ABORT)
	if (ieeer .ne. 0) print *, ' ieee trapping not supported'
	r = 14.2
	s = 0.0
	r = r/s
c
	print *, 'you should not see this; system should abort'
c
	end

Miscellaneous

sigfpe -- Trapping Integer Exceptions

The previous section showed examples of using ieee_handler. In general, when there is a choice between using ieee_handler or sigfpe, the former is recommended.

Note - SIGFPE is available only in the Solaris operating environment.
(SPARC) There are instances, such as trapping integer arithmetic exceptions, when sigfpe is the handler to be used. The following example traps on integer division by zero: 

/* Generate the integer division by zero exception */

#include <siginfo.h>
#include <ucontext.h>
#include <signal.h>

void int_handler(int sig, siginfo_t *sip, ucontext_t *uap);

main() {
	int a, b, c;

/*
 * Use sigfpe(3) to establish "int_handler" as the signal handler
 * to use on integer division by zero
 */

/*
 * Integer division-by-zero aborts unless a signal
 * handler for integer division by zero is set up
 */

	sigfpe(FPE_INTDIV, int_handler);

	a = 4;
	b = 0;
	c = a / b;
	printf("%d / %d = %d\n\n", a, b, c);

	exit(0);
}

void int_handler(int sig, siginfo_t *sip, ucontext_t *uap) {
	printf("Signal %d, code %d, at addr %x\n",
	   sig, sip->si_code, sip->_data._fault._addr);

/*
 * Increment program counter; ieee_handler does this by default,
 * but here we have to use sigfpe to set up the signal handler
 * for int div by 0
 */
	uap->uc_mcontext.gregs[REG_PC] =
	   uap->uc_mcontext.gregs[REG_nPC];
}

Calling FORTRAN from C

Here is a simple example of a C driver calling FORTRAN subroutines. Refer to the appropriate C and FORTRAN manuals for more information on working with C and FORTRAN. The following is the C driver (save it in a file named driver.c): 

/*
 * a demo program that shows:
 * 
 * 1. how to call f77 subroutine from C, passing an array argument
 * 2. how to call single precision f77 function from C
 * 3. how to call double precision f77 function from C
 */

extern int      demo_one_(double *);
extern float    demo_two_(float *);
extern double   demo_three_(double *);

main()
{
	double          array[3][4];
	float           f, g;
	double          x, y;
	int             i, j;

	for (i = 0; i < 3; i++)
		for (j = 0; j < 4; j++)
			array[i][j] = i + 2*j;
	g = 1.5;
	y = g;

	/* pass an array to a fortran function (print the array) */
	demo_one_(&array[0][0]);
	printf(" from the driver\n");
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 4; j++)
			printf("    array[%d][%d] = %e\n",
			   i, j, array[i][j]);
		printf("\n");
	}

	/* call a single precision fortran function */
	f = demo_two_(&g);
	printf(
	   " f = sin(g) from a single precision fortran function\n");
	printf("    f, g: %8.7e, %8.7e\n", f, g);
	printf("\n");

	/* call a double precision fortran function */
	x = demo_three_(&y);
	printf(
	   " x = sin(y) from a double precision fortran function\n");
	printf("    x, y: %18.17e, %18.17e\n", x, y);

	ieee_retrospective_();
	exit(0);
}

Save the FORTRAN subroutines in a file named drivee.f: 

      subroutine demo_one(array)
      double precision array(4,3)
      print *, 'from the fortran routine:'
      do 10 i =1,4
         do 20 j = 1,3
            print *, '   array[', i, '][', j, '] = ', array(i,j)
 20      continue
	 print *
 10   continue
      return
      end

      real function demo_two(number)
      real number
      demo_two = sin(number)
      return
      end
     
      double precision function demo_three(number)
      double precision number
      demo_three = sin(number)
      return 
      end

Then, perform the compilation and linking:

cc -c driver.c
f77 -c drivee.f
	demo_one:
	demo_two:
	demo_three:
f77 -o driver driver.o drivee.o

The output looks like this:

 from the fortran routine:
    array[  1][  1] =   0.
    array[  1][  2] =     1.0000000000000
    array[  1][  3] =     2.0000000000000

    array[  2][  1] =     2.0000000000000
    array[  2][  2] =     3.0000000000000
    array[  2][  3] =     4.0000000000000

    array[  3][  1] =     4.0000000000000
    array[  3][  2] =     5.0000000000000
    array[  3][  3] =     6.0000000000000

    array[  4][  1] =     6.0000000000000
    array[  4][  2] =     7.0000000000000
    array[  4][  3] =     8.0000000000000

 from the driver
    array[0][0] = 0.000000e+00
    array[0][1] = 2.000000e+00
    array[0][2] = 4.000000e+00
    array[0][3] = 6.000000e+00

    array[1][0] = 1.000000e+00
    array[1][1] = 3.000000e+00
    array[1][2] = 5.000000e+00
    array[1][3] = 7.000000e+00

    array[2][0] = 2.000000e+00
    array[2][1] = 4.000000e+00
    array[2][2] = 6.000000e+00
    array[2][3] = 8.000000e+00

 f = sin(g) from a single precision fortran function
    f, g: 9.9749500e-01, 1.5000000e+00

 x = sin(y) from a double precision fortran function
    x, y: 9.97494986604054446e-01, 1.50000000000000000e+00

A Few Useful Debugging Commands

Table A-1 shows examples of debugging commands for the SPARC architecture

Table  A-1 Some Debugging Commands (SPARC)  

Action
dbx
adb

Set breakpoint

at function

at line number

at absolute address

at relative address

stop in myfunct

stop at 29

myfunct:b

23a8:b

main+0x40:b

Run until breakpoint met

run

:r

Examine source code

list

<pc,10?ia

Examine a fp register

IEEE single precision

decimal equivalent

IEEE double precision

decimal equivalent

examine &$f0/X

print $f0

examine &$f0/2X

print $f0f1

<f0=X

<f0=f

<f0=X; <f1=X

<f0=F

Examine all fp registers

examine &$f0/32X

$x for f0-f15

$X for f16-f31

Examine all registers

examine &$g0/64X

$r; $x; $X

Examine fp status register

examine &$fsr/X

<fsr=X

Put single precision 1.0 in f0

Put double prec 1.0 in f0/f1

assign $f0=1.0

assign $f0f1=1.0

3f800000>f0

3ff00000>f0; 0>f1

Continue execution

cont

:c

Single step

step (or next)

:s

Exit the debugger

quit

$q

Table A-2 shows examples of debugging commands for the Intel architecture.

Table  A-2 Some Debugging Commands (Intel)  

Action
dbx
adb

Set breakpoint

at function

at line number

at absolute address

at relative address

stop in myfunct

stop at 29

myfunct:b

23a8:b

main+0x40:b

Run until breakpoint met

run

:r

Examine source code

list

<pc,10?ia

Examine fp registers

print $st0

...

print $st7

$x

Examine all registers

examine &$gs/19X

$r

Examine fp status register

examine &$fstat/X

<fstat=X

or $x

Continue execution

cont

:c

Single step

step (or next)

:s

Exit the debugger

quit

$q

Table A-3 shows examples of debugging commands for the PowerPC architecture.

Table  A-3 Some Debugging Commands (PowerPC)  

Action
dbx
adb

Set breakpoint

at function

at line number

at absolute address

at relative address

stop in myfunct

stop at 29

myfunct:b

23a8:b

main+0x40:b

Run until breakpoint met

run

:r

Examine source code

list

<pc,10?ia

Examine a fp register

IEEE double precision

decimal equivalent

examine &$f0/2X

print $f0

<f0=F

Examine all fp registers

examine &$f0/64X

$x

Examine all registers

examine &$r0/159X

$r; $x

Examine fp status register

examine &$fpscr/X

<fpscr=X

Put single precision 1.0 in f0

assign $f0 = 1.0

Continue execution

cont

:c

Single step

step (or next)

:s

Exit the debugger

quit

$q

The following examples show two ways to set a breakpoint at the beginning of the code corresponding to a routine myfunction in adb. First you can say: 

myfunction:b

Second, you can determine the absolute address that corresponds to the beginning of the piece of code corresponding to myfunction, and then set a break at that absolute address: 

myfunction=X 
			23a8 
23a8:b

The main in FORTRAN programs is known as MAIN_ to adb. To set a breakpoint at a FORTRAN main in adb: 

   MAIN_:b
When examining the contents of floating-point registers, the hex value shown by the dbx command x &$f0/X is the base-16 representation, not the number's decimal representation. For SPARC, the adb commands $x and $X display both the hexadecimal representation, and the decimal value. For Intel and PowerPC, the adb command $x displays only the decimal value. For SPARC, the double precision values show the decimal value next to the odd-numbered register.

Because the operating system disables the floating-point unit until it is first used by a process, you cannot modify the floating-point registers until they have been accessed by the program being debugged.

(SPARC) When displaying floating point numbers, you should keep in mind that the size of registers is 32 bits, a single precision floating-point number occupies 32 bits (hence it fits in one register), and double precision floating-point numbers occupy 64 bits (therefore two registers are used to hold a double precision number). In the hexadecimal representation 32 bits correspond to 8 digit numbers. In the following snapshot of FPU registers displayed with adb, the display is organized as follows:

<name of fpu register> <IEEE hex value> <single precision> <double precision>

(SPARC) The third column holds the single precision decimal interpretation of the hexadecimal pattern shown in the second column. The fourth column interprets pairs of registers. For example, the fourth column of the f11 line interprets f10 and f11 as a 64-bit IEEE double precision number.

(SPARC) Because f10 and f11 are used to hold a double precision value, the interpretation (on the f10 line) of the first 32 bits of that value, 7ff00000, as +NaN, is irrelevant. The interpretation of all 64 bits, 7ff00000 00000000, as +Infinity, happens to be the meaningful translation.

(SPARC) The adb command $x, that was used to display the first 16 floating-point data registers, also displayed fsr (the floating-point status register): 

$x 
fsr	   40020
f0	400921fb	     +2.1426990e+00 
f1	54442d18	     +3.3702806e+12	     +3.1415926535897931e+00 
f2	       2	     +2.8025969e-45 
f3	       0	     +0.0000000e+00	     +4.2439915819305446e-314 
f4	40000000	     +2.0000000e+00 
f5	       0	     +0.0000000e+00	     +2.0000000000000000e+00 
f6	3de0b460	     +1.0971904e-01 
f7	       0	     +0.0000000e+00	     +1.2154188766544394e-10 
f8	3de0b460	     +1.0971904e-01 
f9	       0	     +0.0000000e+00	     +1.2154188766544394e-10 
f10	7ff00000	     +NaN 
f11	       0	     +0.0000000e+00	     +Infinity 
f12	ffffffff	     -NaN 
f13	ffffffff	     -NaN	                -NaN 
f14	ffffffff	     -NaN 
f15	ffffffff	     -NaN	                -NaN

(Intel) The corresponding output on Intel looks like: 

$x
80387 chip is present.
cw      0x137f
sw      0x3920
cssel 0x17  ipoff 0x2d93                datasel 0x1f  dataoff 0x5740
 
 st[0]  +3.24999988079071044921875 e-1            VALID
 st[1]  +5.6539133243479549034419688 e73          EMPTY
 st[2]  +2.0000000000000008881784197              EMPTY
 st[3]  +1.8073218308070440556016047 e-1          EMPTY
 st[4]  +7.9180300235748291015625 e-1             EMPTY
 st[5]  +4.201639036693904927233234 e-13          EMPTY
 st[6]  +4.201639036693904927233234 e-13          EMPTY
 st[7]  +2.7224999213218694649185636              EMPTY

Note - (Intel) cw is the control word; sw is the status word.


Previous Next Contents Index Doc Set Home