Objects/floatobject.c

   1
   2 /* Float object implementation */
   3
   4 /* XXX There should be overflow checks here, but it's hard to check
   5    for any kind of float exception without losing portability. */
   6
   7 #include "Python.h"
   8 #include "structseq.h"
   9
  10 #include <ctype.h>
  11 #include <float.h>
  12
  13 #undef MAX
  14 #undef MIN
  15 #define MAX(x, y) ((x) < (y) ? (y) : (x))
  16 #define MIN(x, y) ((x) < (y) ? (x) : (y))
  17
  18 #ifdef HAVE_IEEEFP_H
  19 #include <ieeefp.h>
  20 #endif
  21
  22
  23 #ifdef _OSF_SOURCE
  24 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
  25 extern int finite(double);
  26 #endif
  27
  28 /* Special free list -- see comments for same code in intobject.c. */
  29 #define BLOCK_SIZE      1000    /* 1K less typical malloc overhead */
  30 #define BHEAD_SIZE      8       /* Enough for a 64-bit pointer */
  31 #define N_FLOATOBJECTS  ((BLOCK_SIZE - BHEAD_SIZE) / sizeof(PyFloatObject))
  32
  33 struct _floatblock {
  34     struct _floatblock *next;
  35     PyFloatObject objects[N_FLOATOBJECTS];
  36 };
  37
  38 typedef struct _floatblock PyFloatBlock;
  39
  40 static PyFloatBlock *block_list = NULL;
  41 static PyFloatObject *free_list = NULL;
  42
  43 static PyFloatObject *
  44 fill_free_list(void)
  45 {
  46     PyFloatObject *p, *q;
  47     /* XXX Float blocks escape the object heap. Use PyObject_MALLOC ??? */
  48     p = (PyFloatObject *) PyMem_MALLOC(sizeof(PyFloatBlock));
  49     if (p == NULL)
  50         return (PyFloatObject *) PyErr_NoMemory();
  51     ((PyFloatBlock *)p)->next = block_list;
  52     block_list = (PyFloatBlock *)p;
  53     p = &((PyFloatBlock *)p)->objects[0];
  54     q = p + N_FLOATOBJECTS;
  55     while (--q > p)
  56         Py_TYPE(q) = (struct _typeobject *)(q-1);
  57     Py_TYPE(q) = NULL;
  58     return p + N_FLOATOBJECTS - 1;
  59 }
  60
  61 double
  62 PyFloat_GetMax(void)
  63 {
  64     return DBL_MAX;
  65 }
  66
  67 double
  68 PyFloat_GetMin(void)
  69 {
  70     return DBL_MIN;
  71 }
  72
  73 static PyTypeObject FloatInfoType;
  74
  75 PyDoc_STRVAR(floatinfo__doc__,
  76 "sys.floatinfo\n\
  77 \n\
  78 A structseq holding information about the float type. It contains low level\n\
  79 information about the precision and internal representation. Please study\n\
  80 your system's :file:`float.h` for more information.");
  81
  82 static PyStructSequence_Field floatinfo_fields[] = {
  83     {"max",             "DBL_MAX -- maximum representable finite float"},
  84     {"max_exp",         "DBL_MAX_EXP -- maximum int e such that radix**(e-1) "
  85                     "is representable"},
  86     {"max_10_exp",      "DBL_MAX_10_EXP -- maximum int e such that 10**e "
  87                     "is representable"},
  88     {"min",             "DBL_MIN -- Minimum positive normalizer float"},
  89     {"min_exp",         "DBL_MIN_EXP -- minimum int e such that radix**(e-1) "
  90                     "is a normalized float"},
  91     {"min_10_exp",      "DBL_MIN_10_EXP -- minimum int e such that 10**e is "
  92                     "a normalized"},
  93     {"dig",             "DBL_DIG -- digits"},
  94     {"mant_dig",        "DBL_MANT_DIG -- mantissa digits"},
  95     {"epsilon",         "DBL_EPSILON -- Difference between 1 and the next "
  96                     "representable float"},
  97     {"radix",           "FLT_RADIX -- radix of exponent"},
  98     {"rounds",          "FLT_ROUNDS -- addition rounds"},
  99     {0}
 100 };
 101
 102 static PyStructSequence_Desc floatinfo_desc = {
 103     "sys.floatinfo",            /* name */
 104     floatinfo__doc__,           /* doc */
 105     floatinfo_fields,           /* fields */
 106     11
 107 };
 108
 109 PyObject *
 110 PyFloat_GetInfo(void)
 111 {
 112     PyObject* floatinfo;
 113     int pos = 0;
 114
 115     floatinfo = PyStructSequence_New(&FloatInfoType);
 116     if (floatinfo == NULL) {
 117         return NULL;
 118     }
 119
 120 #define SetIntFlag(flag) \
 121     PyStructSequence_SET_ITEM(floatinfo, pos++, PyLong_FromLong(flag))
 122 #define SetDblFlag(flag) \
 123     PyStructSequence_SET_ITEM(floatinfo, pos++, PyFloat_FromDouble(flag))
 124
 125     SetDblFlag(DBL_MAX);
 126     SetIntFlag(DBL_MAX_EXP);
 127     SetIntFlag(DBL_MAX_10_EXP);
 128     SetDblFlag(DBL_MIN);
 129     SetIntFlag(DBL_MIN_EXP);
 130     SetIntFlag(DBL_MIN_10_EXP);
 131     SetIntFlag(DBL_DIG);
 132     SetIntFlag(DBL_MANT_DIG);
 133     SetDblFlag(DBL_EPSILON);
 134     SetIntFlag(FLT_RADIX);
 135     SetIntFlag(FLT_ROUNDS);
 136 #undef SetIntFlag
 137 #undef SetDblFlag
 138
 139     if (PyErr_Occurred()) {
 140         Py_CLEAR(floatinfo);
 141         return NULL;
 142     }
 143     return floatinfo;
 144 }
 145
 146 PyObject *
 147 PyFloat_FromDouble(double fval)
 148 {
 149     register PyFloatObject *op;
 150     if (free_list == NULL) {
 151         if ((free_list = fill_free_list()) == NULL)
 152             return NULL;
 153     }
 154     /* Inline PyObject_New */
 155     op = free_list;
 156     free_list = (PyFloatObject *)Py_TYPE(op);
 157     PyObject_INIT(op, &PyFloat_Type);
 158     op->ob_fval = fval;
 159     return (PyObject *) op;
 160 }
 161
 162 PyObject *
 163 PyFloat_FromString(PyObject *v)
 164 {
 165     const char *s, *last, *end;
 166     double x;
 167     char buffer[256]; /* for errors */
 168     char *s_buffer = NULL;
 169     Py_ssize_t len;
 170     PyObject *result = NULL;
 171
 172     if (PyUnicode_Check(v)) {
 173         s_buffer = (char *)PyMem_MALLOC(PyUnicode_GET_SIZE(v)+1);
 174         if (s_buffer == NULL)
 175             return PyErr_NoMemory();
 176         if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
 177                                     PyUnicode_GET_SIZE(v),
 178                                     s_buffer,
 179                                     NULL))
 180             goto error;
 181         s = s_buffer;
 182         len = strlen(s);
 183     }
 184     else if (PyObject_AsCharBuffer(v, &s, &len)) {
 185         PyErr_SetString(PyExc_TypeError,
 186                         "float() argument must be a string or a number");
 187         return NULL;
 188     }
 189     last = s + len;
 190
 191     while (Py_ISSPACE(*s))
 192         s++;
 193     /* We don't care about overflow or underflow.  If the platform
 194      * supports them, infinities and signed zeroes (on underflow) are
 195      * fine. */
 196     x = PyOS_string_to_double(s, (char **)&end, NULL);
 197     if (x == -1.0 && PyErr_Occurred())
 198         goto error;
 199     while (Py_ISSPACE(*end))
 200         end++;
 201     if (end == last)
 202         result = PyFloat_FromDouble(x);
 203     else {
 204         PyOS_snprintf(buffer, sizeof(buffer),
 205                       "invalid literal for float(): %.200s", s);
 206         PyErr_SetString(PyExc_ValueError, buffer);
 207         result = NULL;
 208     }
 209
 210   error:
 211     if (s_buffer)
 212         PyMem_FREE(s_buffer);
 213     return result;
 214 }
 215
 216 static void
 217 float_dealloc(PyFloatObject *op)
 218 {
 219     if (PyFloat_CheckExact(op)) {
 220         Py_TYPE(op) = (struct _typeobject *)free_list;
 221         free_list = op;
 222     }
 223     else
 224         Py_TYPE(op)->tp_free((PyObject *)op);
 225 }
 226
 227 double
 228 PyFloat_AsDouble(PyObject *op)
 229 {
 230     PyNumberMethods *nb;
 231     PyFloatObject *fo;
 232     double val;
 233
 234     if (op && PyFloat_Check(op))
 235         return PyFloat_AS_DOUBLE((PyFloatObject*) op);
 236
 237     if (op == NULL) {
 238         PyErr_BadArgument();
 239         return -1;
 240     }
 241
 242     if ((nb = Py_TYPE(op)->tp_as_number) == NULL || nb->nb_float == NULL) {
 243         PyErr_SetString(PyExc_TypeError, "a float is required");
 244         return -1;
 245     }
 246
 247     fo = (PyFloatObject*) (*nb->nb_float) (op);
 248     if (fo == NULL)
 249         return -1;
 250     if (!PyFloat_Check(fo)) {
 251         PyErr_SetString(PyExc_TypeError,
 252                         "nb_float should return float object");
 253         return -1;
 254     }
 255
 256     val = PyFloat_AS_DOUBLE(fo);
 257     Py_DECREF(fo);
 258
 259     return val;
 260 }
 261
 262 /* Macro and helper that convert PyObject obj to a C double and store
 263    the value in dbl.  If conversion to double raises an exception, obj is
 264    set to NULL, and the function invoking this macro returns NULL.  If
 265    obj is not of float, int or long type, Py_NotImplemented is incref'ed,
 266    stored in obj, and returned from the function invoking this macro.
 267 */
 268 #define CONVERT_TO_DOUBLE(obj, dbl)                     \
 269     if (PyFloat_Check(obj))                             \
 270         dbl = PyFloat_AS_DOUBLE(obj);                   \
 271     else if (convert_to_double(&(obj), &(dbl)) < 0)     \
 272         return obj;
 273
 274 /* Methods */
 275
 276 static int
 277 convert_to_double(PyObject **v, double *dbl)
 278 {
 279     register PyObject *obj = *v;
 280
 281     if (PyLong_Check(obj)) {
 282         *dbl = PyLong_AsDouble(obj);
 283         if (*dbl == -1.0 && PyErr_Occurred()) {
 284             *v = NULL;
 285             return -1;
 286         }
 287     }
 288     else {
 289         Py_INCREF(Py_NotImplemented);
 290         *v = Py_NotImplemented;
 291         return -1;
 292     }
 293     return 0;
 294 }
 295
 296 static PyObject *
 297 float_str_or_repr(PyFloatObject *v, int precision, char format_code)
 298 {
 299     PyObject *result;
 300     char *buf = PyOS_double_to_string(PyFloat_AS_DOUBLE(v),
 301                                       format_code, precision,
 302                                       Py_DTSF_ADD_DOT_0,
 303                                       NULL);
 304     if (!buf)
 305     return PyErr_NoMemory();
 306     result = PyUnicode_FromString(buf);
 307     PyMem_Free(buf);
 308     return result;
 309 }
 310
 311 static PyObject *
 312 float_repr(PyFloatObject *v)
 313 {
 314     return float_str_or_repr(v, 0, 'r');
 315 }
 316
 317 static PyObject *
 318 float_str(PyFloatObject *v)
 319 {
 320     return float_str_or_repr(v, PyFloat_STR_PRECISION, 'g');
 321 }
 322
 323 /* Comparison is pretty much a nightmare.  When comparing float to float,
 324  * we do it as straightforwardly (and long-windedly) as conceivable, so
 325  * that, e.g., Python x == y delivers the same result as the platform
 326  * C x == y when x and/or y is a NaN.
 327  * When mixing float with an integer type, there's no good *uniform* approach.
 328  * Converting the double to an integer obviously doesn't work, since we
 329  * may lose info from fractional bits.  Converting the integer to a double
 330  * also has two failure modes:  (1) a long int may trigger overflow (too
 331  * large to fit in the dynamic range of a C double); (2) even a C long may have
 332  * more bits than fit in a C double (e.g., on a a 64-bit box long may have
 333  * 63 bits of precision, but a C double probably has only 53), and then
 334  * we can falsely claim equality when low-order integer bits are lost by
 335  * coercion to double.  So this part is painful too.
 336  */
 337
 338 static PyObject*
 339 float_richcompare(PyObject *v, PyObject *w, int op)
 340 {
 341     double i, j;
 342     int r = 0;
 343
 344     assert(PyFloat_Check(v));
 345     i = PyFloat_AS_DOUBLE(v);
 346
 347     /* Switch on the type of w.  Set i and j to doubles to be compared,
 348      * and op to the richcomp to use.
 349      */
 350     if (PyFloat_Check(w))
 351         j = PyFloat_AS_DOUBLE(w);
 352
 353     else if (!Py_IS_FINITE(i)) {
 354         if (PyLong_Check(w))
 355             /* If i is an infinity, its magnitude exceeds any
 356              * finite integer, so it doesn't matter which int we
 357              * compare i with.  If i is a NaN, similarly.
 358              */
 359             j = 0.0;
 360         else
 361             goto Unimplemented;
 362     }
 363
 364     else if (PyLong_Check(w)) {
 365         int vsign = i == 0.0 ? 0 : i < 0.0 ? -1 : 1;
 366         int wsign = _PyLong_Sign(w);
 367         size_t nbits;
 368         int exponent;
 369
 370         if (vsign != wsign) {
 371             /* Magnitudes are irrelevant -- the signs alone
 372              * determine the outcome.
 373              */
 374             i = (double)vsign;
 375             j = (double)wsign;
 376             goto Compare;
 377         }
 378         /* The signs are the same. */
 379         /* Convert w to a double if it fits.  In particular, 0 fits. */
 380         nbits = _PyLong_NumBits(w);
 381         if (nbits == (size_t)-1 && PyErr_Occurred()) {
 382             /* This long is so large that size_t isn't big enough
 383              * to hold the # of bits.  Replace with little doubles
 384              * that give the same outcome -- w is so large that
 385              * its magnitude must exceed the magnitude of any
 386              * finite float.
 387              */
 388             PyErr_Clear();
 389             i = (double)vsign;
 390             assert(wsign != 0);
 391             j = wsign * 2.0;
 392             goto Compare;
 393         }
 394         if (nbits <= 48) {
 395             j = PyLong_AsDouble(w);
 396             /* It's impossible that <= 48 bits overflowed. */
 397             assert(j != -1.0 || ! PyErr_Occurred());
 398             goto Compare;
 399         }
 400         assert(wsign != 0); /* else nbits was 0 */
 401         assert(vsign != 0); /* if vsign were 0, then since wsign is
 402                              * not 0, we would have taken the
 403                              * vsign != wsign branch at the start */
 404         /* We want to work with non-negative numbers. */
 405         if (vsign < 0) {
 406             /* "Multiply both sides" by -1; this also swaps the
 407              * comparator.
 408              */
 409             i = -i;
 410             op = _Py_SwappedOp[op];
 411         }
 412         assert(i > 0.0);
 413         (void) frexp(i, &exponent);
 414         /* exponent is the # of bits in v before the radix point;
 415          * we know that nbits (the # of bits in w) > 48 at this point
 416          */
 417         if (exponent < 0 || (size_t)exponent < nbits) {
 418             i = 1.0;
 419             j = 2.0;
 420             goto Compare;
 421         }
 422         if ((size_t)exponent > nbits) {
 423             i = 2.0;
 424             j = 1.0;
 425             goto Compare;
 426         }
 427         /* v and w have the same number of bits before the radix
 428          * point.  Construct two longs that have the same comparison
 429          * outcome.
 430          */
 431         {
 432             double fracpart;
 433             double intpart;
 434             PyObject *result = NULL;
 435             PyObject *one = NULL;
 436             PyObject *vv = NULL;
 437             PyObject *ww = w;
 438
 439             if (wsign < 0) {
 440                 ww = PyNumber_Negative(w);
 441                 if (ww == NULL)
 442                     goto Error;
 443             }
 444             else
 445                 Py_INCREF(ww);
 446
 447             fracpart = modf(i, &intpart);
 448             vv = PyLong_FromDouble(intpart);
 449             if (vv == NULL)
 450                 goto Error;
 451
 452             if (fracpart != 0.0) {
 453                 /* Shift left, and or a 1 bit into vv
 454                  * to represent the lost fraction.
 455                  */
 456                 PyObject *temp;
 457
 458                 one = PyLong_FromLong(1);
 459                 if (one == NULL)
 460                     goto Error;
 461
 462                 temp = PyNumber_Lshift(ww, one);
 463                 if (temp == NULL)
 464                     goto Error;
 465                 Py_DECREF(ww);
 466                 ww = temp;
 467
 468                 temp = PyNumber_Lshift(vv, one);
 469                 if (temp == NULL)
 470                     goto Error;
 471                 Py_DECREF(vv);
 472                 vv = temp;
 473
 474                 temp = PyNumber_Or(vv, one);
 475                 if (temp == NULL)
 476                     goto Error;
 477                 Py_DECREF(vv);
 478                 vv = temp;
 479             }
 480
 481             r = PyObject_RichCompareBool(vv, ww, op);
 482             if (r < 0)
 483                 goto Error;
 484             result = PyBool_FromLong(r);
 485          Error:
 486             Py_XDECREF(vv);
 487             Py_XDECREF(ww);
 488             Py_XDECREF(one);
 489             return result;
 490         }
 491     } /* else if (PyLong_Check(w)) */
 492
 493     else        /* w isn't float, int, or long */
 494         goto Unimplemented;
 495
 496  Compare:
 497     PyFPE_START_PROTECT("richcompare", return NULL)
 498     switch (op) {
 499     case Py_EQ:
 500         r = i == j;
 501         break;
 502     case Py_NE:
 503         r = i != j;
 504         break;
 505     case Py_LE:
 506         r = i <= j;
 507         break;
 508     case Py_GE:
 509         r = i >= j;
 510         break;
 511     case Py_LT:
 512         r = i < j;
 513         break;
 514     case Py_GT:
 515         r = i > j;
 516         break;
 517     }
 518     PyFPE_END_PROTECT(r)
 519     return PyBool_FromLong(r);
 520
 521  Unimplemented:
 522     Py_INCREF(Py_NotImplemented);
 523     return Py_NotImplemented;
 524 }
 525
 526 static long
 527 float_hash(PyFloatObject *v)
 528 {
 529     return _Py_HashDouble(v->ob_fval);
 530 }
 531
 532 static PyObject *
 533 float_add(PyObject *v, PyObject *w)
 534 {
 535     double a,b;
 536     CONVERT_TO_DOUBLE(v, a);
 537     CONVERT_TO_DOUBLE(w, b);
 538     PyFPE_START_PROTECT("add", return 0)
 539     a = a + b;
 540     PyFPE_END_PROTECT(a)
 541     return PyFloat_FromDouble(a);
 542 }
 543
 544 static PyObject *
 545 float_sub(PyObject *v, PyObject *w)
 546 {
 547     double a,b;
 548     CONVERT_TO_DOUBLE(v, a);
 549     CONVERT_TO_DOUBLE(w, b);
 550     PyFPE_START_PROTECT("subtract", return 0)
 551     a = a - b;
 552     PyFPE_END_PROTECT(a)
 553     return PyFloat_FromDouble(a);
 554 }
 555
 556 static PyObject *
 557 float_mul(PyObject *v, PyObject *w)
 558 {
 559     double a,b;
 560     CONVERT_TO_DOUBLE(v, a);
 561     CONVERT_TO_DOUBLE(w, b);
 562     PyFPE_START_PROTECT("multiply", return 0)
 563     a = a * b;
 564     PyFPE_END_PROTECT(a)
 565     return PyFloat_FromDouble(a);
 566 }
 567
 568 static PyObject *
 569 float_div(PyObject *v, PyObject *w)
 570 {
 571     double a,b;
 572     CONVERT_TO_DOUBLE(v, a);
 573     CONVERT_TO_DOUBLE(w, b);
 574 #ifdef Py_NAN
 575     if (b == 0.0) {
 576         PyErr_SetString(PyExc_ZeroDivisionError,
 577                         "float division");
 578         return NULL;
 579     }
 580 #endif
 581     PyFPE_START_PROTECT("divide", return 0)
 582     a = a / b;
 583     PyFPE_END_PROTECT(a)
 584     return PyFloat_FromDouble(a);
 585 }
 586
 587 static PyObject *
 588 float_rem(PyObject *v, PyObject *w)
 589 {
 590     double vx, wx;
 591     double mod;
 592     CONVERT_TO_DOUBLE(v, vx);
 593     CONVERT_TO_DOUBLE(w, wx);
 594 #ifdef Py_NAN
 595     if (wx == 0.0) {
 596         PyErr_SetString(PyExc_ZeroDivisionError,
 597                         "float modulo");
 598         return NULL;
 599     }
 600 #endif
 601     PyFPE_START_PROTECT("modulo", return 0)
 602     mod = fmod(vx, wx);
 603     /* note: checking mod*wx < 0 is incorrect -- underflows to
 604        0 if wx < sqrt(smallest nonzero double) */
 605     if (mod && ((wx < 0) != (mod < 0))) {
 606         mod += wx;
 607     }
 608     PyFPE_END_PROTECT(mod)
 609     return PyFloat_FromDouble(mod);
 610 }
 611
 612 static PyObject *
 613 float_divmod(PyObject *v, PyObject *w)
 614 {
 615     double vx, wx;
 616     double div, mod, floordiv;
 617     CONVERT_TO_DOUBLE(v, vx);
 618     CONVERT_TO_DOUBLE(w, wx);
 619     if (wx == 0.0) {
 620         PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");
 621         return NULL;
 622     }
 623     PyFPE_START_PROTECT("divmod", return 0)
 624     mod = fmod(vx, wx);
 625     /* fmod is typically exact, so vx-mod is *mathematically* an
 626        exact multiple of wx.  But this is fp arithmetic, and fp
 627        vx - mod is an approximation; the result is that div may
 628        not be an exact integral value after the division, although
 629        it will always be very close to one.
 630     */
 631     div = (vx - mod) / wx;
 632     if (mod) {
 633         /* ensure the remainder has the same sign as the denominator */
 634         if ((wx < 0) != (mod < 0)) {
 635             mod += wx;
 636             div -= 1.0;
 637         }
 638     }
 639     else {
 640         /* the remainder is zero, and in the presence of signed zeroes
 641            fmod returns different results across platforms; ensure
 642            it has the same sign as the denominator; we'd like to do
 643            "mod = wx * 0.0", but that may get optimized away */
 644         mod *= mod;  /* hide "mod = +0" from optimizer */
 645         if (wx < 0.0)
 646             mod = -mod;
 647     }
 648     /* snap quotient to nearest integral value */
 649     if (div) {
 650         floordiv = floor(div);
 651         if (div - floordiv > 0.5)
 652             floordiv += 1.0;
 653     }
 654     else {
 655         /* div is zero - get the same sign as the true quotient */
 656         div *= div;             /* hide "div = +0" from optimizers */
 657         floordiv = div * vx / wx; /* zero w/ sign of vx/wx */
 658     }
 659     PyFPE_END_PROTECT(floordiv)
 660     return Py_BuildValue("(dd)", floordiv, mod);
 661 }
 662
 663 static PyObject *
 664 float_floor_div(PyObject *v, PyObject *w)
 665 {
 666     PyObject *t, *r;
 667
 668     t = float_divmod(v, w);
 669     if (t == NULL || t == Py_NotImplemented)
 670         return t;
 671     assert(PyTuple_CheckExact(t));
 672     r = PyTuple_GET_ITEM(t, 0);
 673     Py_INCREF(r);
 674     Py_DECREF(t);
 675     return r;
 676 }
 677
 678 static PyObject *
 679 float_pow(PyObject *v, PyObject *w, PyObject *z)
 680 {
 681     double iv, iw, ix;
 682
 683     if ((PyObject *)z != Py_None) {
 684         PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "
 685             "allowed unless all arguments are integers");
 686         return NULL;
 687     }
 688
 689     CONVERT_TO_DOUBLE(v, iv);
 690     CONVERT_TO_DOUBLE(w, iw);
 691
 692     /* Sort out special cases here instead of relying on pow() */
 693     if (iw == 0) {              /* v**0 is 1, even 0**0 */
 694         return PyFloat_FromDouble(1.0);
 695     }
 696     if (iv == 0.0) {  /* 0**w is error if w<0, else 1 */
 697         if (iw < 0.0) {
 698             PyErr_SetString(PyExc_ZeroDivisionError,
 699                             "0.0 cannot be raised to a negative power");
 700             return NULL;
 701         }
 702         return PyFloat_FromDouble(0.0);
 703     }
 704     if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
 705         return PyFloat_FromDouble(1.0);
 706     }
 707     if (iv < 0.0) {
 708         /* Whether this is an error is a mess, and bumps into libm
 709          * bugs so we have to figure it out ourselves.
 710          */
 711         if (iw != floor(iw)) {
 712             /* Negative numbers raised to fractional powers
 713              * become complex.
 714              */
 715             return PyComplex_Type.tp_as_number->nb_power(v, w, z);
 716         }
 717         /* iw is an exact integer, albeit perhaps a very large one.
 718          * -1 raised to an exact integer should never be exceptional.
 719          * Alas, some libms (chiefly glibc as of early 2003) return
 720          * NaN and set EDOM on pow(-1, large_int) if the int doesn't
 721          * happen to be representable in a *C* integer.  That's a
 722          * bug; we let that slide in math.pow() (which currently
 723          * reflects all platform accidents), but not for Python's **.
 724          */
 725          if (iv == -1.0 && Py_IS_FINITE(iw)) {
 726             /* Return 1 if iw is even, -1 if iw is odd; there's
 727              * no guarantee that any C integral type is big
 728              * enough to hold iw, so we have to check this
 729              * indirectly.
 730              */
 731             ix = floor(iw * 0.5) * 2.0;
 732             return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0);
 733         }
 734         /* Else iv != -1.0, and overflow or underflow are possible.
 735          * Unless we're to write pow() ourselves, we have to trust
 736          * the platform to do this correctly.
 737          */
 738     }
 739     errno = 0;
 740     PyFPE_START_PROTECT("pow", return NULL)
 741     ix = pow(iv, iw);
 742     PyFPE_END_PROTECT(ix)
 743     Py_ADJUST_ERANGE1(ix);
 744     if (errno != 0) {
 745         /* We don't expect any errno value other than ERANGE, but
 746          * the range of libm bugs appears unbounded.
 747          */
 748         PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 749                              PyExc_ValueError);
 750         return NULL;
 751     }
 752     return PyFloat_FromDouble(ix);
 753 }
 754
 755 static PyObject *
 756 float_neg(PyFloatObject *v)
 757 {
 758     return PyFloat_FromDouble(-v->ob_fval);
 759 }
 760
 761 static PyObject *
 762 float_abs(PyFloatObject *v)
 763 {
 764     return PyFloat_FromDouble(fabs(v->ob_fval));
 765 }
 766
 767 static int
 768 float_bool(PyFloatObject *v)
 769 {
 770     return v->ob_fval != 0.0;
 771 }
 772
 773 static PyObject *
 774 float_is_integer(PyObject *v)
 775 {
 776     double x = PyFloat_AsDouble(v);
 777     PyObject *o;
 778
 779     if (x == -1.0 && PyErr_Occurred())
 780         return NULL;
 781     if (!Py_IS_FINITE(x))
 782         Py_RETURN_FALSE;
 783     errno = 0;
 784     PyFPE_START_PROTECT("is_integer", return NULL)
 785     o = (floor(x) == x) ? Py_True : Py_False;
 786     PyFPE_END_PROTECT(x)
 787     if (errno != 0) {
 788         PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 789                              PyExc_ValueError);
 790         return NULL;
 791     }
 792     Py_INCREF(o);
 793     return o;
 794 }
 795
 796 #if 0
 797 static PyObject *
 798 float_is_inf(PyObject *v)
 799 {
 800     double x = PyFloat_AsDouble(v);
 801     if (x == -1.0 && PyErr_Occurred())
 802         return NULL;
 803     return PyBool_FromLong((long)Py_IS_INFINITY(x));
 804 }
 805
 806 static PyObject *
 807 float_is_nan(PyObject *v)
 808 {
 809     double x = PyFloat_AsDouble(v);
 810     if (x == -1.0 && PyErr_Occurred())
 811         return NULL;
 812     return PyBool_FromLong((long)Py_IS_NAN(x));
 813 }
 814
 815 static PyObject *
 816 float_is_finite(PyObject *v)
 817 {
 818     double x = PyFloat_AsDouble(v);
 819     if (x == -1.0 && PyErr_Occurred())
 820         return NULL;
 821     return PyBool_FromLong((long)Py_IS_FINITE(x));
 822 }
 823 #endif
 824
 825 static PyObject *
 826 float_trunc(PyObject *v)
 827 {
 828     double x = PyFloat_AsDouble(v);
 829     double wholepart;           /* integral portion of x, rounded toward 0 */
 830
 831     (void)modf(x, &wholepart);
 832     /* Try to get out cheap if this fits in a Python int.  The attempt
 833      * to cast to long must be protected, as C doesn't define what
 834      * happens if the double is too big to fit in a long.  Some rare
 835      * systems raise an exception then (RISCOS was mentioned as one,
 836      * and someone using a non-default option on Sun also bumped into
 837      * that).  Note that checking for >= and <= LONG_{MIN,MAX} would
 838      * still be vulnerable:  if a long has more bits of precision than
 839      * a double, casting MIN/MAX to double may yield an approximation,
 840      * and if that's rounded up, then, e.g., wholepart=LONG_MAX+1 would
 841      * yield true from the C expression wholepart<=LONG_MAX, despite
 842      * that wholepart is actually greater than LONG_MAX.
 843      */
 844     if (LONG_MIN < wholepart && wholepart < LONG_MAX) {
 845         const long aslong = (long)wholepart;
 846         return PyLong_FromLong(aslong);
 847     }
 848     return PyLong_FromDouble(wholepart);
 849 }
 850
 851 /* double_round: rounds a finite double to the closest multiple of
 852    10**-ndigits; here ndigits is within reasonable bounds (typically, -308 <=
 853    ndigits <= 323).  Returns a Python float, or sets a Python error and
 854    returns NULL on failure (OverflowError and memory errors are possible). */
 855
 856 #ifndef PY_NO_SHORT_FLOAT_REPR
 857 /* version of double_round that uses the correctly-rounded string<->double
 858    conversions from Python/dtoa.c */
 859
 860 static PyObject *
 861 double_round(double x, int ndigits) {
 862
 863     double rounded;
 864     Py_ssize_t buflen, mybuflen=100;
 865     char *buf, *buf_end, shortbuf[100], *mybuf=shortbuf;
 866     int decpt, sign;
 867     PyObject *result = NULL;
 868
 869     /* round to a decimal string */
 870     buf = _Py_dg_dtoa(x, 3, ndigits, &decpt, &sign, &buf_end);
 871     if (buf == NULL) {
 872         PyErr_NoMemory();
 873         return NULL;
 874     }
 875
 876     /* Get new buffer if shortbuf is too small.  Space needed <= buf_end -
 877     buf + 8: (1 extra for '0', 1 for sign, 5 for exp, 1 for '\0').  */
 878     buflen = buf_end - buf;
 879     if (buflen + 8 > mybuflen) {
 880         mybuflen = buflen+8;
 881         mybuf = (char *)PyMem_Malloc(mybuflen);
 882         if (mybuf == NULL) {
 883             PyErr_NoMemory();
 884             goto exit;
 885         }
 886     }
 887     /* copy buf to mybuf, adding exponent, sign and leading 0 */
 888     PyOS_snprintf(mybuf, mybuflen, "%s0%se%d", (sign ? "-" : ""),
 889                   buf, decpt - (int)buflen);
 890
 891     /* and convert the resulting string back to a double */
 892     errno = 0;
 893     rounded = _Py_dg_strtod(mybuf, NULL);
 894     if (errno == ERANGE && fabs(rounded) >= 1.)
 895         PyErr_SetString(PyExc_OverflowError,
 896                         "rounded value too large to represent");
 897     else
 898         result = PyFloat_FromDouble(rounded);
 899
 900     /* done computing value;  now clean up */
 901     if (mybuf != shortbuf)
 902         PyMem_Free(mybuf);
 903   exit:
 904     _Py_dg_freedtoa(buf);
 905     return result;
 906 }
 907
 908 #else /* PY_NO_SHORT_FLOAT_REPR */
 909
 910 /* fallback version, to be used when correctly rounded binary<->decimal
 911    conversions aren't available */
 912
 913 static PyObject *
 914 double_round(double x, int ndigits) {
 915     double pow1, pow2, y, z;
 916     if (ndigits >= 0) {
 917         if (ndigits > 22) {
 918             /* pow1 and pow2 are each safe from overflow, but
 919                pow1*pow2 ~= pow(10.0, ndigits) might overflow */
 920             pow1 = pow(10.0, (double)(ndigits-22));
 921             pow2 = 1e22;
 922         }
 923         else {
 924             pow1 = pow(10.0, (double)ndigits);
 925             pow2 = 1.0;
 926         }
 927         y = (x*pow1)*pow2;
 928         /* if y overflows, then rounded value is exactly x */
 929         if (!Py_IS_FINITE(y))
 930             return PyFloat_FromDouble(x);
 931     }
 932     else {
 933         pow1 = pow(10.0, (double)-ndigits);
 934         pow2 = 1.0; /* unused; silences a gcc compiler warning */
 935         y = x / pow1;
 936     }
 937
 938     z = round(y);
 939     if (fabs(y-z) == 0.5)
 940         /* halfway between two integers; use round-half-even */
 941         z = 2.0*round(y/2.0);
 942
 943     if (ndigits >= 0)
 944         z = (z / pow2) / pow1;
 945     else
 946         z *= pow1;
 947
 948     /* if computation resulted in overflow, raise OverflowError */
 949     if (!Py_IS_FINITE(z)) {
 950         PyErr_SetString(PyExc_OverflowError,
 951                         "overflow occurred during round");
 952         return NULL;
 953     }
 954
 955     return PyFloat_FromDouble(z);
 956 }
 957
 958 #endif /* PY_NO_SHORT_FLOAT_REPR */
 959
 960 /* round a Python float v to the closest multiple of 10**-ndigits */
 961
 962 static PyObject *
 963 float_round(PyObject *v, PyObject *args)
 964 {
 965     double x, rounded;
 966     PyObject *o_ndigits = NULL;
 967     Py_ssize_t ndigits;
 968
 969     x = PyFloat_AsDouble(v);
 970     if (!PyArg_ParseTuple(args, "|O", &o_ndigits))
 971         return NULL;
 972     if (o_ndigits == NULL) {
 973         /* single-argument round: round to nearest integer */
 974         rounded = round(x);
 975         if (fabs(x-rounded) == 0.5)
 976             /* halfway case: round to even */
 977             rounded = 2.0*round(x/2.0);
 978         return PyLong_FromDouble(rounded);
 979     }
 980
 981     /* interpret second argument as a Py_ssize_t; clips on overflow */
 982     ndigits = PyNumber_AsSsize_t(o_ndigits, NULL);
 983     if (ndigits == -1 && PyErr_Occurred())
 984         return NULL;
 985
 986     /* nans and infinities round to themselves */
 987     if (!Py_IS_FINITE(x))
 988         return PyFloat_FromDouble(x);
 989
 990     /* Deal with extreme values for ndigits. For ndigits > NDIGITS_MAX, x
 991        always rounds to itself.  For ndigits < NDIGITS_MIN, x always
 992        rounds to +-0.0.  Here 0.30103 is an upper bound for log10(2). */
 993 #define NDIGITS_MAX ((int)((DBL_MANT_DIG-DBL_MIN_EXP) * 0.30103))
 994 #define NDIGITS_MIN (-(int)((DBL_MAX_EXP + 1) * 0.30103))
 995     if (ndigits > NDIGITS_MAX)
 996         /* return x */
 997         return PyFloat_FromDouble(x);
 998     else if (ndigits < NDIGITS_MIN)
 999         /* return 0.0, but with sign of x */
1000         return PyFloat_FromDouble(0.0*x);
1001     else
1002         /* finite x, and ndigits is not unreasonably large */
1003         return double_round(x, (int)ndigits);
1004 #undef NDIGITS_MAX
1005 #undef NDIGITS_MIN
1006 }
1007
1008 static PyObject *
1009 float_float(PyObject *v)
1010 {
1011     if (PyFloat_CheckExact(v))
1012         Py_INCREF(v);
1013     else
1014         v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
1015     return v;
1016 }
1017
1018 /* turn ASCII hex characters into integer values and vice versa */
1019
1020 static char
1021 char_from_hex(int x)
1022 {
1023     assert(0 <= x && x < 16);
1024     return "0123456789abcdef"[x];
1025 }
1026
1027 static int
1028 hex_from_char(char c) {
1029     int x;
1030     switch(c) {
1031     case '0':
1032         x = 0;
1033         break;
1034     case '1':
1035         x = 1;
1036         break;
1037     case '2':
1038         x = 2;
1039         break;
1040     case '3':
1041         x = 3;
1042         break;
1043     case '4':
1044         x = 4;
1045         break;
1046     case '5':
1047         x = 5;
1048         break;
1049     case '6':
1050         x = 6;
1051         break;
1052     case '7':
1053         x = 7;
1054         break;
1055     case '8':
1056         x = 8;
1057         break;
1058     case '9':
1059         x = 9;
1060         break;
1061     case 'a':
1062     case 'A':
1063         x = 10;
1064         break;
1065     case 'b':
1066     case 'B':
1067         x = 11;
1068         break;
1069     case 'c':
1070     case 'C':
1071         x = 12;
1072         break;
1073     case 'd':
1074     case 'D':
1075         x = 13;
1076         break;
1077     case 'e':
1078     case 'E':
1079         x = 14;
1080         break;
1081     case 'f':
1082     case 'F':
1083         x = 15;
1084         break;
1085     default:
1086         x = -1;
1087         break;
1088     }
1089     return x;
1090 }
1091
1092 /* convert a float to a hexadecimal string */
1093
1094 /* TOHEX_NBITS is DBL_MANT_DIG rounded up to the next integer
1095    of the form 4k+1. */
1096 #define TOHEX_NBITS DBL_MANT_DIG + 3 - (DBL_MANT_DIG+2)%4
1097
1098 static PyObject *
1099 float_hex(PyObject *v)
1100 {
1101     double x, m;
1102     int e, shift, i, si, esign;
1103     /* Space for 1+(TOHEX_NBITS-1)/4 digits, a decimal point, and the
1104        trailing NUL byte. */
1105     char s[(TOHEX_NBITS-1)/4+3];
1106
1107     CONVERT_TO_DOUBLE(v, x);
1108
1109     if (Py_IS_NAN(x) || Py_IS_INFINITY(x))
1110         return float_str((PyFloatObject *)v);
1111
1112     if (x == 0.0) {
1113         if (copysign(1.0, x) == -1.0)
1114             return PyUnicode_FromString("-0x0.0p+0");
1115         else
1116             return PyUnicode_FromString("0x0.0p+0");
1117     }
1118
1119     m = frexp(fabs(x), &e);
1120     shift = 1 - MAX(DBL_MIN_EXP - e, 0);
1121     m = ldexp(m, shift);
1122     e -= shift;
1123
1124     si = 0;
1125     s[si] = char_from_hex((int)m);
1126     si++;
1127     m -= (int)m;
1128     s[si] = '.';
1129     si++;
1130     for (i=0; i < (TOHEX_NBITS-1)/4; i++) {
1131         m *= 16.0;
1132         s[si] = char_from_hex((int)m);
1133         si++;
1134         m -= (int)m;
1135     }
1136     s[si] = '\0';
1137
1138     if (e < 0) {
1139         esign = (int)'-';
1140         e = -e;
1141     }
1142     else
1143         esign = (int)'+';
1144
1145     if (x < 0.0)
1146         return PyUnicode_FromFormat("-0x%sp%c%d", s, esign, e);
1147     else
1148         return PyUnicode_FromFormat("0x%sp%c%d", s, esign, e);
1149 }
1150
1151 PyDoc_STRVAR(float_hex_doc,
1152 "float.hex() -> string\n\
1153 \n\
1154 Return a hexadecimal representation of a floating-point number.\n\
1155 >>> (-0.1).hex()\n\
1156 '-0x1.999999999999ap-4'\n\
1157 >>> 3.14159.hex()\n\
1158 '0x1.921f9f01b866ep+1'");
1159
1160 /* Convert a hexadecimal string to a float. */
1161
1162 static PyObject *
1163 float_fromhex(PyObject *cls, PyObject *arg)
1164 {
1165     PyObject *result_as_float, *result;
1166     double x;
1167     long exp, top_exp, lsb, key_digit;
1168     char *s, *coeff_start, *s_store, *coeff_end, *exp_start, *s_end;
1169     int half_eps, digit, round_up, negate=0;
1170     Py_ssize_t length, ndigits, fdigits, i;
1171
1172     /*
1173      * For the sake of simplicity and correctness, we impose an artificial
1174      * limit on ndigits, the total number of hex digits in the coefficient
1175      * The limit is chosen to ensure that, writing exp for the exponent,
1176      *
1177      *   (1) if exp > LONG_MAX/2 then the value of the hex string is
1178      *   guaranteed to overflow (provided it's nonzero)
1179      *
1180      *   (2) if exp < LONG_MIN/2 then the value of the hex string is
1181      *   guaranteed to underflow to 0.
1182      *
1183      *   (3) if LONG_MIN/2 <= exp <= LONG_MAX/2 then there's no danger of
1184      *   overflow in the calculation of exp and top_exp below.
1185      *
1186      * More specifically, ndigits is assumed to satisfy the following
1187      * inequalities:
1188      *
1189      *   4*ndigits <= DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2
1190      *   4*ndigits <= LONG_MAX/2 + 1 - DBL_MAX_EXP
1191      *
1192      * If either of these inequalities is not satisfied, a ValueError is
1193      * raised.  Otherwise, write x for the value of the hex string, and
1194      * assume x is nonzero.  Then
1195      *
1196      *   2**(exp-4*ndigits) <= |x| < 2**(exp+4*ndigits).
1197      *
1198      * Now if exp > LONG_MAX/2 then:
1199      *
1200      *   exp - 4*ndigits >= LONG_MAX/2 + 1 - (LONG_MAX/2 + 1 - DBL_MAX_EXP)
1201      *                    = DBL_MAX_EXP
1202      *
1203      * so |x| >= 2**DBL_MAX_EXP, which is too large to be stored in C
1204      * double, so overflows.  If exp < LONG_MIN/2, then
1205      *
1206      *   exp + 4*ndigits <= LONG_MIN/2 - 1 + (
1207      *                      DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2)
1208      *                    = DBL_MIN_EXP - DBL_MANT_DIG - 1
1209      *
1210      * and so |x| < 2**(DBL_MIN_EXP-DBL_MANT_DIG-1), hence underflows to 0
1211      * when converted to a C double.
1212      *
1213      * It's easy to show that if LONG_MIN/2 <= exp <= LONG_MAX/2 then both
1214      * exp+4*ndigits and exp-4*ndigits are within the range of a long.
1215      */
1216
1217     s = _PyUnicode_AsStringAndSize(arg, &length);
1218     if (s == NULL)
1219         return NULL;
1220     s_end = s + length;
1221
1222     /********************
1223      * Parse the string *
1224      ********************/
1225
1226     /* leading whitespace */
1227     while (Py_ISSPACE(*s))
1228         s++;
1229
1230     /* infinities and nans */
1231     x = _Py_parse_inf_or_nan(s, &coeff_end);
1232     if (coeff_end != s) {
1233         s = coeff_end;
1234         goto finished;
1235     }
1236
1237     /* optional sign */
1238     if (*s == '-') {
1239         s++;
1240         negate = 1;
1241     }
1242     else if (*s == '+')
1243         s++;
1244
1245     /* [0x] */
1246     s_store = s;
1247     if (*s == '0') {
1248         s++;
1249         if (*s == 'x' || *s == 'X')
1250             s++;
1251         else
1252             s = s_store;
1253     }
1254
1255     /* coefficient: <integer> [. <fraction>] */
1256     coeff_start = s;
1257     while (hex_from_char(*s) >= 0)
1258         s++;
1259     s_store = s;
1260     if (*s == '.') {
1261         s++;
1262         while (hex_from_char(*s) >= 0)
1263             s++;
1264         coeff_end = s-1;
1265     }
1266     else
1267         coeff_end = s;
1268
1269     /* ndigits = total # of hex digits; fdigits = # after point */
1270     ndigits = coeff_end - coeff_start;
1271     fdigits = coeff_end - s_store;
1272     if (ndigits == 0)
1273         goto parse_error;
1274     if (ndigits > MIN(DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2,
1275                       LONG_MAX/2 + 1 - DBL_MAX_EXP)/4)
1276         goto insane_length_error;
1277
1278     /* [p <exponent>] */
1279     if (*s == 'p' || *s == 'P') {
1280         s++;
1281         exp_start = s;
1282         if (*s == '-' || *s == '+')
1283             s++;
1284         if (!('0' <= *s && *s <= '9'))
1285             goto parse_error;
1286         s++;
1287         while ('0' <= *s && *s <= '9')
1288             s++;
1289         exp = strtol(exp_start, NULL, 10);
1290     }
1291     else
1292         exp = 0;
1293
1294 /* for 0 <= j < ndigits, HEX_DIGIT(j) gives the jth most significant digit */
1295 #define HEX_DIGIT(j) hex_from_char(*((j) < fdigits ?            \
1296                      coeff_end-(j) :                                    \
1297                      coeff_end-1-(j)))
1298
1299     /*******************************************
1300      * Compute rounded value of the hex string *
1301      *******************************************/
1302
1303     /* Discard leading zeros, and catch extreme overflow and underflow */
1304     while (ndigits > 0 && HEX_DIGIT(ndigits-1) == 0)
1305         ndigits--;
1306     if (ndigits == 0 || exp < LONG_MIN/2) {
1307         x = 0.0;
1308         goto finished;
1309     }
1310     if (exp > LONG_MAX/2)
1311         goto overflow_error;
1312
1313     /* Adjust exponent for fractional part. */
1314     exp = exp - 4*((long)fdigits);
1315
1316     /* top_exp = 1 more than exponent of most sig. bit of coefficient */
1317     top_exp = exp + 4*((long)ndigits - 1);
1318     for (digit = HEX_DIGIT(ndigits-1); digit != 0; digit /= 2)
1319         top_exp++;
1320
1321     /* catch almost all nonextreme cases of overflow and underflow here */
1322     if (top_exp < DBL_MIN_EXP - DBL_MANT_DIG) {
1323         x = 0.0;
1324         goto finished;
1325     }
1326     if (top_exp > DBL_MAX_EXP)
1327         goto overflow_error;
1328
1329     /* lsb = exponent of least significant bit of the *rounded* value.
1330        This is top_exp - DBL_MANT_DIG unless result is subnormal. */
1331     lsb = MAX(top_exp, (long)DBL_MIN_EXP) - DBL_MANT_DIG;
1332
1333     x = 0.0;
1334     if (exp >= lsb) {
1335         /* no rounding required */
1336         for (i = ndigits-1; i >= 0; i--)
1337             x = 16.0*x + HEX_DIGIT(i);
1338         x = ldexp(x, (int)(exp));
1339         goto finished;
1340     }
1341     /* rounding required.  key_digit is the index of the hex digit
1342        containing the first bit to be rounded away. */
1343     half_eps = 1 << (int)((lsb - exp - 1) % 4);
1344     key_digit = (lsb - exp - 1) / 4;
1345     for (i = ndigits-1; i > key_digit; i--)
1346         x = 16.0*x + HEX_DIGIT(i);
1347     digit = HEX_DIGIT(key_digit);
1348     x = 16.0*x + (double)(digit & (16-2*half_eps));
1349
1350     /* round-half-even: round up if bit lsb-1 is 1 and at least one of
1351        bits lsb, lsb-2, lsb-3, lsb-4, ... is 1. */
1352     if ((digit & half_eps) != 0) {
1353         round_up = 0;
1354         if ((digit & (3*half_eps-1)) != 0 ||
1355             (half_eps == 8 && (HEX_DIGIT(key_digit+1) & 1) != 0))
1356             round_up = 1;
1357         else
1358             for (i = key_digit-1; i >= 0; i--)
1359                 if (HEX_DIGIT(i) != 0) {
1360                     round_up = 1;
1361                     break;
1362                 }
1363         if (round_up == 1) {
1364             x += 2*half_eps;
1365             if (top_exp == DBL_MAX_EXP &&
1366                 x == ldexp((double)(2*half_eps), DBL_MANT_DIG))
1367                 /* overflow corner case: pre-rounded value <
1368                    2**DBL_MAX_EXP; rounded=2**DBL_MAX_EXP. */
1369                 goto overflow_error;
1370         }
1371     }
1372     x = ldexp(x, (int)(exp+4*key_digit));
1373
1374   finished:
1375     /* optional trailing whitespace leading to the end of the string */
1376     while (Py_ISSPACE(*s))
1377         s++;
1378     if (s != s_end)
1379         goto parse_error;
1380     result_as_float = Py_BuildValue("(d)", negate ? -x : x);
1381     if (result_as_float == NULL)
1382         return NULL;
1383     result = PyObject_CallObject(cls, result_as_float);
1384     Py_DECREF(result_as_float);
1385     return result;
1386
1387   overflow_error:
1388     PyErr_SetString(PyExc_OverflowError,
1389                     "hexadecimal value too large to represent as a float");
1390     return NULL;
1391
1392   parse_error:
1393     PyErr_SetString(PyExc_ValueError,
1394                     "invalid hexadecimal floating-point string");
1395     return NULL;
1396
1397   insane_length_error:
1398     PyErr_SetString(PyExc_ValueError,
1399                     "hexadecimal string too long to convert");
1400     return NULL;
1401 }
1402
1403 PyDoc_STRVAR(float_fromhex_doc,
1404 "float.fromhex(string) -> float\n\
1405 \n\
1406 Create a floating-point number from a hexadecimal string.\n\
1407 >>> float.fromhex('0x1.ffffp10')\n\
1408 2047.984375\n\
1409 >>> float.fromhex('-0x1p-1074')\n\
1410 -4.9406564584124654e-324");
1411
1412
1413 static PyObject *
1414 float_as_integer_ratio(PyObject *v, PyObject *unused)
1415 {
1416     double self;
1417     double float_part;
1418     int exponent;
1419     int i;
1420
1421     PyObject *prev;
1422     PyObject *py_exponent = NULL;
1423     PyObject *numerator = NULL;
1424     PyObject *denominator = NULL;
1425     PyObject *result_pair = NULL;
1426     PyNumberMethods *long_methods = PyLong_Type.tp_as_number;
1427
1428 #define INPLACE_UPDATE(obj, call) \
1429     prev = obj; \
1430     obj = call; \
1431     Py_DECREF(prev); \
1432
1433     CONVERT_TO_DOUBLE(v, self);
1434
1435     if (Py_IS_INFINITY(self)) {
1436       PyErr_SetString(PyExc_OverflowError,
1437                       "Cannot pass infinity to float.as_integer_ratio.");
1438       return NULL;
1439     }
1440 #ifdef Py_NAN
1441     if (Py_IS_NAN(self)) {
1442       PyErr_SetString(PyExc_ValueError,
1443                       "Cannot pass NaN to float.as_integer_ratio.");
1444       return NULL;
1445     }
1446 #endif
1447
1448     PyFPE_START_PROTECT("as_integer_ratio", goto error);
1449     float_part = frexp(self, &exponent);        /* self == float_part * 2**exponent exactly */
1450     PyFPE_END_PROTECT(float_part);
1451
1452     for (i=0; i<300 && float_part != floor(float_part) ; i++) {
1453         float_part *= 2.0;
1454         exponent--;
1455     }
1456     /* self == float_part * 2**exponent exactly and float_part is integral.
1457        If FLT_RADIX != 2, the 300 steps may leave a tiny fractional part
1458        to be truncated by PyLong_FromDouble(). */
1459
1460     numerator = PyLong_FromDouble(float_part);
1461     if (numerator == NULL) goto error;
1462
1463     /* fold in 2**exponent */
1464     denominator = PyLong_FromLong(1);
1465     py_exponent = PyLong_FromLong(labs((long)exponent));
1466     if (py_exponent == NULL) goto error;
1467     INPLACE_UPDATE(py_exponent,
1468                    long_methods->nb_lshift(denominator, py_exponent));
1469     if (py_exponent == NULL) goto error;
1470     if (exponent > 0) {
1471         INPLACE_UPDATE(numerator,
1472                        long_methods->nb_multiply(numerator, py_exponent));
1473         if (numerator == NULL) goto error;
1474     }
1475     else {
1476         Py_DECREF(denominator);
1477         denominator = py_exponent;
1478         py_exponent = NULL;
1479     }
1480
1481     result_pair = PyTuple_Pack(2, numerator, denominator);
1482
1483 #undef INPLACE_UPDATE
1484 error:
1485     Py_XDECREF(py_exponent);
1486     Py_XDECREF(denominator);
1487     Py_XDECREF(numerator);
1488     return result_pair;
1489 }
1490
1491 PyDoc_STRVAR(float_as_integer_ratio_doc,
1492 "float.as_integer_ratio() -> (int, int)\n"
1493 "\n"
1494 "Returns a pair of integers, whose ratio is exactly equal to the original\n"
1495 "float and with a positive denominator.\n"
1496 "Raises OverflowError on infinities and a ValueError on NaNs.\n"
1497 "\n"
1498 ">>> (10.0).as_integer_ratio()\n"
1499 "(10, 1)\n"
1500 ">>> (0.0).as_integer_ratio()\n"
1501 "(0, 1)\n"
1502 ">>> (-.25).as_integer_ratio()\n"
1503 "(-1, 4)");
1504
1505
1506 static PyObject *
1507 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
1508
1509 static PyObject *
1510 float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1511 {
1512     PyObject *x = Py_False; /* Integer zero */
1513     static char *kwlist[] = {"x", 0};
1514
1515     if (type != &PyFloat_Type)
1516         return float_subtype_new(type, args, kwds); /* Wimp out */
1517     if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
1518         return NULL;
1519     /* If it's a string, but not a string subclass, use
1520        PyFloat_FromString. */
1521     if (PyUnicode_CheckExact(x))
1522         return PyFloat_FromString(x);
1523     return PyNumber_Float(x);
1524 }
1525
1526 /* Wimpy, slow approach to tp_new calls for subtypes of float:
1527    first create a regular float from whatever arguments we got,
1528    then allocate a subtype instance and initialize its ob_fval
1529    from the regular float.  The regular float is then thrown away.
1530 */
1531 static PyObject *
1532 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1533 {
1534     PyObject *tmp, *newobj;
1535
1536     assert(PyType_IsSubtype(type, &PyFloat_Type));
1537     tmp = float_new(&PyFloat_Type, args, kwds);
1538     if (tmp == NULL)
1539         return NULL;
1540     assert(PyFloat_CheckExact(tmp));
1541     newobj = type->tp_alloc(type, 0);
1542     if (newobj == NULL) {
1543         Py_DECREF(tmp);
1544         return NULL;
1545     }
1546     ((PyFloatObject *)newobj)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
1547     Py_DECREF(tmp);
1548     return newobj;
1549 }
1550
1551 static PyObject *
1552 float_getnewargs(PyFloatObject *v)
1553 {
1554     return Py_BuildValue("(d)", v->ob_fval);
1555 }
1556
1557 /* this is for the benefit of the pack/unpack routines below */
1558
1559 typedef enum {
1560     unknown_format, ieee_big_endian_format, ieee_little_endian_format
1561 } float_format_type;
1562
1563 static float_format_type double_format, float_format;
1564 static float_format_type detected_double_format, detected_float_format;
1565
1566 static PyObject *
1567 float_getformat(PyTypeObject *v, PyObject* arg)
1568 {
1569     char* s;
1570     float_format_type r;
1571
1572     if (!PyUnicode_Check(arg)) {
1573         PyErr_Format(PyExc_TypeError,
1574          "__getformat__() argument must be string, not %.500s",
1575                          Py_TYPE(arg)->tp_name);
1576         return NULL;
1577     }
1578     s = _PyUnicode_AsString(arg);
1579     if (s == NULL)
1580         return NULL;
1581     if (strcmp(s, "double") == 0) {
1582         r = double_format;
1583     }
1584     else if (strcmp(s, "float") == 0) {
1585         r = float_format;
1586     }
1587     else {
1588         PyErr_SetString(PyExc_ValueError,
1589                         "__getformat__() argument 1 must be "
1590                         "'double' or 'float'");
1591         return NULL;
1592     }
1593
1594     switch (r) {
1595     case unknown_format:
1596         return PyUnicode_FromString("unknown");
1597     case ieee_little_endian_format:
1598         return PyUnicode_FromString("IEEE, little-endian");
1599     case ieee_big_endian_format:
1600         return PyUnicode_FromString("IEEE, big-endian");
1601     default:
1602         Py_FatalError("insane float_format or double_format");
1603         return NULL;
1604     }
1605 }
1606
1607 PyDoc_STRVAR(float_getformat_doc,
1608 "float.__getformat__(typestr) -> string\n"
1609 "\n"
1610 "You probably don't want to use this function.  It exists mainly to be\n"
1611 "used in Python's test suite.\n"
1612 "\n"
1613 "typestr must be 'double' or 'float'.  This function returns whichever of\n"
1614 "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1615 "format of floating point numbers used by the C type named by typestr.");
1616
1617 static PyObject *
1618 float_setformat(PyTypeObject *v, PyObject* args)
1619 {
1620     char* typestr;
1621     char* format;
1622     float_format_type f;
1623     float_format_type detected;
1624     float_format_type *p;
1625
1626     if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format))
1627         return NULL;
1628
1629     if (strcmp(typestr, "double") == 0) {
1630         p = &double_format;
1631         detected = detected_double_format;
1632     }
1633     else if (strcmp(typestr, "float") == 0) {
1634         p = &float_format;
1635         detected = detected_float_format;
1636     }
1637     else {
1638         PyErr_SetString(PyExc_ValueError,
1639                         "__setformat__() argument 1 must "
1640                         "be 'double' or 'float'");
1641         return NULL;
1642     }
1643
1644     if (strcmp(format, "unknown") == 0) {
1645         f = unknown_format;
1646     }
1647     else if (strcmp(format, "IEEE, little-endian") == 0) {
1648         f = ieee_little_endian_format;
1649     }
1650     else if (strcmp(format, "IEEE, big-endian") == 0) {
1651         f = ieee_big_endian_format;
1652     }
1653     else {
1654         PyErr_SetString(PyExc_ValueError,
1655                         "__setformat__() argument 2 must be "
1656                         "'unknown', 'IEEE, little-endian' or "
1657                         "'IEEE, big-endian'");
1658         return NULL;
1659
1660     }
1661
1662     if (f != unknown_format && f != detected) {
1663         PyErr_Format(PyExc_ValueError,
1664                      "can only set %s format to 'unknown' or the "
1665                      "detected platform value", typestr);
1666         return NULL;
1667     }
1668
1669     *p = f;
1670     Py_RETURN_NONE;
1671 }
1672
1673 PyDoc_STRVAR(float_setformat_doc,
1674 "float.__setformat__(typestr, fmt) -> None\n"
1675 "\n"
1676 "You probably don't want to use this function.  It exists mainly to be\n"
1677 "used in Python's test suite.\n"
1678 "\n"
1679 "typestr must be 'double' or 'float'.  fmt must be one of 'unknown',\n"
1680 "'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1681 "one of the latter two if it appears to match the underlying C reality.\n"
1682 "\n"
1683 "Overrides the automatic determination of C-level floating point type.\n"
1684 "This affects how floats are converted to and from binary strings.");
1685
1686 static PyObject *
1687 float_getzero(PyObject *v, void *closure)
1688 {
1689     return PyFloat_FromDouble(0.0);
1690 }
1691
1692 static PyObject *
1693 float__format__(PyObject *self, PyObject *args)
1694 {
1695     PyObject *format_spec;
1696
1697     if (!PyArg_ParseTuple(args, "U:__format__", &format_spec))
1698         return NULL;
1699     return _PyFloat_FormatAdvanced(self,
1700                                    PyUnicode_AS_UNICODE(format_spec),
1701                                    PyUnicode_GET_SIZE(format_spec));
1702 }
1703
1704 PyDoc_STRVAR(float__format__doc,
1705 "float.__format__(format_spec) -> string\n"
1706 "\n"
1707 "Formats the float according to format_spec.");
1708
1709
1710 static PyMethodDef float_methods[] = {
1711     {"conjugate",       (PyCFunction)float_float,       METH_NOARGS,
1712      "Returns self, the complex conjugate of any float."},
1713     {"__trunc__",       (PyCFunction)float_trunc, METH_NOARGS,
1714      "Returns the Integral closest to x between 0 and x."},
1715     {"__round__",       (PyCFunction)float_round, METH_VARARGS,
1716      "Returns the Integral closest to x, rounding half toward even.\n"
1717      "When an argument is passed, works like built-in round(x, ndigits)."},
1718     {"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
1719      float_as_integer_ratio_doc},
1720     {"fromhex", (PyCFunction)float_fromhex,
1721      METH_O|METH_CLASS, float_fromhex_doc},
1722     {"hex", (PyCFunction)float_hex,
1723      METH_NOARGS, float_hex_doc},
1724     {"is_integer",      (PyCFunction)float_is_integer,  METH_NOARGS,
1725      "Returns True if the float is an integer."},
1726 #if 0
1727     {"is_inf",          (PyCFunction)float_is_inf,      METH_NOARGS,
1728      "Returns True if the float is positive or negative infinite."},
1729     {"is_finite",       (PyCFunction)float_is_finite,   METH_NOARGS,
1730      "Returns True if the float is finite, neither infinite nor NaN."},
1731     {"is_nan",          (PyCFunction)float_is_nan,      METH_NOARGS,
1732      "Returns True if the float is not a number (NaN)."},
1733 #endif
1734     {"__getnewargs__",          (PyCFunction)float_getnewargs,  METH_NOARGS},
1735     {"__getformat__",           (PyCFunction)float_getformat,
1736      METH_O|METH_CLASS,                 float_getformat_doc},
1737     {"__setformat__",           (PyCFunction)float_setformat,
1738      METH_VARARGS|METH_CLASS,           float_setformat_doc},
1739     {"__format__",          (PyCFunction)float__format__,
1740      METH_VARARGS,                  float__format__doc},
1741     {NULL,              NULL}           /* sentinel */
1742 };
1743
1744 static PyGetSetDef float_getset[] = {
1745     {"real",
1746      (getter)float_float, (setter)NULL,
1747      "the real part of a complex number",
1748      NULL},
1749     {"imag",
1750      (getter)float_getzero, (setter)NULL,
1751      "the imaginary part of a complex number",
1752      NULL},
1753     {NULL}  /* Sentinel */
1754 };
1755
1756 PyDoc_STRVAR(float_doc,
1757 "float(x) -> floating point number\n\
1758 \n\
1759 Convert a string or number to a floating point number, if possible.");
1760
1761
1762 static PyNumberMethods float_as_number = {
1763     float_add,          /*nb_add*/
1764     float_sub,          /*nb_subtract*/
1765     float_mul,          /*nb_multiply*/
1766     float_rem,          /*nb_remainder*/
1767     float_divmod,       /*nb_divmod*/
1768     float_pow,          /*nb_power*/
1769     (unaryfunc)float_neg, /*nb_negative*/
1770     (unaryfunc)float_float, /*nb_positive*/
1771     (unaryfunc)float_abs, /*nb_absolute*/
1772     (inquiry)float_bool, /*nb_bool*/
1773     0,                  /*nb_invert*/
1774     0,                  /*nb_lshift*/
1775     0,                  /*nb_rshift*/
1776     0,                  /*nb_and*/
1777     0,                  /*nb_xor*/
1778     0,                  /*nb_or*/
1779     float_trunc,        /*nb_int*/
1780     0,                  /*nb_reserved*/
1781     float_float,        /*nb_float*/
1782     0,                  /* nb_inplace_add */
1783     0,                  /* nb_inplace_subtract */
1784     0,                  /* nb_inplace_multiply */
1785     0,                  /* nb_inplace_remainder */
1786     0,                  /* nb_inplace_power */
1787     0,                  /* nb_inplace_lshift */
1788     0,                  /* nb_inplace_rshift */
1789     0,                  /* nb_inplace_and */
1790     0,                  /* nb_inplace_xor */
1791     0,                  /* nb_inplace_or */
1792     float_floor_div, /* nb_floor_divide */
1793     float_div,          /* nb_true_divide */
1794     0,                  /* nb_inplace_floor_divide */
1795     0,                  /* nb_inplace_true_divide */
1796 };
1797
1798 PyTypeObject PyFloat_Type = {
1799     PyVarObject_HEAD_INIT(&PyType_Type, 0)
1800     "float",
1801     sizeof(PyFloatObject),
1802     0,
1803     (destructor)float_dealloc,                  /* tp_dealloc */
1804     0,                                          /* tp_print */
1805     0,                                          /* tp_getattr */
1806     0,                                          /* tp_setattr */
1807     0,                                          /* tp_reserved */
1808     (reprfunc)float_repr,                       /* tp_repr */
1809     &float_as_number,                           /* tp_as_number */
1810     0,                                          /* tp_as_sequence */
1811     0,                                          /* tp_as_mapping */
1812     (hashfunc)float_hash,                       /* tp_hash */
1813     0,                                          /* tp_call */
1814     (reprfunc)float_str,                        /* tp_str */
1815     PyObject_GenericGetAttr,                    /* tp_getattro */
1816     0,                                          /* tp_setattro */
1817     0,                                          /* tp_as_buffer */
1818     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags */
1819     float_doc,                                  /* tp_doc */
1820     0,                                          /* tp_traverse */
1821     0,                                          /* tp_clear */
1822     float_richcompare,                          /* tp_richcompare */
1823     0,                                          /* tp_weaklistoffset */
1824     0,                                          /* tp_iter */
1825     0,                                          /* tp_iternext */
1826     float_methods,                              /* tp_methods */
1827     0,                                          /* tp_members */
1828     float_getset,                               /* tp_getset */
1829     0,                                          /* tp_base */
1830     0,                                          /* tp_dict */
1831     0,                                          /* tp_descr_get */
1832     0,                                          /* tp_descr_set */
1833     0,                                          /* tp_dictoffset */
1834     0,                                          /* tp_init */
1835     0,                                          /* tp_alloc */
1836     float_new,                                  /* tp_new */
1837 };
1838
1839 void
1840 _PyFloat_Init(void)
1841 {
1842     /* We attempt to determine if this machine is using IEEE
1843        floating point formats by peering at the bits of some
1844        carefully chosen values.  If it looks like we are on an
1845        IEEE platform, the float packing/unpacking routines can
1846        just copy bits, if not they resort to arithmetic & shifts
1847        and masks.  The shifts & masks approach works on all finite
1848        values, but what happens to infinities, NaNs and signed
1849        zeroes on packing is an accident, and attempting to unpack
1850        a NaN or an infinity will raise an exception.
1851
1852        Note that if we're on some whacked-out platform which uses
1853        IEEE formats but isn't strictly little-endian or big-
1854        endian, we will fall back to the portable shifts & masks
1855        method. */
1856
1857 #if SIZEOF_DOUBLE == 8
1858     {
1859         double x = 9006104071832581.0;
1860         if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
1861             detected_double_format = ieee_big_endian_format;
1862         else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
1863             detected_double_format = ieee_little_endian_format;
1864         else
1865             detected_double_format = unknown_format;
1866     }
1867 #else
1868     detected_double_format = unknown_format;
1869 #endif
1870
1871 #if SIZEOF_FLOAT == 4
1872     {
1873         float y = 16711938.0;
1874         if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
1875             detected_float_format = ieee_big_endian_format;
1876         else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
1877             detected_float_format = ieee_little_endian_format;
1878         else
1879             detected_float_format = unknown_format;
1880     }
1881 #else
1882     detected_float_format = unknown_format;
1883 #endif
1884
1885     double_format = detected_double_format;
1886     float_format = detected_float_format;
1887
1888     /* Init float info */
1889     if (FloatInfoType.tp_name == 0)
1890         PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);
1891 }
1892
1893 int
1894 PyFloat_ClearFreeList(void)
1895 {
1896     PyFloatObject *p;
1897     PyFloatBlock *list, *next;
1898     int i;
1899     int u;                      /* remaining unfreed floats per block */
1900     int freelist_size = 0;
1901
1902     list = block_list;
1903     block_list = NULL;
1904     free_list = NULL;
1905     while (list != NULL) {
1906         u = 0;
1907         for (i = 0, p = &list->objects[0];
1908              i < N_FLOATOBJECTS;
1909              i++, p++) {
1910             if (PyFloat_CheckExact(p) && Py_REFCNT(p) != 0)
1911                 u++;
1912         }
1913         next = list->next;
1914         if (u) {
1915             list->next = block_list;
1916             block_list = list;
1917             for (i = 0, p = &list->objects[0];
1918                  i < N_FLOATOBJECTS;
1919                  i++, p++) {
1920                 if (!PyFloat_CheckExact(p) ||
1921                     Py_REFCNT(p) == 0) {
1922                     Py_TYPE(p) = (struct _typeobject *)
1923                         free_list;
1924                     free_list = p;
1925                 }
1926             }
1927         }
1928         else {
1929             PyMem_FREE(list);
1930         }
1931         freelist_size += u;
1932         list = next;
1933     }
1934     return freelist_size;
1935 }
1936
1937 void
1938 PyFloat_Fini(void)
1939 {
1940     PyFloatObject *p;
1941     PyFloatBlock *list;
1942     int i;
1943     int u;                      /* total unfreed floats per block */
1944
1945     u = PyFloat_ClearFreeList();
1946
1947     if (!Py_VerboseFlag)
1948         return;
1949     fprintf(stderr, "# cleanup floats");
1950     if (!u) {
1951         fprintf(stderr, "\n");
1952     }
1953     else {
1954         fprintf(stderr,
1955             ": %d unfreed float%s\n",
1956             u, u == 1 ? "" : "s");
1957     }
1958     if (Py_VerboseFlag > 1) {
1959         list = block_list;
1960         while (list != NULL) {
1961             for (i = 0, p = &list->objects[0];
1962                  i < N_FLOATOBJECTS;
1963                  i++, p++) {
1964                 if (PyFloat_CheckExact(p) &&
1965                     Py_REFCNT(p) != 0) {
1966                     char *buf = PyOS_double_to_string(
1967                         PyFloat_AS_DOUBLE(p), 'r',
1968                         0, 0, NULL);
1969                     if (buf) {
1970                         /* XXX(twouters) cast
1971                            refcount to long
1972                            until %zd is
1973                            universally
1974                            available
1975                         */
1976                         fprintf(stderr,
1977                  "#   <float at %p, refcnt=%ld, val=%s>\n",
1978                                     p, (long)Py_REFCNT(p), buf);
1979                                     PyMem_Free(buf);
1980                             }
1981                 }
1982             }
1983             list = list->next;
1984         }
1985     }
1986 }
1987
1988 /*----------------------------------------------------------------------------
1989  * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
1990  */
1991 int
1992 _PyFloat_Pack4(double x, unsigned char *p, int le)
1993 {
1994     if (float_format == unknown_format) {
1995         unsigned char sign;
1996         int e;
1997         double f;
1998         unsigned int fbits;
1999         int incr = 1;
2000
2001         if (le) {
2002             p += 3;
2003             incr = -1;
2004         }
2005
2006         if (x < 0) {
2007             sign = 1;
2008             x = -x;
2009         }
2010         else
2011             sign = 0;
2012
2013         f = frexp(x, &e);
2014
2015         /* Normalize f to be in the range [1.0, 2.0) */
2016         if (0.5 <= f && f < 1.0) {
2017             f *= 2.0;
2018             e--;
2019         }
2020         else if (f == 0.0)
2021             e = 0;
2022         else {
2023             PyErr_SetString(PyExc_SystemError,
2024                             "frexp() result out of range");
2025             return -1;
2026         }
2027
2028         if (e >= 128)
2029             goto Overflow;
2030         else if (e < -126) {
2031             /* Gradual underflow */
2032             f = ldexp(f, 126 + e);
2033             e = 0;
2034         }
2035         else if (!(e == 0 && f == 0.0)) {
2036             e += 127;
2037             f -= 1.0; /* Get rid of leading 1 */
2038         }
2039
2040         f *= 8388608.0; /* 2**23 */
2041         fbits = (unsigned int)(f + 0.5); /* Round */
2042         assert(fbits <= 8388608);
2043         if (fbits >> 23) {
2044             /* The carry propagated out of a string of 23 1 bits. */
2045             fbits = 0;
2046             ++e;
2047             if (e >= 255)
2048                 goto Overflow;
2049         }
2050
2051         /* First byte */
2052         *p = (sign << 7) | (e >> 1);
2053         p += incr;
2054
2055         /* Second byte */
2056         *p = (char) (((e & 1) << 7) | (fbits >> 16));
2057         p += incr;
2058
2059         /* Third byte */
2060         *p = (fbits >> 8) & 0xFF;
2061         p += incr;
2062
2063         /* Fourth byte */
2064         *p = fbits & 0xFF;
2065
2066         /* Done */
2067         return 0;
2068
2069     }
2070     else {
2071         float y = (float)x;
2072         const char *s = (char*)&y;
2073         int i, incr = 1;
2074
2075         if (Py_IS_INFINITY(y) && !Py_IS_INFINITY(x))
2076             goto Overflow;
2077
2078         if ((float_format == ieee_little_endian_format && !le)
2079             || (float_format == ieee_big_endian_format && le)) {
2080             p += 3;
2081             incr = -1;
2082         }
2083
2084         for (i = 0; i < 4; i++) {
2085             *p = *s++;
2086             p += incr;
2087         }
2088         return 0;
2089     }
2090   Overflow:
2091     PyErr_SetString(PyExc_OverflowError,
2092                     "float too large to pack with f format");
2093     return -1;
2094 }
2095
2096 int
2097 _PyFloat_Pack8(double x, unsigned char *p, int le)
2098 {
2099     if (double_format == unknown_format) {
2100         unsigned char sign;
2101         int e;
2102         double f;
2103         unsigned int fhi, flo;
2104         int incr = 1;
2105
2106         if (le) {
2107             p += 7;
2108             incr = -1;
2109         }
2110
2111         if (x < 0) {
2112             sign = 1;
2113             x = -x;
2114         }
2115         else
2116             sign = 0;
2117
2118         f = frexp(x, &e);
2119
2120         /* Normalize f to be in the range [1.0, 2.0) */
2121         if (0.5 <= f && f < 1.0) {
2122             f *= 2.0;
2123             e--;
2124         }
2125         else if (f == 0.0)
2126             e = 0;
2127         else {
2128             PyErr_SetString(PyExc_SystemError,
2129                             "frexp() result out of range");
2130             return -1;
2131         }
2132
2133         if (e >= 1024)
2134             goto Overflow;
2135         else if (e < -1022) {
2136             /* Gradual underflow */
2137             f = ldexp(f, 1022 + e);
2138             e = 0;
2139         }
2140         else if (!(e == 0 && f == 0.0)) {
2141             e += 1023;
2142             f -= 1.0; /* Get rid of leading 1 */
2143         }
2144
2145         /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
2146         f *= 268435456.0; /* 2**28 */
2147         fhi = (unsigned int)f; /* Truncate */
2148         assert(fhi < 268435456);
2149
2150         f -= (double)fhi;
2151         f *= 16777216.0; /* 2**24 */
2152         flo = (unsigned int)(f + 0.5); /* Round */
2153         assert(flo <= 16777216);
2154         if (flo >> 24) {
2155             /* The carry propagated out of a string of 24 1 bits. */
2156             flo = 0;
2157             ++fhi;
2158             if (fhi >> 28) {
2159                 /* And it also progagated out of the next 28 bits. */
2160                 fhi = 0;
2161                 ++e;
2162                 if (e >= 2047)
2163                     goto Overflow;
2164             }
2165         }
2166
2167         /* First byte */
2168         *p = (sign << 7) | (e >> 4);
2169         p += incr;
2170
2171         /* Second byte */
2172         *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
2173         p += incr;
2174
2175         /* Third byte */
2176         *p = (fhi >> 16) & 0xFF;
2177         p += incr;
2178
2179         /* Fourth byte */
2180         *p = (fhi >> 8) & 0xFF;
2181         p += incr;
2182
2183         /* Fifth byte */
2184         *p = fhi & 0xFF;
2185         p += incr;
2186
2187         /* Sixth byte */
2188         *p = (flo >> 16) & 0xFF;
2189         p += incr;
2190
2191         /* Seventh byte */
2192         *p = (flo >> 8) & 0xFF;
2193         p += incr;
2194
2195         /* Eighth byte */
2196         *p = flo & 0xFF;
2197         p += incr;
2198
2199         /* Done */
2200         return 0;
2201
2202       Overflow:
2203         PyErr_SetString(PyExc_OverflowError,
2204                         "float too large to pack with d format");
2205         return -1;
2206     }
2207     else {
2208         const char *s = (char*)&x;
2209         int i, incr = 1;
2210
2211         if ((double_format == ieee_little_endian_format && !le)
2212             || (double_format == ieee_big_endian_format && le)) {
2213             p += 7;
2214             incr = -1;
2215         }
2216
2217         for (i = 0; i < 8; i++) {
2218             *p = *s++;
2219             p += incr;
2220         }
2221         return 0;
2222     }
2223 }
2224
2225 double
2226 _PyFloat_Unpack4(const unsigned char *p, int le)
2227 {
2228     if (float_format == unknown_format) {
2229         unsigned char sign;
2230         int e;
2231         unsigned int f;
2232         double x;
2233         int incr = 1;
2234
2235         if (le) {
2236             p += 3;
2237             incr = -1;
2238         }
2239
2240         /* First byte */
2241         sign = (*p >> 7) & 1;
2242         e = (*p & 0x7F) << 1;
2243         p += incr;
2244
2245         /* Second byte */
2246         e |= (*p >> 7) & 1;
2247         f = (*p & 0x7F) << 16;
2248         p += incr;
2249
2250         if (e == 255) {
2251             PyErr_SetString(
2252                 PyExc_ValueError,
2253                 "can't unpack IEEE 754 special value "
2254                 "on non-IEEE platform");
2255             return -1;
2256         }
2257
2258         /* Third byte */
2259         f |= *p << 8;
2260         p += incr;
2261
2262         /* Fourth byte */
2263         f |= *p;
2264
2265         x = (double)f / 8388608.0;
2266
2267         /* XXX This sadly ignores Inf/NaN issues */
2268         if (e == 0)
2269             e = -126;
2270         else {
2271             x += 1.0;
2272             e -= 127;
2273         }
2274         x = ldexp(x, e);
2275
2276         if (sign)
2277             x = -x;
2278
2279         return x;
2280     }
2281     else {
2282         float x;
2283
2284         if ((float_format == ieee_little_endian_format && !le)
2285             || (float_format == ieee_big_endian_format && le)) {
2286             char buf[4];
2287             char *d = &buf[3];
2288             int i;
2289
2290             for (i = 0; i < 4; i++) {
2291                 *d-- = *p++;
2292             }
2293             memcpy(&x, buf, 4);
2294         }
2295         else {
2296             memcpy(&x, p, 4);
2297         }
2298
2299         return x;
2300     }
2301 }
2302
2303 double
2304 _PyFloat_Unpack8(const unsigned char *p, int le)
2305 {
2306     if (double_format == unknown_format) {
2307         unsigned char sign;
2308         int e;
2309         unsigned int fhi, flo;
2310         double x;
2311         int incr = 1;
2312
2313         if (le) {
2314             p += 7;
2315             incr = -1;
2316         }
2317
2318         /* First byte */
2319         sign = (*p >> 7) & 1;
2320         e = (*p & 0x7F) << 4;
2321
2322         p += incr;
2323
2324         /* Second byte */
2325         e |= (*p >> 4) & 0xF;
2326         fhi = (*p & 0xF) << 24;
2327         p += incr;
2328
2329         if (e == 2047) {
2330             PyErr_SetString(
2331                 PyExc_ValueError,
2332                 "can't unpack IEEE 754 special value "
2333                 "on non-IEEE platform");
2334             return -1.0;
2335         }
2336
2337         /* Third byte */
2338         fhi |= *p << 16;
2339         p += incr;
2340
2341         /* Fourth byte */
2342         fhi |= *p  << 8;
2343         p += incr;
2344
2345         /* Fifth byte */
2346         fhi |= *p;
2347         p += incr;
2348
2349         /* Sixth byte */
2350         flo = *p << 16;
2351         p += incr;
2352
2353         /* Seventh byte */
2354         flo |= *p << 8;
2355         p += incr;
2356
2357         /* Eighth byte */
2358         flo |= *p;
2359
2360         x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
2361         x /= 268435456.0; /* 2**28 */
2362
2363         if (e == 0)
2364             e = -1022;
2365         else {
2366             x += 1.0;
2367             e -= 1023;
2368         }
2369         x = ldexp(x, e);
2370
2371         if (sign)
2372             x = -x;
2373
2374         return x;
2375     }
2376     else {
2377         double x;
2378
2379         if ((double_format == ieee_little_endian_format && !le)
2380             || (double_format == ieee_big_endian_format && le)) {
2381             char buf[8];
2382             char *d = &buf[7];
2383             int i;
2384
2385             for (i = 0; i < 8; i++) {
2386                 *d-- = *p++;
2387             }
2388             memcpy(&x, buf, 8);
2389         }
2390         else {
2391             memcpy(&x, p, 8);
2392         }
2393
2394         return x;
2395     }
2396 }