diff options
author | das <das@FreeBSD.org> | 2008-06-19 22:39:53 +0000 |
---|---|---|
committer | das <das@FreeBSD.org> | 2008-06-19 22:39:53 +0000 |
commit | 9e4d306f6fe863ceecc097ac5554fee53ab1d7e5 (patch) | |
tree | 4708ca4fbfb23bf4b3fdff4cfb4fe046c81a47d5 | |
parent | 4f152d47fa3740ab3f9b056d153678d594931414 (diff) | |
download | FreeBSD-src-9e4d306f6fe863ceecc097ac5554fee53ab1d7e5.zip FreeBSD-src-9e4d306f6fe863ceecc097ac5554fee53ab1d7e5.tar.gz |
Implement fmodl.
Document fmodl and fix some errors in the fmod manpage.
-rw-r--r-- | lib/msun/Makefile | 4 | ||||
-rw-r--r-- | lib/msun/Symbol.map | 1 | ||||
-rw-r--r-- | lib/msun/man/fmod.3 | 46 | ||||
-rw-r--r-- | lib/msun/src/e_fmodl.c | 149 | ||||
-rw-r--r-- | lib/msun/src/math.h | 2 |
5 files changed, 177 insertions, 25 deletions
diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 3ffe7e4..cc97b08 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -79,7 +79,7 @@ SYMBOL_MAPS= ${SYM_MAPS} COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c .if ${LDBL_PREC} != 53 # If long double != double use these; otherwise, we alias the double versions. -COMMON_SRCS+= e_hypotl.c e_remainderl.c e_sqrtl.c \ +COMMON_SRCS+= e_fmodl.c e_hypotl.c e_remainderl.c e_sqrtl.c \ k_cosl.c k_sinl.c k_tanl.c \ s_ceill.c s_cosl.c s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \ s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \ @@ -145,7 +145,7 @@ MLINKS+=floor.3 floorf.3 floor.3 floorl.3 MLINKS+=fma.3 fmaf.3 fma.3 fmal.3 MLINKS+=fmax.3 fmaxf.3 fmax.3 fmaxl.3 \ fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3 -MLINKS+=fmod.3 fmodf.3 +MLINKS+=fmod.3 fmodf.3 fmod.3 fmodl.3 MLINKS+=hypot.3 cabs.3 hypot.3 cabsf.3 hypot.3 cabsl.3 \ hypot.3 hypotf.3 hypot.3 hypotl.3 MLINKS+=ieee_test.3 scalb.3 ieee_test.3 scalbf.3 diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index 3f73c92..526d7e3 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -208,4 +208,5 @@ FBSD_1.1 { csqrtl; remquol; remainderl; + fmodl; }; diff --git a/lib/msun/man/fmod.3 b/lib/msun/man/fmod.3 index 9a73675..670efdb 100644 --- a/lib/msun/man/fmod.3 +++ b/lib/msun/man/fmod.3 @@ -28,12 +28,13 @@ .\" from: @(#)fmod.3 5.1 (Berkeley) 5/2/91 .\" $FreeBSD$ .\" -.Dd May 2, 1991 +.Dd June 19, 2008 .Dt FMOD 3 .Os .Sh NAME .Nm fmod , -.Nm fmodf +.Nm fmodf , +.Nm fmodl .Nd floating-point remainder functions .Sh LIBRARY .Lb libm @@ -43,41 +44,44 @@ .Fn fmod "double x" "double y" .Ft float .Fn fmodf "float x" "float y" +.Ft long double +.Fn fmodl "long double x" "long double y" .Sh DESCRIPTION The -.Fn fmod -and the -.Fn fmodf +.Fn fmod , +.Fn fmodf , +and +.Fn fmodl functions compute the floating-point remainder of .Fa x Ns / Fa y . .Sh RETURN VALUES -The -.Fn fmod -and the -.Fn fmodf +If +.Fa y +is non-zero, the +.Fn fmod , +.Fn fmodf , +and +.Fn fmodl functions return the value .Sm off .Fa x - Em i * Fa y , .Sm on for some integer -.Em i -such that, if -.Fa y -is non-zero, the result has the same sign as +.Em i , +such that the result has the same sign as .Fa x and magnitude less than the magnitude of .Fa y . If .Fa y -is zero, whether a domain error occurs or the -.Fn fmod -and the -.Fn fmodf -function returns zero is implementation-defined. +is zero, a \*(Na is produced. .Sh SEE ALSO .Xr math 3 .Sh STANDARDS The -.Fn fmod -function conforms to -.St -isoC . +.Fn fmod , +.Fn fmodf , +and +.Fn fmodl +functions conform to +.St -isoC-99 . diff --git a/lib/msun/src/e_fmodl.c b/lib/msun/src/e_fmodl.c new file mode 100644 index 0000000..c983000 --- /dev/null +++ b/lib/msun/src/e_fmodl.c @@ -0,0 +1,149 @@ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <float.h> +#include <stdint.h> + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +#define BIAS (LDBL_MAX_EXP - 1) + +#if LDBL_MANL_SIZE > 32 +typedef uint64_t manl_t; +#else +typedef uint32_t manl_t; +#endif + +#if LDBL_MANH_SIZE > 32 +typedef uint64_t manh_t; +#else +typedef uint32_t manh_t; +#endif + +/* + * These macros add and remove an explicit integer bit in front of the + * fractional mantissa, if the architecture doesn't have such a bit by + * default already. + */ +#ifdef LDBL_IMPLICIT_NBIT +#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) +#define HFRAC_BITS LDBL_MANH_SIZE +#else +#define SET_NBIT(hx) (hx) +#define HFRAC_BITS (LDBL_MANH_SIZE - 1) +#endif + +#define MANL_SHIFT (LDBL_MANL_SIZE - 1) + +static const long double one = 1.0, Zero[] = {0.0, -0.0,}; + +/* + * fmodl(x,y) + * Return x mod y in exact arithmetic + * Method: shift and subtract + * + * Assumptions: + * - The low part of the mantissa fits in a manl_t exactly. + * - The high part of the mantissa fits in an int64_t with enough room + * for an explicit integer bit in front of the fractional bits. + */ +long double +fmodl(long double x, long double y) +{ + union IEEEl2bits ux, uy; + int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */ + manh_t hy; + manl_t lx,ly,lz; + int ix,iy,n,sx; + + ux.e = x; + uy.e = y; + sx = ux.bits.sign; + + /* purge off exception values */ + if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */ + (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */ + (uy.bits.exp == BIAS + LDBL_MAX_EXP && + ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */ + return (x*y)/(x*y); + if(ux.bits.exp<=uy.bits.exp) { + if((ux.bits.exp<uy.bits.exp) || + (ux.bits.manh<=uy.bits.manh && + (ux.bits.manh<uy.bits.manh || + ux.bits.manl<uy.bits.manl))) { + return x; /* |x|<|y| return x or x-y */ + } + if(ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl) { + return Zero[sx]; /* |x|=|y| return x*0*/ + } + } + + /* determine ix = ilogb(x) */ + if(ux.bits.exp == 0) { /* subnormal x */ + ux.e *= 0x1.0p512; + ix = ux.bits.exp - (BIAS + 512); + } else { + ix = ux.bits.exp - BIAS; + } + + /* determine iy = ilogb(y) */ + if(uy.bits.exp == 0) { /* subnormal y */ + uy.e *= 0x1.0p512; + iy = uy.bits.exp - (BIAS + 512); + } else { + iy = uy.bits.exp - BIAS; + } + + /* set up {hx,lx}, {hy,ly} and align y to x */ + hx = SET_NBIT(ux.bits.manh); + hy = SET_NBIT(uy.bits.manh); + lx = ux.bits.manl; + ly = uy.bits.manl; + + /* fix point fmod */ + n = ix - iy; + + while(n--) { + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;} + else { + if ((hz|lz)==0) /* return sign(x)*0 */ + return Zero[sx]; + hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; + } + } + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz>=0) {hx=hz;lx=lz;} + + /* convert back to floating value and restore the sign */ + if((hx|lx)==0) /* return sign(x)*0 */ + return Zero[sx]; + while(hx<(1U<<HFRAC_BITS)) { /* normalize x */ + hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx; + iy -= 1; + } + ux.bits.manh = hx; /* The mantissa is truncated here if needed. */ + ux.bits.manl = lx; + if (iy < LDBL_MIN_EXP) { + ux.bits.exp = iy + (BIAS + 512); + ux.e *= 0x1p-512; + } else { + ux.bits.exp = iy + BIAS; + } + x = ux.e * one; /* create necessary signal */ + return x; /* exact output */ +} diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h index c27a2be..97d1a58 100644 --- a/lib/msun/src/math.h +++ b/lib/msun/src/math.h @@ -423,9 +423,7 @@ long double floorl(long double); long double fmal(long double, long double, long double); long double fmaxl(long double, long double) __pure2; long double fminl(long double, long double) __pure2; -#if 0 long double fmodl(long double, long double); -#endif long double frexpl(long double value, int *); /* fundamentally !__pure2 */ long double hypotl(long double, long double); int ilogbl(long double) __pure2; |