text: /*
   * IBM Accurate Mathematical Library
   * written by International Business Machines Corp.
   * Copyright (C) 2001, 2011 Free Software Foundation
   *
   * This program is free software; you can redistribute it and/or modify
   * it under the terms of the GNU Lesser General Public License as published by
   * the Free Software Foundation; either version 2.1 of the License, or
   * (at your option) any later version.
   *
   * This program is distributed in the hope that it will be useful,
   * but WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   * GNU Lesser General Public License for more details.
   *
   * You should have received a copy of the GNU Lesser General Public License
   * along with this program; if not, write to the Free Software
   * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
   */
  /**************************************************************************/
  /*  MODULE_NAME urem.c                                                    */
  /*                                                                        */
  /*  FUNCTION: uremainder                                                  */
  /*                                                                        */
  /* An ultimate remainder routine. Given two IEEE double machine numbers x */
  /* ,y   it computes the correctly rounded (to nearest) value of remainder */
  /* of dividing x by y.                                                    */
  /* Assumption: Machine arithmetic operations are performed in             */
  /* round to nearest mode of IEEE 754 standard.                            */
  /*                                                                        */
  /* ************************************************************************/
  
  #include "endian.h"
  #include "mydefs.h"
  #include "urem.h"
  #include "MathLib.h"
  #include "math_private.h"
  
  /**************************************************************************/
  /* An ultimate remainder routine. Given two IEEE double machine numbers x */
  /* ,y   it computes the correctly rounded (to nearest) value of remainder */
  /**************************************************************************/
  double __ieee754_remainder(double x, double y)
  {
    double z,d,xx;
  #if 0
    double yy;
  #endif
    int4 kx,ky,n,nn,n1,m1,l;
  #if 0
    int4 m;
  #endif
    mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
    u.x=x;
    t.x=y;
    kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
    t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
    ky=t.i[HIGH_HALF];
    /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
    if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
      if (kx+0x00100000<ky) return x;
      if ((kx-0x01500000)<ky) {
        z=x/t.x;
        v.i[HIGH_HALF]=t.i[HIGH_HALF];
        d=(z+big.x)-big.x;
        xx=(x-d*v.x)-d*(t.x-v.x);
        if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
        else {
  	if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
  	else return xx;
        }
      }   /*    (kx<(ky+0x01500000))         */
      else  {
        r.x=1.0/t.x;
        n=t.i[HIGH_HALF];
        nn=(n&0x7ff00000)+0x01400000;
        w.i[HIGH_HALF]=n;
        ww.x=t.x-w.x;
        l=(kx-nn)&0xfff00000;
        n1=ww.i[HIGH_HALF];
        m1=r.i[HIGH_HALF];
        while (l>0) {
  	r.i[HIGH_HALF]=m1-l;
  	z=u.x*r.x;
  	w.i[HIGH_HALF]=n+l;
  	ww.i[HIGH_HALF]=(n1)?n1+l:n1;
  	d=(z+big.x)-big.x;
  	u.x=(u.x-d*w.x)-d*ww.x;
  	l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
        }
        r.i[HIGH_HALF]=m1;
        w.i[HIGH_HALF]=n;
        ww.i[HIGH_HALF]=n1;
        z=u.x*r.x;
        d=(z+big.x)-big.x;
        u.x=(u.x-d*w.x)-d*ww.x;
        if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
        else
  	if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
  	else
  	{z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
      }
  
    }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
    else {
      if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
        y=ABS(y)*t128.x;
        z=__ieee754_remainder(x,y)*t128.x;
        z=__ieee754_remainder(z,y)*tm128.x;
        return z;
      }
    else {
      if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
        y=ABS(y);
        z=2.0*__ieee754_remainder(0.5*x,y);
        d = ABS(z);
        if (d <= ABS(d-y)) return z;
        else return (z>0)?z-y:z+y;
      }
      else { /* if x is too big */
        if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
  	return x / x;
        if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
  	  (ky==0x7ff00000&&t.i[LOW_HALF]!=0))
  	return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
        else return x;
      }
     }
    }
  }
  strong_alias (__ieee754_remainder, __remainder_finite)
  