s_nextafterl.c revision 1.4 
11.37Sad/*	$NetBSD: s_nextafterl.c,v 1.4 2013/07/18 22:31:13 matt Exp $	*/
21.2Sad
31.2Sad/* @(#)s_nextafter.c 5.1 93/09/24 */
41.37Sad/*
51.28Sad * ====================================================
61.2Sad * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
71.2Sad *
81.2Sad * Developed at SunPro, a Sun Microsystems, Inc. business.
91.2Sad * Permission to use, copy, modify, and distribute this
101.2Sad * software is freely granted, provided that this notice
111.2Sad * is preserved.
121.2Sad * ====================================================
131.2Sad */
141.2Sad
151.2Sad#include <sys/cdefs.h>
161.2Sad__RCSID("$NetBSD: s_nextafterl.c,v 1.4 2013/07/18 22:31:13 matt Exp $");
171.2Sad
181.2Sad#include <float.h>
191.2Sad#include <math.h>
201.2Sad#include <machine/ieee.h>
211.2Sad
221.2Sad#ifdef __HAVE_LONG_DOUBLE
231.2Sad
241.2Sad#ifdef EXT_EXP_INFNAN
251.2Sad#if LDBL_MAX_EXP != 0x4000
261.2Sad#error "Unsupported long double format"
271.2Sad#endif
281.2Sad
291.2Sad#ifdef LDBL_IMPLICIT_NBIT
301.2Sad#define	LDBL_NBIT	0
311.2Sad#endif
321.2Sad
331.2Sad/*
341.2Sad * IEEE functions
351.2Sad *      nextafterl(x,y)
361.15Spooka *      return the next machine floating-point number of x in the
371.15Spooka *      direction toward y.
381.15Spooka *   Special cases:
391.2Sad *	If x == y, y shall be returned
401.2Sad *	If x or y is NaN, a NaN shall be returned
411.4Syamt */
421.27Sadlong double
431.2Sadnextafterl(long double x, long double y)
441.2Sad{
451.2Sad	volatile long double t;
461.2Sad	union ieee_ext_u ux, uy;
471.2Sad
481.13Sad	ux.extu_ld = x;
491.37Sad	uy.extu_ld = y;
501.2Sad
511.13Sad	if ((ux.extu_exp == EXT_EXP_NAN &&
521.17Srmind		((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
531.37Sad	    (uy.extu_exp == EXT_EXP_NAN &&
541.37Sad		((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
551.30Sthorpej		return x+y;			/* x or y is nan */
561.37Sad
571.32Sthorpej	if (x == y) return y;			/* x=y, return y */
581.17Srmind
591.2Sad	if (x == 0.0) {
601.23Spooka		ux.extu_frach = 0;		/* return +-minsubnormal */
611.2Sad		ux.extu_fracl = 1;
621.6Sad		ux.extu_sign = uy.extu_sign;
631.6Sad		t = ux.extu_ld * ux.extu_ld;
641.37Sad		if (t == ux.extu_ld)
651.3Syamt			return t;
661.24Schristos		else
671.36Sriastrad			return ux.extu_ld;	/* raise underflow flag */
681.36Sriastrad	}
691.36Sriastrad
701.27Sad	if ((x>0.0) ^ (x<y)) {			/* x -= ulp */
711.27Sad		if (ux.extu_fracl == 0) {
721.27Sad			if ((ux.extu_frach & ~LDBL_NBIT) == 0)
731.27Sad				ux.extu_exp -= 1;
741.27Sad			ux.extu_frach = (ux.extu_frach - 1) |
751.2Sad					(ux.extu_frach & LDBL_NBIT);
761.2Sad		}
771.2Sad		ux.extu_fracl -= 1;
781.2Sad	} else {				/* x += ulp */
791.2Sad		ux.extu_fracl += 1;
801.25Schristos		if (ux.extu_fracl == 0) {
811.6Sad			ux.extu_frach = (ux.extu_frach + 1) |
821.2Sad					(ux.extu_frach & LDBL_NBIT);
831.2Sad			if ((ux.extu_frach & ~LDBL_NBIT) == 0)
841.7Syamt				ux.extu_exp += 1;
851.2Sad		}
861.2Sad	}
871.37Sad
881.2Sad	if (ux.extu_exp == EXT_EXP_INF)
891.33Schristos		return x+x;			/* overflow  */
901.2Sad
911.2Sad	if (ux.extu_exp == 0) {			/* underflow */
92#ifndef LDBL_IMPLICIT_NBIT
93		mask_nbit_l(ux);
94#endif
95		t = ux.extu_ld * ux.extu_ld;
96		if (t != ux.extu_ld)		/* raise underflow flag */
97			return ux.extu_ld;
98	}
99
100	return ux.extu_ld;
101}
102#endif
103
104#endif /* __HAVE_LONG_DOUBLE */
105