sysdeps/ieee754/ldbl-96/e_remainderl.c

   1 /* e_remainderl.c -- long double version of e_remainder.c.
   2  * Conversion to long double by Ulrich Drepper,
   3  * Cygnus Support, drepper@cygnus.com.
   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 /* __ieee754_remainderl(x,p)
  18  * Return :
  19  *      returns  x REM p  =  x - [x/p]*p as if in infinite
  20  *      precise arithmetic, where [x/p] is the (infinite bit)
  21  *      integer nearest x/p (in half way case choose the even one).
  22  * Method :
  23  *      Based on fmod() return x-[x/p]chopped*p exactlp.
  24  */
  25
  26 #include "math.h"
  27 #include "math_private.h"
  28
  29 static const long double zero = 0.0;
  30
  31
  32 long double
  33 __ieee754_remainderl(long double x, long double p)
  34 {
  35         u_int32_t sx,sex,sep,x0,x1,p0,p1;
  36         long double p_half;
  37
  38         GET_LDOUBLE_WORDS(sex,x0,x1,x);
  39         GET_LDOUBLE_WORDS(sep,p0,p1,p);
  40         sx = sex&0x8000;
  41         sep &= 0x7fff;
  42         sex &= 0x7fff;
  43
  44     /* purge off exception values */
  45         if((sep|p0|p1)==0) return (x*p)/(x*p);  /* p = 0 */
  46         if((sex==0x7fff)||                      /* x not finite */
  47           ((sep==0x7fff)&&                      /* p is NaN */
  48            ((p0|p1)!=0)))
  49             return (x*p)/(x*p);
  50
  51
  52         if (sep<0x7ffe) x = __ieee754_fmodl(x,p+p);     /* now x < 2p */
  53         if (((sex-sep)|(x0-p0)|(x1-p1))==0) return zero*x;
  54         x  = fabsl(x);
  55         p  = fabsl(p);
  56         if (sep<0x0002) {
  57             if(x+x>p) {
  58                 x-=p;
  59                 if(x+x>=p) x -= p;
  60             }
  61         } else {
  62             p_half = 0.5*p;
  63             if(x>p_half) {
  64                 x-=p;
  65                 if(x>=p_half) x -= p;
  66             }
  67         }
  68         GET_LDOUBLE_EXP(sex,x);
  69         SET_LDOUBLE_EXP(x,sex^sx);
  70         return x;
  71 }
  72 strong_alias (__ieee754_remainderl, __remainderl_finite)