Branch data     Line data    Source code

       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 "cdefs-compat.h"
      18                 :             : //__FBSDID("$FreeBSD: src/lib/msun/src/e_rem_pio2f.c,v 1.32 2009/06/03 08:16:34 ed Exp $");
      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                 :             : #include <openlibm_math.h>
      29                 :             : 
      30                 :             : #include "math_private.h"
      31                 :             : 
      32                 :             : /*
      33                 :             :  * invpio2:  53 bits of 2/pi
      34                 :             :  * pio2_1:   first  33 bit 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                 :             : __inline int
      44                 :          72 : __ieee754_rem_pio2f(float x, double *y)
      45                 :             : {
      46                 :             :         double w,r,fn;
      47                 :             :         double tx[1],ty[1];
      48                 :             :         float z;
      49                 :             :         int32_t e0,n,ix,hx;
      50                 :             : 
      51                 :          72 :         GET_FLOAT_WORD(hx,x);
      52                 :          72 :         ix = hx&0x7fffffff;
      53                 :             :     /* 33+53 bit pi is good enough for medium size */
      54         [ +  - ]:          72 :         if(ix<0x4dc90fdb) {          /* |x| ~< 2^28*(pi/2), medium size */
      55                 :             :             /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
      56                 :          72 :             STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52);
      57                 :          72 :             fn = fn-0x1.8p52;
      58                 :             : #ifdef HAVE_EFFICIENT_IRINT
      59                 :             :             n  = irint(fn);
      60                 :             : #else
      61                 :          72 :             n  = (int32_t)fn;
      62                 :             : #endif
      63                 :          72 :             r  = x-fn*pio2_1;
      64                 :          72 :             w  = fn*pio2_1t;
      65                 :          72 :             *y = r-w;
      66                 :          72 :             return n;
      67                 :             :         }
      68                 :             :     /*
      69                 :             :      * all other (large) arguments
      70                 :             :      */
      71         [ #  # ]:           0 :         if(ix>=0x7f800000) {         /* x is inf or NaN */
      72                 :           0 :             *y=x-x; return 0;
      73                 :             :         }
      74                 :             :     /* set z = scalbn(|x|,ilogb(|x|)-23) */
      75                 :           0 :         e0 = (ix>>23)-150;                /* e0 = ilogb(|x|)-23; */
      76                 :           0 :         SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
      77                 :           0 :         tx[0] = z;
      78                 :           0 :         n  =  __kernel_rem_pio2(tx,ty,e0,1,0);
      79         [ #  # ]:           0 :         if(hx<0) {*y = -ty[0]; return -n;}
      80                 :           0 :         *y = ty[0]; return n;
      81                 :             : }