Home | History | Annotate | Line # | Download | only in math 
      1 /* e_fmodl.c -- long double version of e_fmod.c.
      2  * Conversion to IEEE quad long double by Jakub Jelinek, jj (at) ultra.linux.cz.
      3  */
      4 /*
      5  * ====================================================
      6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      7  *
      8  * Developed at SunPro, a Sun Microsystems, Inc. business.
      9  * Permission to use, copy, modify, and distribute this
     10  * software is freely granted, provided that this notice
     11  * is preserved.
     12  * ====================================================
     13  */
     14 
     15 /* remainderq(x,p)
     16  * Return :
     17  *	returns  x REM p  =  x - [x/p]*p as if in infinite
     18  *	precise arithmetic, where [x/p] is the (infinite bit)
     19  *	integer nearest x/p (in half way case choose the even one).
     20  * Method :
     21  *	Based on fmodl() return x-[x/p]chopped*p exactlp.
     22  */
     23 
     24 #include "quadmath-imp.h"
     25 
     26 static const __float128 zero = 0;
     27 
     28 
     29 __float128
     30 remainderq(__float128 x, __float128 p)
     31 {
     32 	int64_t hx,hp;
     33 	uint64_t sx,lx,lp;
     34 	__float128 p_half;
     35 
     36 	GET_FLT128_WORDS64(hx,lx,x);
     37 	GET_FLT128_WORDS64(hp,lp,p);
     38 	sx = hx&0x8000000000000000ULL;
     39 	hp &= 0x7fffffffffffffffLL;
     40 	hx &= 0x7fffffffffffffffLL;
     41 
     42     /* purge off exception values */
     43 	if((hp|lp)==0) return (x*p)/(x*p);	/* p = 0 */
     44 	if((hx>=0x7fff000000000000LL)||			/* x not finite */
     45 	  ((hp>=0x7fff000000000000LL)&&			/* p is NaN */
     46 	  (((hp-0x7fff000000000000LL)|lp)!=0)))
     47 	    return (x*p)/(x*p);
     48 
     49 
     50 	if (hp<=0x7ffdffffffffffffLL) x = fmodq(x,p+p);	/* now x < 2p */
     51 	if (((hx-hp)|(lx-lp))==0) return zero*x;
     52 	x  = fabsq(x);
     53 	p  = fabsq(p);
     54 	if (hp<0x0002000000000000LL) {
     55 	    if(x+x>p) {
     56 		x-=p;
     57 		if(x+x>=p) x -= p;
     58 	    }
     59 	} else {
     60 	    p_half = 0.5Q*p;
     61 	    if(x>p_half) {
     62 		x-=p;
     63 		if(x>=p_half) x -= p;
     64 	    }
     65 	}
     66 	GET_FLT128_MSW64(hx,x);
     67 	SET_FLT128_MSW64(x,hx^sx);
     68 	return x;
     69 }
     70