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