nsz Git - musl/blob - src/math/remainderf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_remainderf.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4  */
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #include "libm.h"
  17
  18 static const float zero = 0.0;
  19
  20 float remainderf(float x, float p)
  21 {
  22         int32_t hx,hp;
  23         uint32_t sx;
  24         float p_half;
  25
  26         GET_FLOAT_WORD(hx, x);
  27         GET_FLOAT_WORD(hp, p);
  28         sx = hx & 0x80000000;
  29         hp &= 0x7fffffff;
  30         hx &= 0x7fffffff;
  31
  32         /* purge off exception values */
  33         if (hp == 0)  /* p = 0 */
  34                 return (x*p)/(x*p);
  35         if (hx >= 0x7f800000 || hp > 0x7f800000)  /* x not finite, p is NaN */
  36                 // FIXME: why long double?
  37                 return ((long double)x*p)/((long double)x*p);
  38
  39         if (hp <= 0x7effffff)
  40                 x = fmodf(x, p + p);  /* now x < 2p */
  41         if (hx - hp == 0)
  42                 return zero*x;
  43         x = fabsf(x);
  44         p = fabsf(p);
  45         if (hp < 0x01000000) {
  46                 if (x + x > p) {
  47                         x -= p;
  48                         if (x + x >= p)
  49                                 x -= p;
  50                 }
  51         } else {
  52                 p_half = 0.5f*p;
  53                 if (x > p_half) {
  54                         x -= p;
  55                         if (x >= p_half)
  56                                 x -= p;
  57                 }
  58         }
  59         GET_FLOAT_WORD(hx, x);
  60         if ((hx & 0x7fffffff) == 0)
  61                 hx = 0;
  62         SET_FLOAT_WORD(x, hx ^ sx);
  63         return x;
  64 }