source: trunk/sys/libm/s_nextafter.c @ 9

Last change on this file since 9 was 1, checked in by alain, 7 years ago

First import

File size: 2.2 KB
Line 
1
2/* @(#)s_nextafter.c 5.1 93/09/24 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 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/* IEEE functions
15 *      nextafter(x,y)
16 *      return the next machine floating-point number of x in the
17 *      direction toward y.
18 *   Special cases:
19 */
20
21#include <libm/fdlibm.h>
22
23#ifdef __STDC__
24static const double one=1.0;
25#else
26static double one=1.0;
27#endif
28
29#ifdef __STDC__
30        double nextafter(double x, double y)
31#else
32        double nextafter(x,y)
33        double x,y;
34#endif
35{
36        int     n0,n1,hx,hy,ix,iy;
37        unsigned lx,ly;
38
39        n0 = ((*(int*)&one)>>29)^1;     /* index of high word */
40        n1 = 1-n0;
41        hx = *( n0 + (int*)&x);         /* high word of x */
42        lx = *( n1 + (int*)&x);         /* low  word of x */
43        hy = *( n0 + (int*)&y);         /* high word of y */
44        ly = *( n1 + (int*)&y);         /* low  word of y */
45        ix = hx&0x7fffffff;             /* |x| */
46        iy = hy&0x7fffffff;             /* |y| */
47
48        if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */ 
49           ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */ 
50           return x+y;                         
51        if(x==y) return x;              /* x=y, return x */
52        if((ix|lx)==0) {                        /* x == 0 */
53            *(n0+(int*)&x) = hy&0x80000000;     /* return +-minsubnormal */
54            *(n1+(int*)&x) = 1;
55            y = x*x;
56            if(y==x) return y; else return x;   /* raise underflow flag */
57        } 
58        if(hx>=0) {                             /* x > 0 */
59            if(hx>hy||((hx==hy)&&(lx>ly))) {    /* x > y, x -= ulp */
60                if(lx==0) hx -= 1;
61                lx -= 1;
62            } else {                            /* x < y, x += ulp */
63                lx += 1;
64                if(lx==0) hx += 1;
65            }
66        } else {                                /* x < 0 */
67            if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
68                if(lx==0) hx -= 1;
69                lx -= 1;
70            } else {                            /* x > y, x += ulp */
71                lx += 1;
72                if(lx==0) hx += 1;
73            }
74        }
75        hy = hx&0x7ff00000;
76        if(hy>=0x7ff00000) return x+x;  /* overflow  */
77        if(hy<0x00100000) {             /* underflow */
78            y = x*x;
79            if(y!=x) {          /* raise underflow flag */
80                *(n0+(int*)&y) = hx; *(n1+(int*)&y) = lx;
81                return y;
82            }
83        }
84        *(n0+(int*)&x) = hx; *(n1+(int*)&x) = lx;
85        return x;
86}
Note: See TracBrowser for help on using the repository browser.