sysdeps/libm-ieee754/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 #if defined(LIBM_SCCS) && !defined(lint)
  18 static char rcsid[] = "$NetBSD: $";
  19 #endif
  20
  21 /* __ieee754_remainderl(x,p)
  22  * Return :
  23  *      returns  x REM p  =  x - [x/p]*p as if in infinite
  24  *      precise arithmetic, where [x/p] is the (infinite bit)
  25  *      integer nearest x/p (in half way case choose the even one).
  26  * Method :
  27  *      Based on fmod() return x-[x/p]chopped*p exactlp.
  28  */
  29
  30 #include "math.h"
  31 #include "math_private.h"
  32
  33 #ifdef __STDC__
  34 static const long double zero = 0.0;
  35 #else
  36 static long double zero = 0.0;
  37 #endif
  38
  39
  40 #ifdef __STDC__
  41         long double __ieee754_remainderl(long double x, long double p)
  42 #else
  43         long double __ieee754_remainderl(x,p)
  44         long double x,p;
  45 #endif
  46 {
  47         u_int32_t sx,sex,sep,x0,x1,p0,p1;
  48         long double p_half;
  49
  50         GET_LDOUBLE_WORDS(sex,x0,x1,x);
  51         GET_LDOUBLE_WORDS(sep,p0,p1,p);
  52         sx = sex&0x8000;
  53         sep &= 0x7fff;
  54         sex &= 0x7fff;
  55
  56     /* purge off exception values */
  57         if((sep|p0|p1)==0) return (x*p)/(x*p);  /* p = 0 */
  58         if((sex==0x7fff)||                      /* x not finite */
  59           ((sep==0x7fff)&&                      /* p is NaN */
  60            ((p0|p1)!=0)))
  61             return (x*p)/(x*p);
  62
  63
  64         if (sep<0x7ffe) x = __ieee754_fmodl(x,p+p);     /* now x < 2p */
  65         if (((sex-sep)|(x0-p0)|(x1-p1))==0) return zero*x;
  66         x  = fabsl(x);
  67         p  = fabsl(p);
  68         if (sep<0x0002) {
  69             if(x+x>p) {
  70                 x-=p;
  71                 if(x+x>=p) x -= p;
  72             }
  73         } else {
  74             p_half = 0.5*p;
  75             if(x>p_half) {
  76                 x-=p;
  77                 if(x>=p_half) x -= p;
  78             }
  79         }
  80         GET_LDOUBLE_EXP(sex,x);
  81         SET_LDOUBLE_EXP(x,sex^sx);
  82         return x;
  83 }