Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /* ef_fmod.c -- float version of e_fmod.c.
  2.  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  3.  */
  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. /*
  17.  * __ieee754_fmodf(x,y)
  18.  * Return x mod y in exact arithmetic
  19.  * Method: shift and subtract
  20.  */
  21.  
  22. #include "fdlibm.h"
  23.  
  24. #ifdef __STDC__
  25. static const float one = 1.0, Zero[] = {0.0, -0.0,};
  26. #else
  27. static float one = 1.0, Zero[] = {0.0, -0.0,};
  28. #endif
  29.  
  30. #ifdef __STDC__
  31.         float __ieee754_fmodf(float x, float y)
  32. #else
  33.         float __ieee754_fmodf(x,y)
  34.         float x,y ;
  35. #endif
  36. {
  37.         __int32_t n,hx,hy,hz,ix,iy,sx,i;
  38.  
  39.         GET_FLOAT_WORD(hx,x);
  40.         GET_FLOAT_WORD(hy,y);
  41.         sx = hx&0x80000000;             /* sign of x */
  42.         hx ^=sx;                /* |x| */
  43.         hy &= 0x7fffffff;       /* |y| */
  44.  
  45.     /* purge off exception values */
  46.         if(FLT_UWORD_IS_ZERO(hy)||
  47.            !FLT_UWORD_IS_FINITE(hx)||
  48.            FLT_UWORD_IS_NAN(hy))
  49.             return (x*y)/(x*y);
  50.         if(hx<hy) return x;                     /* |x|<|y| return x */
  51.         if(hx==hy)
  52.             return Zero[(__uint32_t)sx>>31];    /* |x|=|y| return x*0*/
  53.  
  54.     /* Note: y cannot be zero if we reach here. */
  55.  
  56.     /* determine ix = ilogb(x) */
  57.         if(FLT_UWORD_IS_SUBNORMAL(hx)) {        /* subnormal x */
  58.             for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
  59.         } else ix = (hx>>23)-127;
  60.  
  61.     /* determine iy = ilogb(y) */
  62.         if(FLT_UWORD_IS_SUBNORMAL(hy)) {        /* subnormal y */
  63.             for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
  64.         } else iy = (hy>>23)-127;
  65.  
  66.     /* set up {hx,lx}, {hy,ly} and align y to x */
  67.         if(ix >= -126)
  68.             hx = 0x00800000|(0x007fffff&hx);
  69.         else {          /* subnormal x, shift x to normal */
  70.             n = -126-ix;
  71.             hx = hx<<n;
  72.         }
  73.         if(iy >= -126)
  74.             hy = 0x00800000|(0x007fffff&hy);
  75.         else {          /* subnormal y, shift y to normal */
  76.             n = -126-iy;
  77.             hy = hy<<n;
  78.         }
  79.  
  80.     /* fix point fmod */
  81.         n = ix - iy;
  82.         while(n--) {
  83.             hz=hx-hy;
  84.             if(hz<0){hx = hx+hx;}
  85.             else {
  86.                 if(hz==0)               /* return sign(x)*0 */
  87.                     return Zero[(__uint32_t)sx>>31];
  88.                 hx = hz+hz;
  89.             }
  90.         }
  91.         hz=hx-hy;
  92.         if(hz>=0) {hx=hz;}
  93.  
  94.     /* convert back to floating value and restore the sign */
  95.         if(hx==0)                       /* return sign(x)*0 */
  96.             return Zero[(__uint32_t)sx>>31];   
  97.         while(hx<0x00800000) {          /* normalize x */
  98.             hx = hx+hx;
  99.             iy -= 1;
  100.         }
  101.         if(iy>= -126) {         /* normalize output */
  102.             hx = ((hx-0x00800000)|((iy+127)<<23));
  103.             SET_FLOAT_WORD(x,hx|sx);
  104.         } else {                /* subnormal output */
  105.             /* If denormals are not supported, this code will generate a
  106.                zero representation.  */
  107.             n = -126 - iy;
  108.             hx >>= n;
  109.             SET_FLOAT_WORD(x,hx|sx);
  110.             x *= one;           /* create necessary signal */
  111.         }
  112.         return x;               /* exact output */
  113. }
  114.