home *** CD-ROM | disk | FTP | other *** search
/ Chip 2001 Mobile / Chip_Mobile_2001.iso / palm / hobby / palmoon / palmoon.EXE / e_remainder.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-08-16  |  2.0 KB  |  83 lines
  1. // 15 August 1997, Rick Huebner:  Small changes made to adapt for MathLib
  2.  
  3. /* @(#)e_remainder.c 5.1 93/09/24 */
  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. #if defined(LIBM_SCCS) && !defined(lint)
  16. static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp $";
  17. #endif
  18.  
  19. /* __ieee754_remainder(x,p)
  20.  * Return :                  
  21.  *     returns  x REM p  =  x - [x/p]*p as if in infinite 
  22.  *     precise arithmetic, where [x/p] is the (infinite bit) 
  23.  *    integer nearest x/p (in half way case choose the even one).
  24.  * Method : 
  25.  *    Based on fmod() return x-[x/p]chopped*p exactlp.
  26.  */
  27.  
  28. #include "math.h"
  29. #include "math_private.h"
  30.  
  31. #ifdef __STDC__
  32. static const double zero = 0.0;
  33. #else
  34. static double zero = 0.0;
  35. #endif
  36.  
  37.  
  38. #ifdef __STDC__
  39.     double __ieee754_remainder(double x, double p)
  40. #else
  41.     double __ieee754_remainder(x,p)
  42.     double x,p;
  43. #endif
  44. {
  45.     int32_t hx,hp;
  46.     u_int32_t sx,lx,lp;
  47.     double p_half;
  48.  
  49.     EXTRACT_WORDS(hx,lx,x);
  50.     EXTRACT_WORDS(hp,lp,p);
  51.     sx = hx&0x80000000;
  52.     hp &= 0x7fffffff;
  53.     hx &= 0x7fffffff;
  54.  
  55.     /* purge off exception values */
  56.     if((hp|lp)==0) return (x*p)/(x*p);     /* p = 0 */
  57.     if((hx>=0x7ff00000)||            /* x not finite */
  58.       ((hp>=0x7ff00000)&&            /* p is NaN */
  59.       (((hp-0x7ff00000)|lp)!=0)))
  60.         return (x*p)/(x*p);
  61.  
  62.  
  63.     if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);    /* now x < 2p */
  64.     if (((hx-hp)|(lx-lp))==0) return zero*x;
  65.     x  = jumpto__fabs(x);
  66.     p  = jumpto__fabs(p);
  67.     if (hp<0x00200000) {
  68.         if(x+x>p) {
  69.         x-=p;
  70.         if(x+x>=p) x -= p;
  71.         }
  72.     } else {
  73.         p_half = 0.5*p;
  74.         if(x>p_half) {
  75.         x-=p;
  76.         if(x>=p_half) x -= p;
  77.         }
  78.     }
  79.     GET_HIGH_WORD(hx,x);
  80.     SET_HIGH_WORD(x,hx^sx);
  81.     return x;
  82. }
  83.