e_remainder.c - vx32 - Local 9vx git repository for patches.
 (HTM) git clone git://r-36.net/vx32
 (DIR) Log
 (DIR) Files
 (DIR) Refs
       ---
       e_remainder.c (1713B)
       ---
            1 /* @(#)e_remainder.c 5.1 93/09/24 */
            2 /*
            3  * ====================================================
            4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
            5  *
            6  * Developed at SunPro, 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 #ifndef lint
           14 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_remainder.c,v 1.9 2003/07/23 04:53:46 peter Exp $";
           15 #endif
           16 
           17 /* __ieee754_remainder(x,p)
           18  * Return :
           19  *         returns  x REM p  =  x - [x/p]*p as if in infinite
           20  *         precise arithmetic, where [x/p] is the (infinite bit)
           21  *        integer nearest x/p (in half way case choose the even one).
           22  * Method :
           23  *        Based on fmod() return x-[x/p]chopped*p exactlp.
           24  */
           25 
           26 #include "math.h"
           27 #include "math_private.h"
           28 
           29 static const double zero = 0.0;
           30 
           31 
           32 double
           33 __ieee754_remainder(double x, double p)
           34 {
           35         int32_t hx,hp;
           36         u_int32_t sx,lx,lp;
           37         double p_half;
           38 
           39         EXTRACT_WORDS(hx,lx,x);
           40         EXTRACT_WORDS(hp,lp,p);
           41         sx = hx&0x80000000;
           42         hp &= 0x7fffffff;
           43         hx &= 0x7fffffff;
           44 
           45     /* purge off exception values */
           46         if((hp|lp)==0) return (x*p)/(x*p);         /* p = 0 */
           47         if((hx>=0x7ff00000)||                        /* x not finite */
           48           ((hp>=0x7ff00000)&&                        /* p is NaN */
           49           (((hp-0x7ff00000)|lp)!=0)))
           50             return (x*p)/(x*p);
           51 
           52 
           53         if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);        /* now x < 2p */
           54         if (((hx-hp)|(lx-lp))==0) return zero*x;
           55         x  = fabs(x);
           56         p  = fabs(p);
           57         if (hp<0x00200000) {
           58             if(x+x>p) {
           59                 x-=p;
           60                 if(x+x>=p) x -= p;
           61             }
           62         } else {
           63             p_half = 0.5*p;
           64             if(x>p_half) {
           65                 x-=p;
           66                 if(x>=p_half) x -= p;
           67             }
           68         }
           69         GET_HIGH_WORD(hx,x);
           70         SET_HIGH_WORD(x,hx^sx);
           71         return x;
           72 }