toolchain/gcc/newlib/libm/mathfp/e_remainder.c

   1
   2 /* @(#)e_remainder.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, 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
  14 /*
  15 FUNCTION
  16 <<remainder>>, <<remainderf>>---round and  remainder
  17 INDEX
  18         remainder
  19 INDEX
  20         remainderf
  21
  22 ANSI_SYNOPSIS
  23         #include <math.h>
  24         double remainder(double <[x]>, double <[y]>);
  25         float remainderf(float <[x]>, float <[y]>);
  26
  27 TRAD_SYNOPSIS
  28         #include <math.h>
  29         double remainder(<[x]>,<[y]>)
  30         double <[x]>, <[y]>;
  31         float remainderf(<[x]>,<[y]>)
  32         float <[x]>, <[y]>;
  33
  34 DESCRIPTION
  35 <<remainder>> and <<remainderf>> find the remainder of
  36 <[x]>/<[y]>; this value is in the range -<[y]>/2 .. +<[y]>/2.
  37
  38 RETURNS
  39 <<remainder>> returns the integer result as a double.
  40
  41 PORTABILITY
  42 <<remainder>> is a System V release 4.
  43 <<remainderf>> is an extension.
  44
  45 */
  46
  47 /* remainder(x,p)
  48  * Return :
  49  *      returns  x REM p  =  x - [x/p]*p as if in infinite
  50  *      precise arithmetic, where [x/p] is the (infinite bit)
  51  *      integer nearest x/p (in half way case choose the even one).
  52  * Method :
  53  *      Based on fmod() return x-[x/p]chopped*p exactlp.
  54  */
  55
  56 #include "fdlibm.h"
  57
  58 #ifndef _DOUBLE_IS_32BITS
  59
  60 #ifdef __STDC__
  61 static const double zero = 0.0;
  62 #else
  63 static double zero = 0.0;
  64 #endif
  65
  66
  67 #ifdef __STDC__
  68         double remainder(double x, double p)
  69 #else
  70         double remainder(x,p)
  71         double x,p;
  72 #endif
  73 {
  74         __int32_t hx,hp;
  75         __uint32_t sx,lx,lp;
  76         double p_half;
  77
  78         EXTRACT_WORDS(hx,lx,x);
  79         EXTRACT_WORDS(hp,lp,p);
  80         sx = hx&0x80000000;
  81         hp &= 0x7fffffff;
  82         hx &= 0x7fffffff;
  83
  84     /* purge off exception values */
  85         if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
  86         if((hx>=0x7ff00000)||                   /* x not finite */
  87           ((hp>=0x7ff00000)&&                   /* p is NaN */
  88           (((hp-0x7ff00000)|lp)!=0)))
  89             return (x*p)/(x*p);
  90
  91
  92         if (hp<=0x7fdfffff) x = fmod(x,p+p);    /* now x < 2p */
  93         if (((hx-hp)|(lx-lp))==0) return zero*x;
  94         x  = fabs(x);
  95         p  = fabs(p);
  96         if (hp<0x00200000) {
  97             if(x+x>p) {
  98                 x-=p;
  99                 if(x+x>=p) x -= p;
 100             }
 101         } else {
 102             p_half = 0.5*p;
 103             if(x>p_half) {
 104                 x-=p;
 105                 if(x>=p_half) x -= p;
 106             }
 107         }
 108         GET_HIGH_WORD(hx,x);
 109         SET_HIGH_WORD(x,hx^sx);
 110         return x;
 111 }
 112
 113 #endif /* defined(_DOUBLE_IS_32BITS) */