Home | History | Annotate | Line # | Download | only in src 
      1 /*	$NetBSD: s_nexttoward.c,v 1.3 2024/05/11 02:07:53 riastradh Exp $	*/
      2 
      3 /* @(#)s_nextafter.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 #include <sys/cdefs.h>
     16 __RCSID("$NetBSD: s_nexttoward.c,v 1.3 2024/05/11 02:07:53 riastradh Exp $");
     17 
     18 /*
     19  * We assume that a long double has a 15-bit exponent.  On systems
     20  * where long double is the same as double, nexttoward() is an alias
     21  * for nextafter(), so we don't use this routine.
     22  */
     23 #include <float.h>
     24 
     25 #include <machine/ieee.h>
     26 #include "math.h"
     27 #include "math_private.h"
     28 
     29 #if LDBL_MAX_EXP != 0x4000
     30 #error "Unsupported long double format"
     31 #endif
     32 
     33 #ifdef LDBL_IMPLICIT_NBIT
     34 #define	LDBL_NBIT	0
     35 #endif
     36 
     37 /*
     38  * The nexttoward() function is equivalent to nextafter() function,
     39  * except that the second parameter shall have type long double and
     40  * the functions shall return y converted to the type of the function
     41  * if x equals y.
     42  *
     43  * Special cases: XXX
     44  */
     45 double
     46 nexttoward(double x, long double y)
     47 {
     48 	union ieee_ext_u uy;
     49 	volatile double t;
     50 	int32_t hx, ix;
     51 	uint32_t lx;
     52 
     53 	EXTRACT_WORDS(hx, lx, x);
     54 	ix = hx & 0x7fffffff;			/* |x| */
     55 	uy.extu_ld = y;
     56 
     57 	if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
     58 	    (uy.extu_exp == 0x7fff &&
     59 		((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
     60 		return x+y;			/* x or y is nan */
     61 
     62 	if (x == y)
     63 		return (double)y;		/* x=y, return y */
     64 
     65 	if (x == 0.0) {
     66 		INSERT_WORDS(x, uy.extu_sign<<31, 1);	/* return +-minsubnormal */
     67 		t = x*x;
     68 		if (t == x)
     69 			return t;
     70 		else
     71 			return x;		/* raise underflow flag */
     72 	}
     73 
     74 	if ((hx >= 0) ^ (x < y)) {		/* x -= ulp */
     75 		if (lx == 0) hx -= 1;
     76 		lx -= 1;
     77 	} else {				/* x += ulp */
     78 		lx += 1;
     79 		if (lx == 0) hx += 1;
     80 	}
     81 	ix = hx & 0x7ff00000;
     82 	if (ix >= 0x7ff00000) return x+x;	/* overflow  */
     83 	if (ix <  0x00100000) {			/* underflow */
     84 		t = x*x;
     85 		if (t != x) {			/* raise underflow flag */
     86 			INSERT_WORDS(y, hx, lx);
     87 			return y;
     88 		}
     89 	}
     90 	INSERT_WORDS(x, hx, lx);
     91 
     92 	return x;
     93 }
     94