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

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_remquof.c */
   2 /*-
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /*
  13  * Return the IEEE remainder and set *quo to the last n bits of the
  14  * quotient, rounded to the nearest integer.  We choose n=31 because
  15  * we wind up computing all the integer bits of the quotient anyway as
  16  * a side-effect of computing the remainder by the shift and subtract
  17  * method.  In practice, this is far more bits than are needed to use
  18  * remquo in reduction algorithms.
  19  */
  20
  21 #include "libm.h"
  22
  23 static const float Zero[] = {0.0, -0.0,};
  24
  25 float remquof(float x, float y, int *quo)
  26 {
  27         int32_t n,hx,hy,hz,ix,iy,sx,i;
  28         uint32_t q,sxy;
  29
  30         GET_FLOAT_WORD(hx, x);
  31         GET_FLOAT_WORD(hy, y);
  32         sxy = (hx ^ hy) & 0x80000000;
  33         sx = hx & 0x80000000;   /* sign of x */
  34         hx ^= sx;               /* |x| */
  35         hy &= 0x7fffffff;       /* |y| */
  36
  37         /* purge off exception values */
  38         if (hy == 0 || hx >= 0x7f800000 || hy > 0x7f800000) /* y=0,NaN;or x not finite */
  39                 return (x*y)/(x*y);
  40         if (hx < hy) {       /* |x| < |y| return x or x-y */
  41                 q = 0;
  42                 goto fixup;
  43         } else if(hx==hy) {  /* |x| = |y| return x*0*/
  44                 *quo = 1;
  45                 return Zero[(uint32_t)sx>>31];
  46         }
  47
  48         /* determine ix = ilogb(x) */
  49         if (hx < 0x00800000) {  /* subnormal x */
  50                 for (ix = -126, i=hx<<8; i>0; i<<=1) ix--;
  51         } else
  52                 ix = (hx>>23) - 127;
  53
  54         /* determine iy = ilogb(y) */
  55         if (hy < 0x00800000) {  /* subnormal y */
  56                 for (iy = -126, i=hy<<8; i>0; i<<=1) iy--;
  57         } else
  58                 iy = (hy>>23) - 127;
  59
  60         /* set up {hx,lx}, {hy,ly} and align y to x */
  61         if (ix >= -126)
  62                 hx = 0x00800000|(0x007fffff&hx);
  63         else {  /* subnormal x, shift x to normal */
  64                 n = -126 - ix;
  65                 hx <<= n;
  66         }
  67         if (iy >= -126)
  68                 hy = 0x00800000|(0x007fffff&hy);
  69         else {  /* subnormal y, shift y to normal */
  70                 n = -126 - iy;
  71                 hy <<= n;
  72         }
  73
  74         /* fix point fmod */
  75         n = ix - iy;
  76         q = 0;
  77         while (n--) {
  78                 hz = hx - hy;
  79                 if (hz < 0)
  80                         hx = hx << 1;
  81                 else {
  82                         hx = hz << 1;
  83                         q++;
  84                 }
  85                 q <<= 1;
  86         }
  87         hz = hx - hy;
  88         if (hz >= 0) {
  89                 hx = hz;
  90                 q++;
  91         }
  92
  93         /* convert back to floating value and restore the sign */
  94         if (hx == 0) {                             /* return sign(x)*0 */
  95                 *quo = sxy ? -q : q;
  96                 return Zero[(uint32_t)sx>>31];
  97         }
  98         while (hx < 0x00800000) {  /* normalize x */
  99                 hx <<= 1;
 100                 iy--;
 101         }
 102         if (iy >= -126) {          /* normalize output */
 103                 hx = (hx-0x00800000)|((iy+127)<<23);
 104         } else {                   /* subnormal output */
 105                 n = -126 - iy;
 106                 hx >>= n;
 107         }
 108 fixup:
 109         SET_FLOAT_WORD(x,hx);
 110         y = fabsf(y);
 111         if (y < 0x1p-125f) {
 112                 if (x + x > y || (x + x == y && (q & 1))) {
 113                         q++;
 114                         x -= y;
 115                 }
 116         } else if (x > 0.5f*y || (x == 0.5f*y && (q & 1))) {
 117                 q++;
 118                 x -= y;
 119         }
 120         GET_FLOAT_WORD(hx, x);
 121         SET_FLOAT_WORD(x, hx ^ sx);
 122         q &= 0x7fffffff;
 123         *quo = sxy ? -q : q;
 124         return x;
 125 }