rjw | 1f88458 | 2022-01-06 17:20:42 +0800 | [diff] [blame^] | 1 | /* e_rem_pio2f.c -- float version of e_rem_pio2.c |
| 2 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
| 3 | * Debugged and optimized by Bruce D. Evans. |
| 4 | */ |
| 5 | |
| 6 | /* |
| 7 | * ==================================================== |
| 8 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 9 | * |
| 10 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 11 | * Permission to use, copy, modify, and distribute this |
| 12 | * software is freely granted, provided that this notice |
| 13 | * is preserved. |
| 14 | * ==================================================== |
| 15 | */ |
| 16 | |
| 17 | #include <sys/cdefs.h> |
| 18 | __FBSDID("$FreeBSD$"); |
| 19 | |
| 20 | /* __ieee754_rem_pio2f(x,y) |
| 21 | * |
| 22 | * return the remainder of x rem pi/2 in *y |
| 23 | * use double precision for everything except passing x |
| 24 | * use __kernel_rem_pio2() for large x |
| 25 | */ |
| 26 | |
| 27 | #include <float.h> |
| 28 | |
| 29 | #include "math.h" |
| 30 | #include "math_private.h" |
| 31 | |
| 32 | /* |
| 33 | * invpio2: 53 bits of 2/pi |
| 34 | * pio2_1: first 25 bits of pi/2 |
| 35 | * pio2_1t: pi/2 - pio2_1 |
| 36 | */ |
| 37 | |
| 38 | static const double |
| 39 | invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ |
| 40 | pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ |
| 41 | pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ |
| 42 | |
| 43 | #ifdef INLINE_REM_PIO2F |
| 44 | static __inline __always_inline |
| 45 | #endif |
| 46 | int |
| 47 | __ieee754_rem_pio2f(float x, double *y) |
| 48 | { |
| 49 | double w,r,fn; |
| 50 | double tx[1],ty[1]; |
| 51 | float z; |
| 52 | int32_t e0,n,ix,hx; |
| 53 | |
| 54 | GET_FLOAT_WORD(hx,x); |
| 55 | ix = hx&0x7fffffff; |
| 56 | /* 33+53 bit pi is good enough for medium size */ |
| 57 | if (ix<0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ |
| 58 | /* Use a specialized rint() to get fn. Assume round-to-nearest. */ |
| 59 | STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); |
| 60 | fn = fn-0x1.8p52; |
| 61 | #ifdef HAVE_EFFICIENT_IRINT |
| 62 | n = irint(fn); |
| 63 | #else |
| 64 | n = (int32_t)fn; |
| 65 | #endif |
| 66 | r = x-fn*pio2_1; |
| 67 | w = fn*pio2_1t; |
| 68 | *y = r-w; |
| 69 | return n; |
| 70 | } |
| 71 | /* |
| 72 | * all other (large) arguments |
| 73 | */ |
| 74 | if (ix>=0x7f800000) { /* x is inf or NaN */ |
| 75 | *y=x-x; |
| 76 | return 0; |
| 77 | } |
| 78 | /* set z = scalbn(|x|,ilogb(|x|)-23) */ |
| 79 | e0 = (ix>>23)-150; /* e0 = ilogb(|x|)-23; */ |
| 80 | SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23))); |
| 81 | tx[0] = z; |
| 82 | n = __kernel_rem_pio2(tx,ty,e0,1,0); |
| 83 | if (hx<0) {*y = -ty[0]; return -n;} |
| 84 | *y = ty[0]; |
| 85 | return n; |
| 86 | } |