n_fmod.c revision 1.3 
1/*	$NetBSD: n_fmod.c,v 1.3 1998/10/20 02:26:11 matt Exp $	*/
2/*
3 * Copyright (c) 1989, 1993
4 *	The Regents of the University of California.  All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 *    notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 *    notice, this list of conditions and the following disclaimer in the
13 *    documentation and/or other materials provided with the distribution.
14 * 3. All advertising materials mentioning features or use of this software
15 *    must display the following acknowledgement:
16 *	This product includes software developed by the University of
17 *	California, Berkeley and its contributors.
18 * 4. Neither the name of the University nor the names of its contributors
19 *    may be used to endorse or promote products derived from this software
20 *    without specific prior written permission.
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32 * SUCH DAMAGE.
33 */
34
35#ifndef lint
36#if 0
37static char sccsid[] = "@(#)fmod.c	8.1 (Berkeley) 6/4/93";
38#endif
39#endif /* not lint */
40
41#include "mathimpl.h"
42
43/* fmod.c
44 *
45 * SYNOPSIS
46 *
47 *    #include <math.h>
48 *    double fmod(double x, double y)
49 *
50 * DESCRIPTION
51 *
52 *    The fmod function computes the floating-point remainder of x/y.
53 *
54 * RETURNS
55 *
56 *    The fmod function returns the value x-i*y, for some integer i
57 * such that, if y is nonzero, the result has the same sign as x and
58 * magnitude less than the magnitude of y.
59 *
60 * On a VAX or CCI,
61 *
62 *    fmod(x,0) traps/faults on floating-point divided-by-zero.
63 *
64 * On IEEE-754 conforming machines with "isnan()" primitive,
65 *
66 *    fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned.
67 *
68 */
69#if !defined(__vax__) && !defined(tahoe)
70extern int isnan(),finite();
71#endif	/* !defined(__vax__) && !defined(tahoe) */
72
73#ifdef TEST_FMOD
74static double
75_fmod(x,y)
76#else	/* TEST_FMOD */
77double
78fmod(x,y)
79#endif	/* TEST_FMOD */
80double x,y;
81{
82	int ir,iy;
83	double r,w;
84
85	if (y == (double)0
86#if !defined(__vax__) && !defined(tahoe)	/* per "fmod" manual entry, SunOS 4.0 */
87		|| isnan(y) || !finite(x)
88#endif	/* !defined(__vax__) && !defined(tahoe) */
89	    )
90	    return (x*y)/(x*y);
91
92	r = fabs(x);
93	y = fabs(y);
94	(void)frexp(y,&iy);
95	while (r >= y) {
96		(void)frexp(r,&ir);
97		w = ldexp(y,ir-iy);
98		r -= w <= r ? w : w*(double)0.5;
99	}
100	return x >= (double)0 ? r : -r;
101}
102
103#ifdef TEST_FMOD
104extern long random();
105extern double fmod();
106
107#define	NTEST	10000
108#define	NCASES	3
109
110static int nfail = 0;
111
112static void
113doit(x,y)
114double x,y;
115{
116	double ro = fmod(x,y),rn = _fmod(x,y);
117	if (ro != rn) {
118		(void)printf(" x    = 0x%08.8x %08.8x (%24.16e)\n",x,x);
119		(void)printf(" y    = 0x%08.8x %08.8x (%24.16e)\n",y,y);
120		(void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro);
121		(void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn);
122		(void)printf("\n");
123	}
124}
125
126main()
127{
128	register int i,cases;
129	double x,y;
130
131	srandom(12345);
132	for (i = 0; i < NTEST; i++) {
133		x = (double)random();
134		y = (double)random();
135		for (cases = 0; cases < NCASES; cases++) {
136			switch (cases) {
137			case 0:
138				break;
139			case 1:
140				y = (double)1/y; break;
141			case 2:
142				x = (double)1/x; break;
143			default:
144				abort(); break;
145			}
146			doit(x,y);
147			doit(x,-y);
148			doit(-x,y);
149			doit(-x,-y);
150		}
151	}
152	if (nfail)
153		(void)printf("Number of failures: %d (out of a total of %d)\n",
154			nfail,NTEST*NCASES*4);
155	else
156		(void)printf("No discrepancies were found\n");
157	exit(0);
158}
159#endif	/* TEST_FMOD */
160