1  /* e_fmodl.c -- long double version of e_fmod.c.
       2   */
       3  /*
       4   * ====================================================
       5   * Copyright (C) 1993, 2011 by Sun Microsystems, Inc. All rights reserved.
       6   *
       7   * Developed at SunPro, a Sun Microsystems, Inc. business.
       8   * Permission to use, copy, modify, and distribute this
       9   * software is freely granted, provided that this notice
      10   * is preserved.
      11   * ====================================================
      12   */
      13  
      14  /*
      15   * __ieee754_fmodl(x,y)
      16   * Return x mod y in exact arithmetic
      17   * Method: shift and subtract
      18   */
      19  
      20  #include <math.h>
      21  #include <math_private.h>
      22  #include <libm-alias-finite.h>
      23  
      24  static const _Float128 one = 1.0, Zero[] = {0.0, -0.0,};
      25  
      26  _Float128
      27  __ieee754_fmodl (_Float128 x, _Float128 y)
      28  {
      29  	int64_t n,hx,hy,hz,ix,iy,sx,i;
      30  	uint64_t lx,ly,lz;
      31  
      32  	GET_LDOUBLE_WORDS64(hx,lx,x);
      33  	GET_LDOUBLE_WORDS64(hy,ly,y);
      34  	sx = hx&0x8000000000000000ULL;		/* sign of x */
      35  	hx ^=sx;				/* |x| */
      36  	hy &= 0x7fffffffffffffffLL;		/* |y| */
      37  
      38      /* purge off exception values */
      39  	if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
      40  	  ((hy|((ly|-ly)>>63))>0x7fff000000000000LL))	/* or y is NaN */
      41  	    return (x*y)/(x*y);
      42  	if(hx<=hy) {
      43  	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
      44  	    if(lx==ly)
      45  		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
      46  	}
      47  
      48      /* determine ix = ilogb(x) */
      49  	if(hx<0x0001000000000000LL) {	/* subnormal x */
      50  	    if(hx==0) {
      51  		for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
      52  	    } else {
      53  		for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
      54  	    }
      55  	} else ix = (hx>>48)-0x3fff;
      56  
      57      /* determine iy = ilogb(y) */
      58  	if(hy<0x0001000000000000LL) {	/* subnormal y */
      59  	    if(hy==0) {
      60  		for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
      61  	    } else {
      62  		for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
      63  	    }
      64  	} else iy = (hy>>48)-0x3fff;
      65  
      66      /* set up {hx,lx}, {hy,ly} and align y to x */
      67  	if(ix >= -16382)
      68  	    hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
      69  	else {		/* subnormal x, shift x to normal */
      70  	    n = -16382-ix;
      71  	    if(n<=63) {
      72  		hx = (hx<<n)|(lx>>(64-n));
      73  		lx <<= n;
      74  	    } else {
      75  		hx = lx<<(n-64);
      76  		lx = 0;
      77  	    }
      78  	}
      79  	if(iy >= -16382)
      80  	    hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
      81  	else {		/* subnormal y, shift y to normal */
      82  	    n = -16382-iy;
      83  	    if(n<=63) {
      84  		hy = (hy<<n)|(ly>>(64-n));
      85  		ly <<= n;
      86  	    } else {
      87  		hy = ly<<(n-64);
      88  		ly = 0;
      89  	    }
      90  	}
      91  
      92      /* fix point fmod */
      93  	n = ix - iy;
      94  	while(n--) {
      95  	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
      96  	    if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
      97  	    else {
      98  		if((hz|lz)==0)		/* return sign(x)*0 */
      99  		    return Zero[(uint64_t)sx>>63];
     100  		hx = hz+hz+(lz>>63); lx = lz+lz;
     101  	    }
     102  	}
     103  	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
     104  	if(hz>=0) {hx=hz;lx=lz;}
     105  
     106      /* convert back to floating value and restore the sign */
     107  	if((hx|lx)==0)			/* return sign(x)*0 */
     108  	    return Zero[(uint64_t)sx>>63];
     109  	while(hx<0x0001000000000000LL) {	/* normalize x */
     110  	    hx = hx+hx+(lx>>63); lx = lx+lx;
     111  	    iy -= 1;
     112  	}
     113  	if(iy>= -16382) {	/* normalize output */
     114  	    hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
     115  	    SET_LDOUBLE_WORDS64(x,hx|sx,lx);
     116  	} else {		/* subnormal output */
     117  	    n = -16382 - iy;
     118  	    if(n<=48) {
     119  		lx = (lx>>n)|((uint64_t)hx<<(64-n));
     120  		hx >>= n;
     121  	    } else if (n<=63) {
     122  		lx = (hx<<(64-n))|(lx>>n); hx = sx;
     123  	    } else {
     124  		lx = hx>>(n-64); hx = sx;
     125  	    }
     126  	    SET_LDOUBLE_WORDS64(x,hx|sx,lx);
     127  	    x *= one;		/* create necessary signal */
     128  	}
     129  	return x;		/* exact output */
     130  }
     131  libm_alias_finite (__ieee754_fmodl, __fmodl)