src/math/fmod.c

   1 #include <math.h>
   2 #include <stdint.h>
   3
   4 double fmod(double x, double y)
   5 {
   6         union {double f; uint64_t i;} ux = {x}, uy = {y};
   7         int ex = ux.i>>52 & 0x7ff;
   8         int ey = uy.i>>52 & 0x7ff;
   9         int sx = ux.i>>63;
  10         uint64_t i;
  11
  12         /* in the followings uxi should be ux.i, but then gcc wrongly adds */
  13         /* float load/store to inner loops ruining performance and code size */
  14         uint64_t uxi = ux.i;
  15
  16         if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
  17                 return (x*y)/(x*y);
  18         if (uxi<<1 <= uy.i<<1) {
  19                 if (uxi<<1 == uy.i<<1)
  20                         return 0*x;
  21                 return x;
  22         }
  23
  24         /* normalize x and y */
  25         if (!ex) {
  26                 for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
  27                 uxi <<= -ex + 1;
  28         } else {
  29                 uxi &= -1ULL >> 12;
  30                 uxi |= 1ULL << 52;
  31         }
  32         if (!ey) {
  33                 for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
  34                 uy.i <<= -ey + 1;
  35         } else {
  36                 uy.i &= -1ULL >> 12;
  37                 uy.i |= 1ULL << 52;
  38         }
  39
  40         /* x mod y */
  41         for (; ex > ey; ex--) {
  42                 i = uxi - uy.i;
  43                 if (i >> 63 == 0) {
  44                         if (i == 0)
  45                                 return 0*x;
  46                         uxi = i;
  47                 }
  48                 uxi <<= 1;
  49         }
  50         i = uxi - uy.i;
  51         if (i >> 63 == 0) {
  52                 if (i == 0)
  53                         return 0*x;
  54                 uxi = i;
  55         }
  56         for (; uxi>>52 == 0; uxi <<= 1, ex--);
  57
  58         /* scale result */
  59         if (ex > 0) {
  60                 uxi -= 1ULL << 52;
  61                 uxi |= (uint64_t)ex << 52;
  62         } else {
  63                 uxi >>= -ex + 1;
  64         }
  65         uxi |= (uint64_t)sx << 63;
  66         ux.i = uxi;
  67         return ux.f;
  68 }