src/math/fmodl.c

   1 #include "libm.h"
   2
   3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4 long double fmodl(long double x, long double y)
   5 {
   6         return fmod(x, y);
   7 }
   8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
   9 long double fmodl(long double x, long double y)
  10 {
  11         union ldshape ux = {x}, uy = {y};
  12         int ex = ux.i.se & 0x7fff;
  13         int ey = uy.i.se & 0x7fff;
  14         int sx = ux.i.se & 0x8000;
  15
  16         if (y == 0 || isnan(y) || ex == 0x7fff)
  17                 return (x*y)/(x*y);
  18         ux.i.se = ex;
  19         uy.i.se = ey;
  20         if (ux.f <= uy.f) {
  21                 if (ux.f == uy.f)
  22                         return 0*x;
  23                 return x;
  24         }
  25
  26         /* normalize x and y */
  27         if (!ex) {
  28                 ux.f *= 0x1p120f;
  29                 ex = ux.i.se - 120;
  30         }
  31         if (!ey) {
  32                 uy.f *= 0x1p120f;
  33                 ey = uy.i.se - 120;
  34         }
  35
  36         /* x mod y */
  37 #if LDBL_MANT_DIG == 64
  38         uint64_t i, mx, my;
  39         mx = ux.i.m;
  40         my = uy.i.m;
  41         for (; ex > ey; ex--) {
  42                 i = mx - my;
  43                 if (mx >= my) {
  44                         if (i == 0)
  45                                 return 0*x;
  46                         mx = 2*i;
  47                 } else if (2*mx < mx) {
  48                         mx = 2*mx - my;
  49                 } else {
  50                         mx = 2*mx;
  51                 }
  52         }
  53         i = mx - my;
  54         if (mx >= my) {
  55                 if (i == 0)
  56                         return 0*x;
  57                 mx = i;
  58         }
  59         for (; mx >> 63 == 0; mx *= 2, ex--);
  60         ux.i.m = mx;
  61 #elif LDBL_MANT_DIG == 113
  62         uint64_t hi, lo, xhi, xlo, yhi, ylo;
  63         xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
  64         yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
  65         xlo = ux.i2.lo;
  66         ylo = ux.i2.lo;
  67         for (; ex > ey; ex--) {
  68                 hi = xhi - yhi;
  69                 lo = xlo - ylo;
  70                 if (xlo < ylo)
  71                         hi -= 1;
  72                 if (hi >> 63 == 0) {
  73                         if ((hi|lo) == 0)
  74                                 return 0*x;
  75                         xhi = 2*hi + (lo>>63);
  76                         xlo = 2*lo;
  77                 } else {
  78                         xhi = 2*xhi + (xlo>>63);
  79                         xlo = 2*xlo;
  80                 }
  81         }
  82         hi = xhi - yhi;
  83         lo = xlo - ylo;
  84         if (xlo < ylo)
  85                 hi -= 1;
  86         if (hi >> 63 == 0) {
  87                 if ((hi|lo) == 0)
  88                         return 0*x;
  89                 xhi = hi;
  90                 xlo = lo;
  91         }
  92         for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
  93         ux.i2.hi = xhi;
  94         ux.i2.lo = xlo;
  95 #endif
  96
  97         /* scale result */
  98         if (ex <= 0) {
  99                 ux.i.se = (ex+120)|sx;
 100                 ux.f *= 0x1p-120f;
 101         } else
 102                 ux.i.se = ex|sx;
 103         return ux.f;
 104 }
 105 #endif