Back to index

glibc  2.9
e_remainderl.c
Go to the documentation of this file.
00001 /* e_fmodl.c -- long double version of e_fmod.c.
00002  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
00003  */
00004 /*
00005  * ====================================================
00006  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00007  *
00008  * Developed at SunPro, a Sun Microsystems, Inc. business.
00009  * Permission to use, copy, modify, and distribute this
00010  * software is freely granted, provided that this notice 
00011  * is preserved.
00012  * ====================================================
00013  */
00014 
00015 /* __ieee754_remainderl(x,p)
00016  * Return :                  
00017  *     returns  x REM p  =  x - [x/p]*p as if in infinite 
00018  *     precise arithmetic, where [x/p] is the (infinite bit) 
00019  *     integer nearest x/p (in half way case choose the even one).
00020  * Method : 
00021  *     Based on fmodl() return x-[x/p]chopped*p exactlp.
00022  */
00023 
00024 #include "math.h"
00025 #include "math_private.h"
00026 
00027 #ifdef __STDC__
00028 static const long double zero = 0.0L;
00029 #else
00030 static long double zero = 0.0L;
00031 #endif
00032 
00033 
00034 #ifdef __STDC__
00035        long double __ieee754_remainderl(long double x, long double p)
00036 #else
00037        long double __ieee754_remainderl(x,p)
00038        long double x,p;
00039 #endif
00040 {
00041        int64_t hx,hp;
00042        u_int64_t sx,lx,lp;
00043        long double p_half;
00044 
00045        GET_LDOUBLE_WORDS64(hx,lx,x);
00046        GET_LDOUBLE_WORDS64(hp,lp,p);
00047        sx = hx&0x8000000000000000ULL;
00048        hp &= 0x7fffffffffffffffLL;
00049        hx &= 0x7fffffffffffffffLL;
00050 
00051     /* purge off exception values */
00052        if((hp|lp)==0) return (x*p)/(x*p);        /* p = 0 */
00053        if((hx>=0x7fff000000000000LL)||                  /* x not finite */
00054          ((hp>=0x7fff000000000000LL)&&                  /* p is NaN */
00055          (((hp-0x7fff000000000000LL)|lp)!=0)))
00056            return (x*p)/(x*p);
00057 
00058 
00059        if (hp<=0x7ffdffffffffffffLL) x = __ieee754_fmodl(x,p+p);      /* now x < 2p */
00060        if (((hx-hp)|(lx-lp))==0) return zero*x;
00061        x  = fabsl(x);
00062        p  = fabsl(p);
00063        if (hp<0x0002000000000000LL) {
00064            if(x+x>p) {
00065               x-=p;
00066               if(x+x>=p) x -= p;
00067            }
00068        } else {
00069            p_half = 0.5L*p;
00070            if(x>p_half) {
00071               x-=p;
00072               if(x>=p_half) x -= p;
00073            }
00074        }
00075        GET_LDOUBLE_MSW64(hx,x);
00076        SET_LDOUBLE_MSW64(x,hx^sx);
00077        return x;
00078 }