Home | History | Annotate | Line # | Download | only in src 
      1 /*	$NetBSD: s_nextafterl.c,v 1.7 2026/04/02 22:02:10 nat 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_nextafterl.c,v 1.7 2026/04/02 22:02:10 nat Exp $");
     17 
     18 #include <float.h>
     19 #include <math.h>
     20 #include <machine/ieee.h>
     21 
     22 #ifdef __HAVE_LONG_DOUBLE
     23 
     24 #ifdef EXT_EXP_INFNAN
     25 #if LDBL_MAX_EXP != 0x4000
     26 #error "Unsupported long double format"
     27 #endif
     28 
     29 #ifdef LDBL_IMPLICIT_NBIT
     30 #define	LDBL_NBIT	0
     31 #endif
     32 
     33 __strong_alias(nexttowardl, nextafterl)
     34 
     35 /*
     36  * IEEE functions
     37  *      nextafterl(x,y)
     38  *      return the next machine floating-point number of x in the
     39  *      direction toward y.
     40  *   Special cases:
     41  *	If x == y, y shall be returned
     42  *	If x or y is NaN, a NaN shall be returned
     43  */
     44 long double
     45 nextafterl(long double x, long double y)
     46 {
     47 	volatile long double t;
     48 	union ieee_ext_u ux, uy;
     49 
     50 	ux.extu_ld = x;
     51 	uy.extu_ld = y;
     52 
     53 	if ((ux.extu_exp == EXT_EXP_INFNAN &&
     54 		((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
     55 	    (uy.extu_exp == EXT_EXP_INFNAN &&
     56 		((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
     57 		return x+y;			/* x or y is nan */
     58 
     59 	if (x == y) return y;			/* x=y, return y */
     60 
     61 	if (x == 0.0) {
     62 		ux.extu_frach = 0;		/* return +-minsubnormal */
     63 		ux.extu_fracl = 1;
     64 		ux.extu_sign = uy.extu_sign;
     65 		t = ux.extu_ld * ux.extu_ld;
     66 		if (t == ux.extu_ld)
     67 			return t;
     68 		else
     69 			return ux.extu_ld;	/* raise underflow flag */
     70 	}
     71 
     72 	uint32_t nbit = (ux.extu_frach & LDBL_NBIT);
     73 	if ((x>0.0) ^ (x<y)) {			/* x -= ulp */
     74 		if (ux.extu_fracl == 0) {
     75 			if ((ux.extu_frach & ~LDBL_NBIT) == 0) {
     76 #ifdef __m68k__
     77 				if (ux.extu_exp == 0)
     78 					nbit = 0;
     79 				else
     80 #endif
     81 					ux.extu_exp -= 1;
     82 			}
     83 			ux.extu_frach = (ux.extu_frach - 1) | nbit;
     84 		}
     85 		ux.extu_fracl -= 1;
     86 	} else {				/* x += ulp */
     87 		ux.extu_fracl += 1;
     88 		if (ux.extu_fracl == 0) {
     89 			ux.extu_frach = (ux.extu_frach + 1) | nbit;
     90 			if ((ux.extu_frach & ~LDBL_NBIT) == 0) {
     91 #ifdef __m68k__
     92 				if (nbit == 0)
     93 					ux.extu_frach |= LDBL_NBIT;
     94 				else
     95 #endif
     96 					ux.extu_exp += 1;
     97 			}
     98 		}
     99 	}
    100 
    101 	if (ux.extu_exp == EXT_EXP_INFNAN)
    102 		return x+x;			/* overflow  */
    103 
    104 #ifdef __m68k__
    105 	if (ux.extu_exp == 0 && (ux.extu_frach & LDBL_NBIT) == 0) {
    106 #else
    107 	if (ux.extu_exp == 0) {			/* underflow */
    108 #endif
    109 #ifndef LDBL_IMPLICIT_NBIT
    110 		mask_nbit_l(ux);
    111 #endif
    112 		t = ux.extu_ld * ux.extu_ld;
    113 		if (t != ux.extu_ld)		/* raise underflow flag */
    114 			return ux.extu_ld;
    115 	}
    116 
    117 	return ux.extu_ld;
    118 }
    119 #endif
    120 
    121 #endif /* __HAVE_LONG_DOUBLE */
    122