Home | History | Annotate | Line # | Download | only in math
fmodq.c revision 1.1 
      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, 2011 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 /*
     16  * fmodq(x,y)
     17  * Return x mod y in exact arithmetic
     18  * Method: shift and subtract
     19  */
     20 
     21 #include "quadmath-imp.h"
     22 
     23 static const __float128 one = 1.0, Zero[] = {0.0, -0.0,};
     24 
     25 __float128
     26 fmodq (__float128 x, __float128 y)
     27 {
     28 	int64_t n,hx,hy,hz,ix,iy,sx,i;
     29 	uint64_t lx,ly,lz;
     30 
     31 	GET_FLT128_WORDS64(hx,lx,x);
     32 	GET_FLT128_WORDS64(hy,ly,y);
     33 	sx = hx&0x8000000000000000ULL;		/* sign of x */
     34 	hx ^=sx;				/* |x| */
     35 	hy &= 0x7fffffffffffffffLL;		/* |y| */
     36 
     37     /* purge off exception values */
     38 	if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
     39 	  ((hy|((ly|-ly)>>63))>0x7fff000000000000LL))	/* or y is NaN */
     40 	    return (x*y)/(x*y);
     41 	if(hx<=hy) {
     42 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
     43 	    if(lx==ly)
     44 		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
     45 	}
     46 
     47     /* determine ix = ilogb(x) */
     48 	if(hx<0x0001000000000000LL) {	/* subnormal x */
     49 	    if(hx==0) {
     50 		for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
     51 	    } else {
     52 		for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
     53 	    }
     54 	} else ix = (hx>>48)-0x3fff;
     55 
     56     /* determine iy = ilogb(y) */
     57 	if(hy<0x0001000000000000LL) {	/* subnormal y */
     58 	    if(hy==0) {
     59 		for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
     60 	    } else {
     61 		for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
     62 	    }
     63 	} else iy = (hy>>48)-0x3fff;
     64 
     65     /* set up {hx,lx}, {hy,ly} and align y to x */
     66 	if(ix >= -16382)
     67 	    hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
     68 	else {		/* subnormal x, shift x to normal */
     69 	    n = -16382-ix;
     70 	    if(n<=63) {
     71 		hx = (hx<<n)|(lx>>(64-n));
     72 		lx <<= n;
     73 	    } else {
     74 		hx = lx<<(n-64);
     75 		lx = 0;
     76 	    }
     77 	}
     78 	if(iy >= -16382)
     79 	    hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
     80 	else {		/* subnormal y, shift y to normal */
     81 	    n = -16382-iy;
     82 	    if(n<=63) {
     83 		hy = (hy<<n)|(ly>>(64-n));
     84 		ly <<= n;
     85 	    } else {
     86 		hy = ly<<(n-64);
     87 		ly = 0;
     88 	    }
     89 	}
     90 
     91     /* fix point fmod */
     92 	n = ix - iy;
     93 	while(n--) {
     94 	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
     95 	    if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
     96 	    else {
     97 		if((hz|lz)==0)		/* return sign(x)*0 */
     98 		    return Zero[(uint64_t)sx>>63];
     99 		hx = hz+hz+(lz>>63); lx = lz+lz;
    100 	    }
    101 	}
    102 	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
    103 	if(hz>=0) {hx=hz;lx=lz;}
    104 
    105     /* convert back to floating value and restore the sign */
    106 	if((hx|lx)==0)			/* return sign(x)*0 */
    107 	    return Zero[(uint64_t)sx>>63];
    108 	while(hx<0x0001000000000000LL) {	/* normalize x */
    109 	    hx = hx+hx+(lx>>63); lx = lx+lx;
    110 	    iy -= 1;
    111 	}
    112 	if(iy>= -16382) {	/* normalize output */
    113 	    hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
    114 	    SET_FLT128_WORDS64(x,hx|sx,lx);
    115 	} else {		/* subnormal output */
    116 	    n = -16382 - iy;
    117 	    if(n<=48) {
    118 		lx = (lx>>n)|((uint64_t)hx<<(64-n));
    119 		hx >>= n;
    120 	    } else if (n<=63) {
    121 		lx = (hx<<(64-n))|(lx>>n); hx = sx;
    122 	    } else {
    123 		lx = hx>>(n-64); hx = sx;
    124 	    }
    125 	    SET_FLT128_WORDS64(x,hx|sx,lx);
    126 	    x *= one;		/* create necessary signal */
    127 	}
    128 	return x;		/* exact output */
    129 }
    130