1143217Sdas/* @(#)s_nextafter.c 5.1 93/09/24 */
2143217Sdas/*
3143217Sdas * ====================================================
4143217Sdas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5143217Sdas *
6143217Sdas * Developed at SunPro, a Sun Microsystems, Inc. business.
7143217Sdas * Permission to use, copy, modify, and distribute this
8143217Sdas * software is freely granted, provided that this notice
9143217Sdas * is preserved.
10143217Sdas * ====================================================
11143217Sdas */
12143217Sdas
13176451Sdas#include <sys/cdefs.h>
14176451Sdas__FBSDID("$FreeBSD$");
15143217Sdas
16143217Sdas/* IEEE functions
17143217Sdas *	nextafter(x,y)
18143217Sdas *	return the next machine floating-point number of x in the
19143217Sdas *	direction toward y.
20143217Sdas *   Special cases:
21143217Sdas */
22143217Sdas
23143217Sdas#include <float.h>
24143217Sdas
25143217Sdas#include "fpmath.h"
26143217Sdas#include "math.h"
27143217Sdas#include "math_private.h"
28143217Sdas
29143217Sdas#if LDBL_MAX_EXP != 0x4000
30143217Sdas#error "Unsupported long double format"
31143217Sdas#endif
32143217Sdas
33143217Sdaslong double
34143217Sdasnextafterl(long double x, long double y)
35143217Sdas{
36143217Sdas	volatile long double t;
37143217Sdas	union IEEEl2bits ux, uy;
38143217Sdas
39143217Sdas	ux.e = x;
40143217Sdas	uy.e = y;
41143217Sdas
42143217Sdas	if ((ux.bits.exp == 0x7fff &&
43143217Sdas	     ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) ||
44143217Sdas	    (uy.bits.exp == 0x7fff &&
45143217Sdas	     ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
46143217Sdas	   return x+y;	/* x or y is nan */
47143217Sdas	if(x==y) return y;		/* x=y, return y */
48143217Sdas	if(x==0.0) {
49143217Sdas	    ux.bits.manh = 0;			/* return +-minsubnormal */
50143217Sdas	    ux.bits.manl = 1;
51143217Sdas	    ux.bits.sign = uy.bits.sign;
52143217Sdas	    t = ux.e*ux.e;
53143217Sdas	    if(t==ux.e) return t; else return ux.e; /* raise underflow flag */
54143217Sdas	}
55143217Sdas	if(x>0.0 ^ x<y) {			/* x -= ulp */
56143217Sdas	    if(ux.bits.manl==0) {
57143217Sdas		if ((ux.bits.manh&~LDBL_NBIT)==0)
58143217Sdas		    ux.bits.exp -= 1;
59143217Sdas		ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
60143217Sdas	    }
61143217Sdas	    ux.bits.manl -= 1;
62143217Sdas	} else {				/* x += ulp */
63143217Sdas	    ux.bits.manl += 1;
64143217Sdas	    if(ux.bits.manl==0) {
65143217Sdas		ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
66143217Sdas		if ((ux.bits.manh&~LDBL_NBIT)==0)
67143217Sdas		    ux.bits.exp += 1;
68143217Sdas	    }
69143217Sdas	}
70143217Sdas	if(ux.bits.exp==0x7fff) return x+x;	/* overflow  */
71143217Sdas	if(ux.bits.exp==0) {			/* underflow */
72143217Sdas	    mask_nbit_l(ux);
73143217Sdas	    t = ux.e * ux.e;
74143217Sdas	    if(t!=ux.e)			/* raise underflow flag */
75143217Sdas		return ux.e;
76143217Sdas	}
77143217Sdas	return ux.e;
78143217Sdas}
79143217Sdas
80143217Sdas__strong_reference(nextafterl, nexttowardl);
81