usr/src/lib/libm/common/m9x/nexttoward.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak nexttoward = __nexttoward
  31
  32 /*
  33  * nexttoward(x, y) delivers the next representable number after x
  34  * in the direction of y.  If x and y are both zero, the result is
  35  * zero with the same sign as y.  If either x or y is NaN, the result
  36  * is NaN.
  37  *
  38  * If x != y and the result is infinite, overflow is raised; if
  39  * x != y and the result is subnormal or zero, underflow is raised.
  40  * (This is wrong, but it's what C99 apparently wants.)
  41  */
  42
  43 #include "libm.h"
  44
  45 #if defined(__sparc)
  46
  47 static union {
  48         unsigned i[2];
  49         double d;
  50 } C[] = {
  51         0x00100000, 0,
  52         0x7fe00000, 0,
  53         0x7fffffff, 0xffffffff
  54 };
  55
  56 #define tiny    C[0].d
  57 #define huge    C[1].d
  58 #define qnan    C[2].d
  59
  60 enum fcc_type {
  61         fcc_equal = 0,
  62         fcc_less = 1,
  63         fcc_greater = 2,
  64         fcc_unordered = 3
  65 };
  66
  67 #ifdef __sparcv9
  68 #define _Q_cmp  _Qp_cmp
  69 #endif
  70
  71 extern enum fcc_type _Q_cmp(const long double *, const long double *);
  72
  73 double
  74 __nexttoward(double x, long double y) {
  75         union {
  76                 unsigned i[2];
  77                 double d;
  78         } xx;
  79         union {
  80                 unsigned i[4];
  81                 long double q;
  82         } yy;
  83         long double lx;
  84         unsigned hx;
  85         volatile double dummy;
  86         enum fcc_type rel;
  87
  88         /*
  89          * It would be somewhat more efficient to check for NaN and
  90          * zero operands before converting x to long double and then
  91          * to code the comparison in line rather than calling _Q_cmp.
  92          * However, since this code probably won't get used much,
  93          * I'm opting in favor of simplicity instead.
  94          */
  95         lx = xx.d = x;
  96         hx = (xx.i[0] & ~0x80000000) | xx.i[1];
  97
  98         /* check for each of four possible orderings */
  99         rel = _Q_cmp(&lx, &y);
 100         if (rel == fcc_unordered)
 101                 return (qnan);
 102
 103         if (rel == fcc_equal) {
 104                 if (hx == 0) {  /* x is zero; return zero with y's sign */
 105                         yy.q = y;
 106                         xx.i[0] = yy.i[0];
 107                         return (xx.d);
 108                 }
 109                 return (x);
 110         }
 111
 112         if (rel == fcc_less) {
 113                 if (hx == 0) {  /* x is zero */
 114                         xx.i[0] = 0;
 115                         xx.i[1] = 0x00000001;
 116                 } else if ((int)xx.i[0] >= 0) { /* x is positive */
 117                         if (++xx.i[1] == 0)
 118                                 xx.i[0]++;
 119                 } else {
 120                         if (xx.i[1]-- == 0)
 121                                 xx.i[0]--;
 122                 }
 123         } else {
 124                 if (hx == 0) {  /* x is zero */
 125                         xx.i[0] = 0x80000000;
 126                         xx.i[1] = 0x00000001;
 127                 } else if ((int)xx.i[0] >= 0) { /* x is positive */
 128                         if (xx.i[1]-- == 0)
 129                                 xx.i[0]--;
 130                 } else {
 131                         if (++xx.i[1] == 0)
 132                                 xx.i[0]++;
 133                 }
 134         }
 135
 136         /* raise exceptions as needed */
 137         hx = xx.i[0] & ~0x80000000;
 138         if (hx == 0x7ff00000) {
 139                 dummy = huge;
 140                 dummy *= huge;
 141         } else if (hx < 0x00100000) {
 142                 dummy = tiny;
 143                 dummy *= tiny;
 144         }
 145
 146         return (xx.d);
 147 }
 148
 149 #elif defined(__x86)
 150
 151 static union {
 152         unsigned i[2];
 153         double d;
 154 } C[] = {
 155         0, 0x00100000,
 156         0, 0x7fe00000,
 157 };
 158
 159 #define tiny    C[0].d
 160 #define huge    C[1].d
 161
 162 double
 163 __nexttoward(double x, long double y) {
 164         union {
 165                 unsigned i[2];
 166                 double d;
 167         } xx;
 168         unsigned hx;
 169         long double lx;
 170         volatile double dummy;
 171
 172         lx = xx.d = x;
 173         hx = (xx.i[1] & ~0x80000000) | xx.i[0];
 174
 175         /* check for each of four possible orderings */
 176         if (isunordered(lx, y))
 177                 return ((double) (lx + y));
 178
 179         if (lx == y)
 180                 return ((double) y);
 181
 182         if (lx < y) {
 183                 if (hx == 0) {  /* x is zero */
 184                         xx.i[0] = 0x00000001;
 185                         xx.i[1] = 0;
 186                 } else if ((int)xx.i[1] >= 0) { /* x is positive */
 187                         if (++xx.i[0] == 0)
 188                                 xx.i[1]++;
 189                 } else {
 190                         if (xx.i[0]-- == 0)
 191                                 xx.i[1]--;
 192                 }
 193         } else {
 194                 if (hx == 0) {  /* x is zero */
 195                         xx.i[0] = 0x00000001;
 196                         xx.i[1] = 0x80000000;
 197                 } else if ((int)xx.i[1] >= 0) { /* x is positive */
 198                         if (xx.i[0]-- == 0)
 199                                 xx.i[1]--;
 200                 } else {
 201                         if (++xx.i[0] == 0)
 202                                 xx.i[1]++;
 203                 }
 204         }
 205
 206         /* raise exceptions as needed */
 207         hx = xx.i[1] & ~0x80000000;
 208         if (hx == 0x7ff00000) {
 209                 dummy = huge;
 210                 dummy *= huge;
 211         } else if (hx < 0x00100000) {
 212                 dummy = tiny;
 213                 dummy *= tiny;
 214         }
 215
 216         return (xx.d);
 217 }
 218
 219 #else
 220 #error Unknown architecture
 221 #endif