packages/language/c/libm/v2_0/src/double/ieee754-core/e_fmod.c

   1 //===========================================================================
   2 //
   3 //      e_fmod.c
   4 //
   5 //      Part of the standard mathematical function library
   6 //
   7 //===========================================================================
   8 //####ECOSGPLCOPYRIGHTBEGIN####
   9 // -------------------------------------------
  10 // This file is part of eCos, the Embedded Configurable Operating System.
  11 // Copyright (C) 1998, 1999, 2000, 2001, 2002 Red Hat, Inc.
  12 //
  13 // eCos is free software; you can redistribute it and/or modify it under
  14 // the terms of the GNU General Public License as published by the Free
  15 // Software Foundation; either version 2 or (at your option) any later version.
  16 //
  17 // eCos is distributed in the hope that it will be useful, but WITHOUT ANY
  18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or
  19 // FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  20 // for more details.
  21 //
  22 // You should have received a copy of the GNU General Public License along
  23 // with eCos; if not, write to the Free Software Foundation, Inc.,
  24 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
  25 //
  26 // As a special exception, if other files instantiate templates or use macros
  27 // or inline functions from this file, or you compile this file and link it
  28 // with other works to produce a work based on this file, this file does not
  29 // by itself cause the resulting work to be covered by the GNU General Public
  30 // License. However the source code for this file must still be made available
  31 // in accordance with section (3) of the GNU General Public License.
  32 //
  33 // This exception does not invalidate any other reasons why a work based on
  34 // this file might be covered by the GNU General Public License.
  35 //
  36 // Alternative licenses for eCos may be arranged by contacting Red Hat, Inc.
  37 // at http://sources.redhat.com/ecos/ecos-license/
  38 // -------------------------------------------
  39 //####ECOSGPLCOPYRIGHTEND####
  40 //===========================================================================
  41 //#####DESCRIPTIONBEGIN####
  42 //
  43 // Author(s):   jlarmour
  44 // Contributors:  jlarmour
  45 // Date:        1998-02-13
  46 // Purpose:
  47 // Description:
  48 // Usage:
  49 //
  50 //####DESCRIPTIONEND####
  51 //
  52 //===========================================================================
  53
  54 // CONFIGURATION
  55
  56 #include <pkgconf/libm.h>   // Configuration header
  57
  58 // Include the Math library?
  59 #ifdef CYGPKG_LIBM
  60
  61 // Derived from code with the following copyright
  62
  63
  64 /* @(#)e_fmod.c 1.3 95/01/18 */
  65 /*
  66  * ====================================================
  67  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  68  *
  69  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  70  * Permission to use, copy, modify, and distribute this
  71  * software is freely granted, provided that this notice
  72  * is preserved.
  73  * ====================================================
  74  */
  75
  76 /*
  77  * __ieee754_fmod(x,y)
  78  * Return x mod y in exact arithmetic
  79  * Method: shift and subtract
  80  */
  81
  82 #include "mathincl/fdlibm.h"
  83
  84 static const double one = 1.0, Zero[] = {0.0, -0.0,};
  85
  86         double __ieee754_fmod(double x, double y)
  87 {
  88         int n,hx,hy,hz,ix,iy,sx,i;
  89         unsigned lx,ly,lz;
  90
  91         hx = CYG_LIBM_HI(x);            /* high word of x */
  92         lx = CYG_LIBM_LO(x);            /* low  word of x */
  93         hy = CYG_LIBM_HI(y);            /* high word of y */
  94         ly = CYG_LIBM_LO(y);            /* low  word of y */
  95         sx = hx&0x80000000;             /* sign of x */
  96         hx ^=sx;                /* |x| */
  97         hy &= 0x7fffffff;       /* |y| */
  98
  99     /* purge off exception values */
 100         if((hy|ly)==0||(hx>=0x7ff00000)||       /* y=0,or x not finite */
 101           ((hy|((ly|-ly)>>31))>0x7ff00000))     /* or y is NaN */
 102             return (x*y)/(x*y);
 103         if(hx<=hy) {
 104             if((hx<hy)||(lx<ly)) return x;      /* |x|<|y| return x */
 105             if(lx==ly)
 106                 return Zero[(unsigned)sx>>31];  /* |x|=|y| return x*0*/
 107         }
 108
 109     /* determine ix = ilogb(x) */
 110         if(hx<0x00100000) {     /* subnormal x */
 111             if(hx==0) {
 112                 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
 113             } else {
 114                 for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
 115             }
 116         } else ix = (hx>>20)-1023;
 117
 118     /* determine iy = ilogb(y) */
 119         if(hy<0x00100000) {     /* subnormal y */
 120             if(hy==0) {
 121                 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
 122             } else {
 123                 for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
 124             }
 125         } else iy = (hy>>20)-1023;
 126
 127     /* set up {hx,lx}, {hy,ly} and align y to x */
 128         if(ix >= -1022)
 129             hx = 0x00100000|(0x000fffff&hx);
 130         else {          /* subnormal x, shift x to normal */
 131             n = -1022-ix;
 132             if(n<=31) {
 133                 hx = (hx<<n)|(lx>>(32-n));
 134                 lx <<= n;
 135             } else {
 136                 hx = lx<<(n-32);
 137                 lx = 0;
 138             }
 139         }
 140         if(iy >= -1022)
 141             hy = 0x00100000|(0x000fffff&hy);
 142         else {          /* subnormal y, shift y to normal */
 143             n = -1022-iy;
 144             if(n<=31) {
 145                 hy = (hy<<n)|(ly>>(32-n));
 146                 ly <<= n;
 147             } else {
 148                 hy = ly<<(n-32);
 149                 ly = 0;
 150             }
 151         }
 152
 153     /* fix point fmod */
 154         n = ix - iy;
 155         while(n--) {
 156             hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 157             if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
 158             else {
 159                 if((hz|lz)==0)          /* return sign(x)*0 */
 160                     return Zero[(unsigned)sx>>31];
 161                 hx = hz+hz+(lz>>31); lx = lz+lz;
 162             }
 163         }
 164         hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 165         if(hz>=0) {hx=hz;lx=lz;}
 166
 167     /* convert back to floating value and restore the sign */
 168         if((hx|lx)==0)                  /* return sign(x)*0 */
 169             return Zero[(unsigned)sx>>31];
 170         while(hx<0x00100000) {          /* normalize x */
 171             hx = hx+hx+(lx>>31); lx = lx+lx;
 172             iy -= 1;
 173         }
 174         if(iy>= -1022) {        /* normalize output */
 175             hx = ((hx-0x00100000)|((iy+1023)<<20));
 176             CYG_LIBM_HI(x) = hx|sx;
 177             CYG_LIBM_LO(x) = lx;
 178         } else {                /* subnormal output */
 179             n = -1022 - iy;
 180             if(n<=20) {
 181                 lx = (lx>>n)|((unsigned)hx<<(32-n));
 182                 hx >>= n;
 183             } else if (n<=31) {
 184                 lx = (hx<<(32-n))|(lx>>n); hx = sx;
 185             } else {
 186                 lx = hx>>(n-32); hx = sx;
 187             }
 188             CYG_LIBM_HI(x) = hx|sx;
 189             CYG_LIBM_LO(x) = lx;
 190             x *= one;           /* create necessary signal */
 191         }
 192         return x;               /* exact output */
 193 }
 194
 195 #endif // ifdef CYGPKG_LIBM
 196
 197 // EOF e_fmod.c