gonzui

    1: /* @(#)s_nextafter.c 5.1 93/09/24 */
    2: /*
    3:  * ====================================================
    4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    5:  *
    6:  * Developed at SunPro, a Sun Microsystems, Inc. business.
    7:  * Permission to use, copy, modify, and distribute this
    8:  * software is freely granted, provided that this notice 
    9:  * is preserved.
   10:  * ====================================================
   11:  */
   12: 
   13: /* LINTLIBRARY */
   14: 
   15: /* IEEE functions
   16:  *      nextafter(x,y)
   17:  *      return the next machine floating-point number of x in the
   18:  *      direction toward y.
   19:  *   Special cases:
   20:  */
   21: 
   22: #include <sys/cdefs.h>
   23: #include <float.h>
   24: #include <math.h>
   25: 
   26: #include "math_private.h"
   27: 
   28: double
   29: nextafter(double x, double y)
   30: {
   31:         int32_t hx,hy,ix,iy;
   32:         u_int32_t lx,ly;
   33: 
   34:         EXTRACT_WORDS(hx,lx,x);
   35:         EXTRACT_WORDS(hy,ly,y);
   36:         ix = hx&0x7fffffff;            /* |x| */
   37:         iy = hy&0x7fffffff;            /* |y| */
   38: 
   39:         if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
   40:            ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */
   41:            return x+y;                         
   42:         if(x==y) return y;             /* x=y, return y */
   43:         if((ix|lx)==0) {                       /* x == 0 */
   44:             INSERT_WORDS(x,hy&0x80000000,1);   /* return +-minsubnormal */
   45:             y = x*x;
   46:             if(y==x) return y; else return x;  /* raise underflow flag */
   47:         }
   48:         if(hx>=0) {                            /* x > 0 */
   49:             if(hx>hy||((hx==hy)&&(lx>ly))) {   /* x > y, x -= ulp */
   50:                 if(lx==0) hx -= 1;
   51:                 lx -= 1;
   52:             } else {                           /* x < y, x += ulp */
   53:                 lx += 1;
   54:                 if(lx==0) hx += 1;
   55:             }
   56:         } else {                               /* x < 0 */
   57:             if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
   58:                 if(lx==0) hx -= 1;
   59:                 lx -= 1;
   60:             } else {                           /* x > y, x += ulp */
   61:                 lx += 1;
   62:                 if(lx==0) hx += 1;
   63:             }
   64:         }
   65:         hy = hx&0x7ff00000;
   66:         if(hy>=0x7ff00000) return x+x; /* overflow  */
   67:         if(hy<0x00100000) {            /* underflow */
   68:             y = x*x;
   69:             if(y!=x) {         /* raise underflow flag */
   70:                 INSERT_WORDS(y,hx,lx);
   71:                 return y;
   72:             }
   73:         }
   74:         INSERT_WORDS(x,hx,lx);
   75:         return x;
   76: }
   77: 
   78: #if     LDBL_MANT_DIG == 53
   79: #ifdef  lint
   80: /* PROTOLIB1 */
   81: long double nextafterl(long double, long double);
   82: double nexttoward(double, long double);
   83: long double nexttowardl(long double, long double);
   84: #else   /* lint */
   85: __weak_alias(nextafterl, nextafter);
   86: __weak_alias(nexttoward, nextafter);
   87: __weak_alias(nexttowardl, nextafter);
   88: #endif  /* lint */
   89: #endif  /* LDBL_MANT_DIG == 53 */