Home | History | Annotate | Line # | Download | only in ld128
e_rem_pio2l.h revision 1.2 
      1 /* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
      2 /*
      3  * ====================================================
      4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      5  * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
      6  *
      7  * Developed at SunSoft, a Sun Microsystems, Inc. business.
      8  * Permission to use, copy, modify, and distribute this
      9  * software is freely granted, provided that this notice
     10  * is preserved.
     11  * ====================================================
     12  *
     13  * Optimized by Bruce D. Evans.
     14  */
     15 
     16 #include <sys/cdefs.h>
     17 #if 0
     18 __FBSDID("$FreeBSD: head/lib/msun/ld128/e_rem_pio2l.h 336545 2018-07-20 12:42:24Z bde $");
     19 #endif
     20 
     21 /* ld128 version of __ieee754_rem_pio2l(x,y)
     22  *
     23  * return the remainder of x rem pi/2 in y[0]+y[1]
     24  * use __kernel_rem_pio2()
     25  */
     26 
     27 #include <float.h>
     28 #include <machine/ieee.h>
     29 
     30 #include "math.h"
     31 #include "math_private.h"
     32 
     33 #define	BIAS	(LDBL_MAX_EXP - 1)
     34 
     35 /*
     36  * XXX need to verify that nonzero integer multiples of pi/2 within the
     37  * range get no closer to a long double than 2**-140, or that
     38  * ilogb(x) + ilogb(min_delta) < 45 - -140.
     39  */
     40 /*
     41  * invpio2:  113 bits of 2/pi
     42  * pio2_1:   first  68 bits of pi/2
     43  * pio2_1t:  pi/2 - pio2_1
     44  * pio2_2:   second 68 bits of pi/2
     45  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
     46  * pio2_3:   third  68 bits of pi/2
     47  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
     48  */
     49 
     50 static const double
     51 zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
     52 two24 =  1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */
     53 
     54 static const long double
     55 invpio2 =  6.3661977236758134307553505349005747e-01L,	/*  0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
     56 pio2_1  =  1.5707963267948966192292994253909555e+00L,	/*  0x1921fb54442d18469800000000000.0p-112 */
     57 pio2_1t =  2.0222662487959507323996846200947577e-21L,	/*  0x13198a2e03707344a4093822299f3.0p-181 */
     58 pio2_2  =  2.0222662487959507323994779168837751e-21L,	/*  0x13198a2e03707344a400000000000.0p-181 */
     59 pio2_2t =  2.0670321098263988236496903051604844e-43L,	/*  0x127044533e63a0105df531d89cd91.0p-254 */
     60 pio2_3  =  2.0670321098263988236499468110329591e-43L,	/*  0x127044533e63a0105e00000000000.0p-254 */
     61 pio2_3t = -2.5650587247459238361625433492959285e-65L;	/* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
     62 
     63 static inline __always_inline int
     64 __ieee754_rem_pio2l(long double x, long double *y)
     65 {
     66 	union ieee_ext_u u,u1;
     67 	long double z,w,t,r,fn;
     68 	double tx[5],ty[3];
     69 	int64_t n;
     70 	int e0,ex,i,j,nx;
     71 	int16_t expsign, expsign1;
     72 
     73 	u.extu_ld = x;
     74 	ex = u.extu_exp;
     75 	expsign = u.extu_exp | (u.extu_sign << 15);
     76 	if (ex < BIAS + 45 || (ex == BIAS + 45 &&
     77 	    u.extu_frach < 0x921fb54442d1LL)) {
     78 	    /* |x| ~< 2^45*(pi/2), medium size */
     79 	    /* TODO: use only double precision for fn, as in expl(). */
     80 	    fn = rnintl(x * invpio2);
     81 	    n  = i64rint(fn);
     82 	    r  = x-fn*pio2_1;
     83 	    w  = fn*pio2_1t;	/* 1st round good to 180 bit */
     84 	    {
     85 		union ieee_ext_u u2;
     86 	        int ex1;
     87 	        j  = ex;
     88 	        y[0] = r-w;
     89 		u2.extu_ld = y[0];
     90 		ex1 = u2.extu_exp;
     91 	        i = j-ex1;
     92 	        if(i>51) {  /* 2nd iteration needed, good to 248 */
     93 		    t  = r;
     94 		    w  = fn*pio2_2;
     95 		    r  = t-w;
     96 		    w  = fn*pio2_2t-((t-r)-w);
     97 		    y[0] = r-w;
     98 		    u2.extu_ld = y[0];
     99 		    ex1 = u2.extu_exp;
    100 		    i = j-ex1;
    101 		    if(i>119) {	/* 3rd iteration need, 316 bits acc */
    102 			t  = r;	/* will cover all possible cases */
    103 			w  = fn*pio2_3;
    104 			r  = t-w;
    105 			w  = fn*pio2_3t-((t-r)-w);
    106 			y[0] = r-w;
    107 		    }
    108 		}
    109 	    }
    110 	    y[1] = (r-y[0])-w;
    111 	    return n;
    112 	}
    113     /*
    114      * all other (large) arguments
    115      */
    116 	if(ex==0x7fff) {		/* x is inf or NaN */
    117 	    y[0]=y[1]=x-x; return 0;
    118 	}
    119     /* set z = scalbn(|x|,ilogb(x)-23) */
    120 	u1.extu_ld = x;
    121 	e0 = ex - BIAS - 23;		/* e0 = ilogb(|x|)-23; */
    122 	expsign1 = ex - e0;
    123 	u1.extu_exp = expsign1;
    124 	u1.extu_sign = expsign1 >> 15;
    125 	z = u1.extu_ld;
    126 	for(i=0;i<4;i++) {
    127 		tx[i] = (double)((int32_t)(z));
    128 		z     = (z-tx[i])*two24;
    129 	}
    130 	tx[4] = z;
    131 	nx = 5;
    132 	while(tx[nx-1]==zero) nx--;	/* skip zero term */
    133 	n  =  __kernel_rem_pio2(tx,ty,e0,nx,3);
    134 	t = (long double)ty[2] + ty[1];
    135 	r = t + ty[0];
    136 	w = ty[0] - (r - t);
    137 	if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
    138 	y[0] = r; y[1] = w; return n;
    139 }
    140