Back to index

lightning-sunbird  0.9+nobinonly
e_fmod.c
Go to the documentation of this file.
00001 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002  *
00003  * ***** BEGIN LICENSE BLOCK *****
00004  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00005  *
00006  * The contents of this file are subject to the Mozilla Public License Version
00007  * 1.1 (the "License"); you may not use this file except in compliance with
00008  * the License. You may obtain a copy of the License at
00009  * http://www.mozilla.org/MPL/
00010  *
00011  * Software distributed under the License is distributed on an "AS IS" basis,
00012  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00013  * for the specific language governing rights and limitations under the
00014  * License.
00015  *
00016  * The Original Code is Mozilla Communicator client code, released
00017  * March 31, 1998.
00018  *
00019  * The Initial Developer of the Original Code is
00020  * Sun Microsystems, Inc.
00021  * Portions created by the Initial Developer are Copyright (C) 1998
00022  * the Initial Developer. All Rights Reserved.
00023  *
00024  * Contributor(s):
00025  *
00026  * Alternatively, the contents of this file may be used under the terms of
00027  * either of the GNU General Public License Version 2 or later (the "GPL"),
00028  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029  * in which case the provisions of the GPL or the LGPL are applicable instead
00030  * of those above. If you wish to allow use of your version of this file only
00031  * under the terms of either the GPL or the LGPL, and not to allow others to
00032  * use your version of this file under the terms of the MPL, indicate your
00033  * decision by deleting the provisions above and replace them with the notice
00034  * and other provisions required by the GPL or the LGPL. If you do not delete
00035  * the provisions above, a recipient may use your version of this file under
00036  * the terms of any one of the MPL, the GPL or the LGPL.
00037  *
00038  * ***** END LICENSE BLOCK ***** */
00039 
00040 /* @(#)e_fmod.c 1.3 95/01/18 */
00041 /*
00042  * ====================================================
00043  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00044  *
00045  * Developed at SunSoft, a Sun Microsystems, Inc. business.
00046  * Permission to use, copy, modify, and distribute this
00047  * software is freely granted, provided that this notice 
00048  * is preserved.
00049  * ====================================================
00050  */
00051 
00052 /* 
00053  * __ieee754_fmod(x,y)
00054  * Return x mod y in exact arithmetic
00055  * Method: shift and subtract
00056  */
00057 
00058 #include "fdlibm.h"
00059 
00060 #ifdef __STDC__
00061 static const double one = 1.0, Zero[] = {0.0, -0.0,};
00062 #else
00063 static double one = 1.0, Zero[] = {0.0, -0.0,};
00064 #endif
00065 
00066 #ifdef __STDC__
00067        double __ieee754_fmod(double x, double y)
00068 #else
00069        double __ieee754_fmod(x,y)
00070        double x,y ;
00071 #endif
00072 {
00073         fd_twoints ux, uy;
00074        int n,hx,hy,hz,ix,iy,sx,i;
00075        unsigned lx,ly,lz;
00076 
00077         ux.d = x; uy.d = y;
00078        hx = __HI(ux);              /* high word of x */
00079        lx = __LO(ux);              /* low  word of x */
00080        hy = __HI(uy);              /* high word of y */
00081        ly = __LO(uy);              /* low  word of y */
00082        sx = hx&0x80000000;         /* sign of x */
00083        hx ^=sx;             /* |x| */
00084        hy &= 0x7fffffff;    /* |y| */
00085 
00086     /* purge off exception values */
00087        if((hy|ly)==0||(hx>=0x7ff00000)||  /* y=0,or x not finite */
00088          ((hy|((ly|-(int)ly)>>31))>0x7ff00000))  /* or y is NaN */
00089            return (x*y)/(x*y);
00090        if(hx<=hy) {
00091            if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
00092            if(lx==ly) 
00093               return Zero[(unsigned)sx>>31];     /* |x|=|y| return x*0*/
00094        }
00095 
00096     /* determine ix = ilogb(x) */
00097        if(hx<0x00100000) {  /* subnormal x */
00098            if(hx==0) {
00099               for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
00100            } else {
00101               for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
00102            }
00103        } else ix = (hx>>20)-1023;
00104 
00105     /* determine iy = ilogb(y) */
00106        if(hy<0x00100000) {  /* subnormal y */
00107            if(hy==0) {
00108               for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
00109            } else {
00110               for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
00111            }
00112        } else iy = (hy>>20)-1023;
00113 
00114     /* set up {hx,lx}, {hy,ly} and align y to x */
00115        if(ix >= -1022) 
00116            hx = 0x00100000|(0x000fffff&hx);
00117        else {        /* subnormal x, shift x to normal */
00118            n = -1022-ix;
00119            if(n<=31) {
00120                hx = (hx<<n)|(lx>>(32-n));
00121                lx <<= n;
00122            } else {
00123               hx = lx<<(n-32);
00124               lx = 0;
00125            }
00126        }
00127        if(iy >= -1022) 
00128            hy = 0x00100000|(0x000fffff&hy);
00129        else {        /* subnormal y, shift y to normal */
00130            n = -1022-iy;
00131            if(n<=31) {
00132                hy = (hy<<n)|(ly>>(32-n));
00133                ly <<= n;
00134            } else {
00135               hy = ly<<(n-32);
00136               ly = 0;
00137            }
00138        }
00139 
00140     /* fix point fmod */
00141        n = ix - iy;
00142        while(n--) {
00143            hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00144            if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
00145            else {
00146               if((hz|lz)==0)              /* return sign(x)*0 */
00147                   return Zero[(unsigned)sx>>31];
00148               hx = hz+hz+(lz>>31); lx = lz+lz;
00149            }
00150        }
00151        hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00152        if(hz>=0) {hx=hz;lx=lz;}
00153 
00154     /* convert back to floating value and restore the sign */
00155        if((hx|lx)==0)                     /* return sign(x)*0 */
00156            return Zero[(unsigned)sx>>31]; 
00157        while(hx<0x00100000) {             /* normalize x */
00158            hx = hx+hx+(lx>>31); lx = lx+lx;
00159            iy -= 1;
00160        }
00161        if(iy>= -1022) {     /* normalize output */
00162            hx = ((hx-0x00100000)|((iy+1023)<<20));
00163             ux.d = x;
00164            __HI(ux) = hx|sx;
00165            __LO(ux) = lx;
00166             x = ux.d;
00167        } else {             /* subnormal output */
00168            n = -1022 - iy;
00169            if(n<=20) {
00170               lx = (lx>>n)|((unsigned)hx<<(32-n));
00171               hx >>= n;
00172            } else if (n<=31) {
00173               lx = (hx<<(32-n))|(lx>>n); hx = sx;
00174            } else {
00175               lx = hx>>(n-32); hx = sx;
00176            }
00177             ux.d = x;
00178            __HI(ux) = hx|sx;
00179            __LO(ux) = lx;
00180             x = ux.d;
00181            x *= one;        /* create necessary signal */
00182        }
00183        return x;            /* exact output */
00184 }