libm/s_nextafterf.c

   1 /* s_nextafterf.c -- float version of s_nextafter.c.
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  */
   4
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #include "math.h"
  17 #include "math_private.h"
  18
  19 float nextafterf(float x, float y)
  20 {
  21         int32_t hx, hy, ix, iy;
  22
  23         GET_FLOAT_WORD(hx, x);
  24         GET_FLOAT_WORD(hy, y);
  25         ix = hx & 0x7fffffff;           /* |x| */
  26         iy = hy & 0x7fffffff;           /* |y| */
  27
  28         /* x is nan or y is nan? */
  29         if ((ix > 0x7f800000) || (iy > 0x7f800000))
  30                 return x + y;
  31
  32         if (x == y)
  33                 return y;
  34
  35         if (ix == 0) { /* x == 0? */
  36 /* glibc 2.4 does not seem to set underflow? */
  37 /*              float u; */
  38                 /* return +-minsubnormal */
  39                 SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);
  40 /*              u = x * x;     raise underflow flag */
  41 /*              math_force_eval(u); */
  42                 return x;
  43         }
  44
  45         if (hx >= 0) { /* x > 0 */
  46                 if (hx > hy) { /* x > y: x -= ulp */
  47                         hx -= 1;
  48                 } else { /* x < y: x += ulp */
  49                         hx += 1;
  50                 }
  51         } else { /* x < 0 */
  52                 if (hy >= 0 || hx > hy) { /* x < y: x -= ulp */
  53                         hx -= 1;
  54                 } else { /* x > y: x += ulp */
  55                         hx += 1;
  56                 }
  57         }
  58         hy = hx & 0x7f800000;
  59         if (hy >= 0x7f800000) {
  60                 x = x + x; /* overflow */
  61                 return x; /* overflow */
  62         }
  63         if (hy < 0x00800000) {
  64                 float u = x * x; /* underflow */
  65                 math_force_eval(u); /* raise underflow flag */
  66         }
  67         SET_FLOAT_WORD(x, hx);
  68         return x;
  69 }
  70
  71 #if 0
  72 /* "testprog N a b"
  73  * calculates a = nextafterf(a, b) and prints a as float
  74  * and as raw bytes; repeats it N times.
  75  */
  76 #include <stdio.h>
  77 #include <stdlib.h>
  78 #include <math.h>
  79 int main(int argc, char **argv)
  80 {
  81         int cnt, i;
  82         float a, b;
  83         cnt = atoi(argv[1]);
  84         a = strtod(argv[2], NULL);
  85         b = strtod(argv[3], NULL);
  86         while (cnt-- > 0) {
  87                 for (i = 0; i < sizeof(a); i++) {
  88                         unsigned char c = ((char*)(&a))[i];
  89                         printf("%x%x", (c >> 4), (c & 0xf));
  90                 }
  91                 printf(" %f\n", a);
  92                 a = nextafterf(a, b);
  93         }
  94         return 0;
  95 }
  96 #endif