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_long(PyObject *v)
1109 {
1110         double x = PyFloat_AsDouble(v);
1111         return PyLong_FromDouble(x);
1112 }
1113
1114 static PyObject *
1115 float_float(PyObject *v)
1116 {
1117         if (PyFloat_CheckExact(v))
1118                 Py_INCREF(v);
1119         else
1120                 v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
1121         return v;
1122 }
1123
1124 /* turn ASCII hex characters into integer values and vice versa */
1125
1126 static char
1127 char_from_hex(int x)
1128 {
1129         assert(0 <= x && x < 16);
1130         return "0123456789abcdef"[x];
1131 }
1132
1133 static int
1134 hex_from_char(char c) {
1135         int x;
1136         switch(c) {
1137         case '0':
1138                 x = 0;
1139                 break;
1140         case '1':
1141                 x = 1;
1142                 break;
1143         case '2':
1144                 x = 2;
1145                 break;
1146         case '3':
1147                 x = 3;
1148                 break;
1149         case '4':
1150                 x = 4;
1151                 break;
1152         case '5':
1153                 x = 5;
1154                 break;
1155         case '6':
1156                 x = 6;
1157                 break;
1158         case '7':
1159                 x = 7;
1160                 break;
1161         case '8':
1162                 x = 8;
1163                 break;
1164         case '9':
1165                 x = 9;
1166                 break;
1167         case 'a':
1168         case 'A':
1169                 x = 10;
1170                 break;
1171         case 'b':
1172         case 'B':
1173                 x = 11;
1174                 break;
1175         case 'c':
1176         case 'C':
1177                 x = 12;
1178                 break;
1179         case 'd':
1180         case 'D':
1181                 x = 13;
1182                 break;
1183         case 'e':
1184         case 'E':
1185                 x = 14;
1186                 break;
1187         case 'f':
1188         case 'F':
1189                 x = 15;
1190                 break;
1191         default:
1192                 x = -1;
1193                 break;
1194         }
1195         return x;
1196 }
1197
1198 /* convert a float to a hexadecimal string */
1199
1200 /* TOHEX_NBITS is DBL_MANT_DIG rounded up to the next integer
1201    of the form 4k+1. */
1202 #define TOHEX_NBITS DBL_MANT_DIG + 3 - (DBL_MANT_DIG+2)%4
1203
1204 static PyObject *
1205 float_hex(PyObject *v)
1206 {
1207         double x, m;
1208         int e, shift, i, si, esign;
1209         /* Space for 1+(TOHEX_NBITS-1)/4 digits, a decimal point, and the
1210            trailing NUL byte. */
1211         char s[(TOHEX_NBITS-1)/4+3];
1212
1213         CONVERT_TO_DOUBLE(v, x);
1214
1215         if (Py_IS_NAN(x) || Py_IS_INFINITY(x))
1216                 return float_str((PyFloatObject *)v);
1217
1218         if (x == 0.0) {
1219                 if(copysign(1.0, x) == -1.0)
1220                         return PyString_FromString("-0x0.0p+0");
1221                 else
1222                         return PyString_FromString("0x0.0p+0");
1223         }
1224
1225         m = frexp(fabs(x), &e);
1226         shift = 1 - MAX(DBL_MIN_EXP - e, 0);
1227         m = ldexp(m, shift);
1228         e -= shift;
1229
1230         si = 0;
1231         s[si] = char_from_hex((int)m);
1232         si++;
1233         m -= (int)m;
1234         s[si] = '.';
1235         si++;
1236         for (i=0; i < (TOHEX_NBITS-1)/4; i++) {
1237                 m *= 16.0;
1238                 s[si] = char_from_hex((int)m);
1239                 si++;
1240                 m -= (int)m;
1241         }
1242         s[si] = '\0';
1243
1244         if (e < 0) {
1245                 esign = (int)'-';
1246                 e = -e;
1247         }
1248         else
1249                 esign = (int)'+';
1250
1251         if (x < 0.0)
1252                 return PyString_FromFormat("-0x%sp%c%d", s, esign, e);
1253         else
1254                 return PyString_FromFormat("0x%sp%c%d", s, esign, e);
1255 }
1256
1257 PyDoc_STRVAR(float_hex_doc,
1258 "float.hex() -> string\n\
1259 \n\
1260 Return a hexadecimal representation of a floating-point number.\n\
1261 >>> (-0.1).hex()\n\
1262 '-0x1.999999999999ap-4'\n\
1263 >>> 3.14159.hex()\n\
1264 '0x1.921f9f01b866ep+1'");
1265
1266 /* Convert a hexadecimal string to a float. */
1267
1268 static PyObject *
1269 float_fromhex(PyObject *cls, PyObject *arg)
1270 {
1271         PyObject *result_as_float, *result;
1272         double x;
1273         long exp, top_exp, lsb, key_digit;
1274         char *s, *coeff_start, *s_store, *coeff_end, *exp_start, *s_end;
1275         int half_eps, digit, round_up, sign=1;
1276         Py_ssize_t length, ndigits, fdigits, i;
1277
1278         /*
1279          * For the sake of simplicity and correctness, we impose an artificial
1280          * limit on ndigits, the total number of hex digits in the coefficient
1281          * The limit is chosen to ensure that, writing exp for the exponent,
1282          *
1283          *   (1) if exp > LONG_MAX/2 then the value of the hex string is
1284          *   guaranteed to overflow (provided it's nonzero)
1285          *
1286          *   (2) if exp < LONG_MIN/2 then the value of the hex string is
1287          *   guaranteed to underflow to 0.
1288          *
1289          *   (3) if LONG_MIN/2 <= exp <= LONG_MAX/2 then there's no danger of
1290          *   overflow in the calculation of exp and top_exp below.
1291          *
1292          * More specifically, ndigits is assumed to satisfy the following
1293          * inequalities:
1294          *
1295          *   4*ndigits <= DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2
1296          *   4*ndigits <= LONG_MAX/2 + 1 - DBL_MAX_EXP
1297          *
1298          * If either of these inequalities is not satisfied, a ValueError is
1299          * raised.  Otherwise, write x for the value of the hex string, and
1300          * assume x is nonzero.  Then
1301          *
1302          *   2**(exp-4*ndigits) <= |x| < 2**(exp+4*ndigits).
1303          *
1304          * Now if exp > LONG_MAX/2 then:
1305          *
1306          *   exp - 4*ndigits >= LONG_MAX/2 + 1 - (LONG_MAX/2 + 1 - DBL_MAX_EXP)
1307          *                    = DBL_MAX_EXP
1308          *
1309          * so |x| >= 2**DBL_MAX_EXP, which is too large to be stored in C
1310          * double, so overflows.  If exp < LONG_MIN/2, then
1311          *
1312          *   exp + 4*ndigits <= LONG_MIN/2 - 1 + (
1313          *                      DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2)
1314          *                    = DBL_MIN_EXP - DBL_MANT_DIG - 1
1315          *
1316          * and so |x| < 2**(DBL_MIN_EXP-DBL_MANT_DIG-1), hence underflows to 0
1317          * when converted to a C double.
1318          *
1319          * It's easy to show that if LONG_MIN/2 <= exp <= LONG_MAX/2 then both
1320          * exp+4*ndigits and exp-4*ndigits are within the range of a long.
1321          */
1322
1323         if (PyString_AsStringAndSize(arg, &s, &length))
1324                 return NULL;
1325         s_end = s + length;
1326
1327         /********************
1328          * Parse the string *
1329          ********************/
1330
1331         /* leading whitespace and optional sign */
1332         while (isspace(*s))
1333                 s++;
1334         if (*s == '-') {
1335                 s++;
1336                 sign = -1;
1337         }
1338         else if (*s == '+')
1339                 s++;
1340
1341         /* infinities and nans */
1342         if (PyOS_strnicmp(s, "nan", 4) == 0) {
1343                 x = Py_NAN;
1344                 goto finished;
1345         }
1346         if (PyOS_strnicmp(s, "inf", 4) == 0 ||
1347             PyOS_strnicmp(s, "infinity", 9) == 0) {
1348                 x = sign*Py_HUGE_VAL;
1349                 goto finished;
1350         }
1351
1352         /* [0x] */
1353         s_store = s;
1354         if (*s == '0') {
1355                 s++;
1356                 if (tolower(*s) == (int)'x')
1357                         s++;
1358                 else
1359                         s = s_store;
1360         }
1361
1362         /* coefficient: <integer> [. <fraction>] */
1363         coeff_start = s;
1364         while (hex_from_char(*s) >= 0)
1365                 s++;
1366         s_store = s;
1367         if (*s == '.') {
1368                 s++;
1369                 while (hex_from_char(*s) >= 0)
1370                         s++;
1371                 coeff_end = s-1;
1372         }
1373         else
1374                 coeff_end = s;
1375
1376         /* ndigits = total # of hex digits; fdigits = # after point */
1377         ndigits = coeff_end - coeff_start;
1378         fdigits = coeff_end - s_store;
1379         if (ndigits == 0)
1380                 goto parse_error;
1381         if (ndigits > MIN(DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2,
1382                           LONG_MAX/2 + 1 - DBL_MAX_EXP)/4)
1383                 goto insane_length_error;
1384
1385         /* [p <exponent>] */
1386         if (tolower(*s) == (int)'p') {
1387                 s++;
1388                 exp_start = s;
1389                 if (*s == '-' || *s == '+')
1390                         s++;
1391                 if (!('0' <= *s && *s <= '9'))
1392                         goto parse_error;
1393                 s++;
1394                 while ('0' <= *s && *s <= '9')
1395                         s++;
1396                 exp = strtol(exp_start, NULL, 10);
1397         }
1398         else
1399                 exp = 0;
1400
1401         /* optional trailing whitespace leading to the end of the string */
1402         while (isspace(*s))
1403                 s++;
1404         if (s != s_end)
1405                 goto parse_error;
1406
1407 /* for 0 <= j < ndigits, HEX_DIGIT(j) gives the jth most significant digit */
1408 #define HEX_DIGIT(j) hex_from_char(*((j) < fdigits ?            \
1409                                      coeff_end-(j) :                    \
1410                                      coeff_end-1-(j)))
1411
1412         /*******************************************
1413          * Compute rounded value of the hex string *
1414          *******************************************/
1415
1416         /* Discard leading zeros, and catch extreme overflow and underflow */
1417         while (ndigits > 0 && HEX_DIGIT(ndigits-1) == 0)
1418                 ndigits--;
1419         if (ndigits == 0 || exp < LONG_MIN/2) {
1420                 x = sign * 0.0;
1421                 goto finished;
1422         }
1423         if (exp > LONG_MAX/2)
1424                 goto overflow_error;
1425
1426         /* Adjust exponent for fractional part. */
1427         exp = exp - 4*((long)fdigits);
1428
1429         /* top_exp = 1 more than exponent of most sig. bit of coefficient */
1430         top_exp = exp + 4*((long)ndigits - 1);
1431         for (digit = HEX_DIGIT(ndigits-1); digit != 0; digit /= 2)
1432                 top_exp++;
1433
1434         /* catch almost all nonextreme cases of overflow and underflow here */
1435         if (top_exp < DBL_MIN_EXP - DBL_MANT_DIG) {
1436                 x = sign * 0.0;
1437                 goto finished;
1438         }
1439         if (top_exp > DBL_MAX_EXP)
1440                 goto overflow_error;
1441
1442         /* lsb = exponent of least significant bit of the *rounded* value.
1443            This is top_exp - DBL_MANT_DIG unless result is subnormal. */
1444         lsb = MAX(top_exp, (long)DBL_MIN_EXP) - DBL_MANT_DIG;
1445
1446         x = 0.0;
1447         if (exp >= lsb) {
1448                 /* no rounding required */
1449                 for (i = ndigits-1; i >= 0; i--)
1450                         x = 16.0*x + HEX_DIGIT(i);
1451                 x = sign * ldexp(x, (int)(exp));
1452                 goto finished;
1453         }
1454         /* rounding required.  key_digit is the index of the hex digit
1455            containing the first bit to be rounded away. */
1456         half_eps = 1 << (int)((lsb - exp - 1) % 4);
1457         key_digit = (lsb - exp - 1) / 4;
1458         for (i = ndigits-1; i > key_digit; i--)
1459                 x = 16.0*x + HEX_DIGIT(i);
1460         digit = HEX_DIGIT(key_digit);
1461         x = 16.0*x + (double)(digit & (16-2*half_eps));
1462
1463         /* round-half-even: round up if bit lsb-1 is 1 and at least one of
1464            bits lsb, lsb-2, lsb-3, lsb-4, ... is 1. */
1465         if ((digit & half_eps) != 0) {
1466                 round_up = 0;
1467                 if ((digit & (3*half_eps-1)) != 0 ||
1468                     (half_eps == 8 && (HEX_DIGIT(key_digit+1) & 1) != 0))
1469                         round_up = 1;
1470                 else
1471                         for (i = key_digit-1; i >= 0; i--)
1472                                 if (HEX_DIGIT(i) != 0) {
1473                                         round_up = 1;
1474                                         break;
1475                                 }
1476                 if (round_up == 1) {
1477                         x += 2*half_eps;
1478                         if (top_exp == DBL_MAX_EXP &&
1479                             x == ldexp((double)(2*half_eps), DBL_MANT_DIG))
1480                                 /* overflow corner case: pre-rounded value <
1481                                    2**DBL_MAX_EXP; rounded=2**DBL_MAX_EXP. */
1482                                 goto overflow_error;
1483                 }
1484         }
1485         x = sign * ldexp(x, (int)(exp+4*key_digit));
1486
1487   finished:
1488         result_as_float = Py_BuildValue("(d)", x);
1489         if (result_as_float == NULL)
1490                 return NULL;
1491         result = PyObject_CallObject(cls, result_as_float);
1492         Py_DECREF(result_as_float);
1493         return result;
1494
1495   overflow_error:
1496         PyErr_SetString(PyExc_OverflowError,
1497                         "hexadecimal value too large to represent as a float");
1498         return NULL;
1499
1500   parse_error:
1501         PyErr_SetString(PyExc_ValueError,
1502                         "invalid hexadecimal floating-point string");
1503         return NULL;
1504
1505   insane_length_error:
1506         PyErr_SetString(PyExc_ValueError,
1507                         "hexadecimal string too long to convert");
1508         return NULL;
1509 }
1510
1511 PyDoc_STRVAR(float_fromhex_doc,
1512 "float.fromhex(string) -> float\n\
1513 \n\
1514 Create a floating-point number from a hexadecimal string.\n\
1515 >>> float.fromhex('0x1.ffffp10')\n\
1516 2047.984375\n\
1517 >>> float.fromhex('-0x1p-1074')\n\
1518 -4.9406564584124654e-324");
1519
1520
1521 static PyObject *
1522 float_as_integer_ratio(PyObject *v, PyObject *unused)
1523 {
1524         double self;
1525         double float_part;
1526         int exponent;
1527         int i;
1528
1529         PyObject *prev;
1530         PyObject *py_exponent = NULL;
1531         PyObject *numerator = NULL;
1532         PyObject *denominator = NULL;
1533         PyObject *result_pair = NULL;
1534         PyNumberMethods *long_methods = PyLong_Type.tp_as_number;
1535
1536 #define INPLACE_UPDATE(obj, call) \
1537         prev = obj; \
1538         obj = call; \
1539         Py_DECREF(prev); \
1540
1541         CONVERT_TO_DOUBLE(v, self);
1542
1543         if (Py_IS_INFINITY(self)) {
1544           PyErr_SetString(PyExc_OverflowError,
1545                           "Cannot pass infinity to float.as_integer_ratio.");
1546           return NULL;
1547         }
1548 #ifdef Py_NAN
1549         if (Py_IS_NAN(self)) {
1550           PyErr_SetString(PyExc_ValueError,
1551                           "Cannot pass NaN to float.as_integer_ratio.");
1552           return NULL;
1553         }
1554 #endif
1555
1556         PyFPE_START_PROTECT("as_integer_ratio", goto error);
1557         float_part = frexp(self, &exponent);    /* self == float_part * 2**exponent exactly */
1558         PyFPE_END_PROTECT(float_part);
1559
1560         for (i=0; i<300 && float_part != floor(float_part) ; i++) {
1561                 float_part *= 2.0;
1562                 exponent--;
1563         }
1564         /* self == float_part * 2**exponent exactly and float_part is integral.
1565            If FLT_RADIX != 2, the 300 steps may leave a tiny fractional part
1566            to be truncated by PyLong_FromDouble(). */
1567
1568         numerator = PyLong_FromDouble(float_part);
1569         if (numerator == NULL) goto error;
1570
1571         /* fold in 2**exponent */
1572         denominator = PyLong_FromLong(1);
1573         py_exponent = PyLong_FromLong(labs((long)exponent));
1574         if (py_exponent == NULL) goto error;
1575         INPLACE_UPDATE(py_exponent,
1576                        long_methods->nb_lshift(denominator, py_exponent));
1577         if (py_exponent == NULL) goto error;
1578         if (exponent > 0) {
1579                 INPLACE_UPDATE(numerator,
1580                                long_methods->nb_multiply(numerator, py_exponent));
1581                 if (numerator == NULL) goto error;
1582         }
1583         else {
1584                 Py_DECREF(denominator);
1585                 denominator = py_exponent;
1586                 py_exponent = NULL;
1587         }
1588
1589         /* Returns ints instead of longs where possible */
1590         INPLACE_UPDATE(numerator, PyNumber_Int(numerator));
1591         if (numerator == NULL) goto error;
1592         INPLACE_UPDATE(denominator, PyNumber_Int(denominator));
1593         if (denominator == NULL) goto error;
1594
1595         result_pair = PyTuple_Pack(2, numerator, denominator);
1596
1597 #undef INPLACE_UPDATE
1598 error:
1599         Py_XDECREF(py_exponent);
1600         Py_XDECREF(denominator);
1601         Py_XDECREF(numerator);
1602         return result_pair;
1603 }
1604
1605 PyDoc_STRVAR(float_as_integer_ratio_doc,
1606 "float.as_integer_ratio() -> (int, int)\n"
1607 "\n"
1608 "Returns a pair of integers, whose ratio is exactly equal to the original\n"
1609 "float and with a positive denominator.\n"
1610 "Raises OverflowError on infinities and a ValueError on NaNs.\n"
1611 "\n"
1612 ">>> (10.0).as_integer_ratio()\n"
1613 "(10, 1)\n"
1614 ">>> (0.0).as_integer_ratio()\n"
1615 "(0, 1)\n"
1616 ">>> (-.25).as_integer_ratio()\n"
1617 "(-1, 4)");
1618
1619
1620 static PyObject *
1621 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
1622
1623 static PyObject *
1624 float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1625 {
1626         PyObject *x = Py_False; /* Integer zero */
1627         static char *kwlist[] = {"x", 0};
1628
1629         if (type != &PyFloat_Type)
1630                 return float_subtype_new(type, args, kwds); /* Wimp out */
1631         if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
1632                 return NULL;
1633         if (PyString_Check(x))
1634                 return PyFloat_FromString(x, NULL);
1635         return PyNumber_Float(x);
1636 }
1637
1638 /* Wimpy, slow approach to tp_new calls for subtypes of float:
1639    first create a regular float from whatever arguments we got,
1640    then allocate a subtype instance and initialize its ob_fval
1641    from the regular float.  The regular float is then thrown away.
1642 */
1643 static PyObject *
1644 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1645 {
1646         PyObject *tmp, *newobj;
1647
1648         assert(PyType_IsSubtype(type, &PyFloat_Type));
1649         tmp = float_new(&PyFloat_Type, args, kwds);
1650         if (tmp == NULL)
1651                 return NULL;
1652         assert(PyFloat_CheckExact(tmp));
1653         newobj = type->tp_alloc(type, 0);
1654         if (newobj == NULL) {
1655                 Py_DECREF(tmp);
1656                 return NULL;
1657         }
1658         ((PyFloatObject *)newobj)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
1659         Py_DECREF(tmp);
1660         return newobj;
1661 }
1662
1663 static PyObject *
1664 float_getnewargs(PyFloatObject *v)
1665 {
1666         return Py_BuildValue("(d)", v->ob_fval);
1667 }
1668
1669 /* this is for the benefit of the pack/unpack routines below */
1670
1671 typedef enum {
1672         unknown_format, ieee_big_endian_format, ieee_little_endian_format
1673 } float_format_type;
1674
1675 static float_format_type double_format, float_format;
1676 static float_format_type detected_double_format, detected_float_format;
1677
1678 static PyObject *
1679 float_getformat(PyTypeObject *v, PyObject* arg)
1680 {
1681         char* s;
1682         float_format_type r;
1683
1684         if (!PyString_Check(arg)) {
1685                 PyErr_Format(PyExc_TypeError,
1686              "__getformat__() argument must be string, not %.500s",
1687                              Py_TYPE(arg)->tp_name);
1688                 return NULL;
1689         }
1690         s = PyString_AS_STRING(arg);
1691         if (strcmp(s, "double") == 0) {
1692                 r = double_format;
1693         }
1694         else if (strcmp(s, "float") == 0) {
1695                 r = float_format;
1696         }
1697         else {
1698                 PyErr_SetString(PyExc_ValueError,
1699                                 "__getformat__() argument 1 must be "
1700                                 "'double' or 'float'");
1701                 return NULL;
1702         }
1703
1704         switch (r) {
1705         case unknown_format:
1706                 return PyString_FromString("unknown");
1707         case ieee_little_endian_format:
1708                 return PyString_FromString("IEEE, little-endian");
1709         case ieee_big_endian_format:
1710                 return PyString_FromString("IEEE, big-endian");
1711         default:
1712                 Py_FatalError("insane float_format or double_format");
1713                 return NULL;
1714         }
1715 }
1716
1717 PyDoc_STRVAR(float_getformat_doc,
1718 "float.__getformat__(typestr) -> string\n"
1719 "\n"
1720 "You probably don't want to use this function.  It exists mainly to be\n"
1721 "used in Python's test suite.\n"
1722 "\n"
1723 "typestr must be 'double' or 'float'.  This function returns whichever of\n"
1724 "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1725 "format of floating point numbers used by the C type named by typestr.");
1726
1727 static PyObject *
1728 float_setformat(PyTypeObject *v, PyObject* args)
1729 {
1730         char* typestr;
1731         char* format;
1732         float_format_type f;
1733         float_format_type detected;
1734         float_format_type *p;
1735
1736         if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format))
1737                 return NULL;
1738
1739         if (strcmp(typestr, "double") == 0) {
1740                 p = &double_format;
1741                 detected = detected_double_format;
1742         }
1743         else if (strcmp(typestr, "float") == 0) {
1744                 p = &float_format;
1745                 detected = detected_float_format;
1746         }
1747         else {
1748                 PyErr_SetString(PyExc_ValueError,
1749                                 "__setformat__() argument 1 must "
1750                                 "be 'double' or 'float'");
1751                 return NULL;
1752         }
1753
1754         if (strcmp(format, "unknown") == 0) {
1755                 f = unknown_format;
1756         }
1757         else if (strcmp(format, "IEEE, little-endian") == 0) {
1758                 f = ieee_little_endian_format;
1759         }
1760         else if (strcmp(format, "IEEE, big-endian") == 0) {
1761                 f = ieee_big_endian_format;
1762         }
1763         else {
1764                 PyErr_SetString(PyExc_ValueError,
1765                                 "__setformat__() argument 2 must be "
1766                                 "'unknown', 'IEEE, little-endian' or "
1767                                 "'IEEE, big-endian'");
1768                 return NULL;
1769
1770         }
1771
1772         if (f != unknown_format && f != detected) {
1773                 PyErr_Format(PyExc_ValueError,
1774                              "can only set %s format to 'unknown' or the "
1775                              "detected platform value", typestr);
1776                 return NULL;
1777         }
1778
1779         *p = f;
1780         Py_RETURN_NONE;
1781 }
1782
1783 PyDoc_STRVAR(float_setformat_doc,
1784 "float.__setformat__(typestr, fmt) -> None\n"
1785 "\n"
1786 "You probably don't want to use this function.  It exists mainly to be\n"
1787 "used in Python's test suite.\n"
1788 "\n"
1789 "typestr must be 'double' or 'float'.  fmt must be one of 'unknown',\n"
1790 "'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1791 "one of the latter two if it appears to match the underlying C reality.\n"
1792 "\n"
1793 "Overrides the automatic determination of C-level floating point type.\n"
1794 "This affects how floats are converted to and from binary strings.");
1795
1796 static PyObject *
1797 float_getzero(PyObject *v, void *closure)
1798 {
1799         return PyFloat_FromDouble(0.0);
1800 }
1801
1802 static PyObject *
1803 float__format__(PyObject *self, PyObject *args)
1804 {
1805         PyObject *format_spec;
1806
1807         if (!PyArg_ParseTuple(args, "O:__format__", &format_spec))
1808                 return NULL;
1809         if (PyBytes_Check(format_spec))
1810                 return _PyFloat_FormatAdvanced(self,
1811                                                PyBytes_AS_STRING(format_spec),
1812                                                PyBytes_GET_SIZE(format_spec));
1813         if (PyUnicode_Check(format_spec)) {
1814                 /* Convert format_spec to a str */
1815                 PyObject *result;
1816                 PyObject *str_spec = PyObject_Str(format_spec);
1817
1818                 if (str_spec == NULL)
1819                         return NULL;
1820
1821                 result = _PyFloat_FormatAdvanced(self,
1822                                                  PyBytes_AS_STRING(str_spec),
1823                                                  PyBytes_GET_SIZE(str_spec));
1824
1825                 Py_DECREF(str_spec);
1826                 return result;
1827         }
1828         PyErr_SetString(PyExc_TypeError, "__format__ requires str or unicode");
1829         return NULL;
1830 }
1831
1832 PyDoc_STRVAR(float__format__doc,
1833 "float.__format__(format_spec) -> string\n"
1834 "\n"
1835 "Formats the float according to format_spec.");
1836
1837
1838 static PyMethodDef float_methods[] = {
1839         {"conjugate",   (PyCFunction)float_float,       METH_NOARGS,
1840          "Returns self, the complex conjugate of any float."},
1841         {"__trunc__",   (PyCFunction)float_trunc, METH_NOARGS,
1842          "Returns the Integral closest to x between 0 and x."},
1843         {"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
1844          float_as_integer_ratio_doc},
1845         {"fromhex", (PyCFunction)float_fromhex,
1846          METH_O|METH_CLASS, float_fromhex_doc},
1847         {"hex", (PyCFunction)float_hex,
1848          METH_NOARGS, float_hex_doc},
1849         {"is_integer",  (PyCFunction)float_is_integer,  METH_NOARGS,
1850          "Returns True if the float is an integer."},
1851 #if 0
1852         {"is_inf",      (PyCFunction)float_is_inf,      METH_NOARGS,
1853          "Returns True if the float is positive or negative infinite."},
1854         {"is_finite",   (PyCFunction)float_is_finite,   METH_NOARGS,
1855          "Returns True if the float is finite, neither infinite nor NaN."},
1856         {"is_nan",      (PyCFunction)float_is_nan,      METH_NOARGS,
1857          "Returns True if the float is not a number (NaN)."},
1858 #endif
1859         {"__getnewargs__",      (PyCFunction)float_getnewargs,  METH_NOARGS},
1860         {"__getformat__",       (PyCFunction)float_getformat,
1861          METH_O|METH_CLASS,             float_getformat_doc},
1862         {"__setformat__",       (PyCFunction)float_setformat,
1863          METH_VARARGS|METH_CLASS,       float_setformat_doc},
1864         {"__format__",          (PyCFunction)float__format__,
1865          METH_VARARGS,                  float__format__doc},
1866         {NULL,          NULL}           /* sentinel */
1867 };
1868
1869 static PyGetSetDef float_getset[] = {
1870     {"real",
1871      (getter)float_float, (setter)NULL,
1872      "the real part of a complex number",
1873      NULL},
1874     {"imag",
1875      (getter)float_getzero, (setter)NULL,
1876      "the imaginary part of a complex number",
1877      NULL},
1878     {NULL}  /* Sentinel */
1879 };
1880
1881 PyDoc_STRVAR(float_doc,
1882 "float(x) -> floating point number\n\
1883 \n\
1884 Convert a string or number to a floating point number, if possible.");
1885
1886
1887 static PyNumberMethods float_as_number = {
1888         float_add,      /*nb_add*/
1889         float_sub,      /*nb_subtract*/
1890         float_mul,      /*nb_multiply*/
1891         float_classic_div, /*nb_divide*/
1892         float_rem,      /*nb_remainder*/
1893         float_divmod,   /*nb_divmod*/
1894         float_pow,      /*nb_power*/
1895         (unaryfunc)float_neg, /*nb_negative*/
1896         (unaryfunc)float_float, /*nb_positive*/
1897         (unaryfunc)float_abs, /*nb_absolute*/
1898         (inquiry)float_nonzero, /*nb_nonzero*/
1899         0,              /*nb_invert*/
1900         0,              /*nb_lshift*/
1901         0,              /*nb_rshift*/
1902         0,              /*nb_and*/
1903         0,              /*nb_xor*/
1904         0,              /*nb_or*/
1905         float_coerce,   /*nb_coerce*/
1906         float_trunc,    /*nb_int*/
1907         float_long,     /*nb_long*/
1908         float_float,    /*nb_float*/
1909         0,              /* nb_oct */
1910         0,              /* nb_hex */
1911         0,              /* nb_inplace_add */
1912         0,              /* nb_inplace_subtract */
1913         0,              /* nb_inplace_multiply */
1914         0,              /* nb_inplace_divide */
1915         0,              /* nb_inplace_remainder */
1916         0,              /* nb_inplace_power */
1917         0,              /* nb_inplace_lshift */
1918         0,              /* nb_inplace_rshift */
1919         0,              /* nb_inplace_and */
1920         0,              /* nb_inplace_xor */
1921         0,              /* nb_inplace_or */
1922         float_floor_div, /* nb_floor_divide */
1923         float_div,      /* nb_true_divide */
1924         0,              /* nb_inplace_floor_divide */
1925         0,              /* nb_inplace_true_divide */
1926 };
1927
1928 PyTypeObject PyFloat_Type = {
1929         PyVarObject_HEAD_INIT(&PyType_Type, 0)
1930         "float",
1931         sizeof(PyFloatObject),
1932         0,
1933         (destructor)float_dealloc,              /* tp_dealloc */
1934         (printfunc)float_print,                 /* tp_print */
1935         0,                                      /* tp_getattr */
1936         0,                                      /* tp_setattr */
1937         0,                                      /* tp_compare */
1938         (reprfunc)float_repr,                   /* tp_repr */
1939         &float_as_number,                       /* tp_as_number */
1940         0,                                      /* tp_as_sequence */
1941         0,                                      /* tp_as_mapping */
1942         (hashfunc)float_hash,                   /* tp_hash */
1943         0,                                      /* tp_call */
1944         (reprfunc)float_str,                    /* tp_str */
1945         PyObject_GenericGetAttr,                /* tp_getattro */
1946         0,                                      /* tp_setattro */
1947         0,                                      /* tp_as_buffer */
1948         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
1949                 Py_TPFLAGS_BASETYPE,            /* tp_flags */
1950         float_doc,                              /* tp_doc */
1951         0,                                      /* tp_traverse */
1952         0,                                      /* tp_clear */
1953         float_richcompare,                      /* tp_richcompare */
1954         0,                                      /* tp_weaklistoffset */
1955         0,                                      /* tp_iter */
1956         0,                                      /* tp_iternext */
1957         float_methods,                          /* tp_methods */
1958         0,                                      /* tp_members */
1959         float_getset,                           /* tp_getset */
1960         0,                                      /* tp_base */
1961         0,                                      /* tp_dict */
1962         0,                                      /* tp_descr_get */
1963         0,                                      /* tp_descr_set */
1964         0,                                      /* tp_dictoffset */
1965         0,                                      /* tp_init */
1966         0,                                      /* tp_alloc */
1967         float_new,                              /* tp_new */
1968 };
1969
1970 void
1971 _PyFloat_Init(void)
1972 {
1973         /* We attempt to determine if this machine is using IEEE
1974            floating point formats by peering at the bits of some
1975            carefully chosen values.  If it looks like we are on an
1976            IEEE platform, the float packing/unpacking routines can
1977            just copy bits, if not they resort to arithmetic & shifts
1978            and masks.  The shifts & masks approach works on all finite
1979            values, but what happens to infinities, NaNs and signed
1980            zeroes on packing is an accident, and attempting to unpack
1981            a NaN or an infinity will raise an exception.
1982
1983            Note that if we're on some whacked-out platform which uses
1984            IEEE formats but isn't strictly little-endian or big-
1985            endian, we will fall back to the portable shifts & masks
1986            method. */
1987
1988 #if SIZEOF_DOUBLE == 8
1989         {
1990                 double x = 9006104071832581.0;
1991                 if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
1992                         detected_double_format = ieee_big_endian_format;
1993                 else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
1994                         detected_double_format = ieee_little_endian_format;
1995                 else
1996                         detected_double_format = unknown_format;
1997         }
1998 #else
1999         detected_double_format = unknown_format;
2000 #endif
2001
2002 #if SIZEOF_FLOAT == 4
2003         {
2004                 float y = 16711938.0;
2005                 if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
2006                         detected_float_format = ieee_big_endian_format;
2007                 else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
2008                         detected_float_format = ieee_little_endian_format;
2009                 else
2010                         detected_float_format = unknown_format;
2011         }
2012 #else
2013         detected_float_format = unknown_format;
2014 #endif
2015
2016         double_format = detected_double_format;
2017         float_format = detected_float_format;
2018
2019         /* Init float info */
2020         if (FloatInfoType.tp_name == 0)
2021                 PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);
2022 }
2023
2024 int
2025 PyFloat_ClearFreeList(void)
2026 {
2027         PyFloatObject *p;
2028         PyFloatBlock *list, *next;
2029         int i;
2030         int u;                  /* remaining unfreed ints per block */
2031         int freelist_size = 0;
2032
2033         list = block_list;
2034         block_list = NULL;
2035         free_list = NULL;
2036         while (list != NULL) {
2037                 u = 0;
2038                 for (i = 0, p = &list->objects[0];
2039                      i < N_FLOATOBJECTS;
2040                      i++, p++) {
2041                         if (PyFloat_CheckExact(p) && Py_REFCNT(p) != 0)
2042                                 u++;
2043                 }
2044                 next = list->next;
2045                 if (u) {
2046                         list->next = block_list;
2047                         block_list = list;
2048                         for (i = 0, p = &list->objects[0];
2049                              i < N_FLOATOBJECTS;
2050                              i++, p++) {
2051                                 if (!PyFloat_CheckExact(p) ||
2052                                     Py_REFCNT(p) == 0) {
2053                                         Py_TYPE(p) = (struct _typeobject *)
2054                                                 free_list;
2055                                         free_list = p;
2056                                 }
2057                         }
2058                 }
2059                 else {
2060                         PyMem_FREE(list);
2061                 }
2062                 freelist_size += u;
2063                 list = next;
2064         }
2065         return freelist_size;
2066 }
2067
2068 void
2069 PyFloat_Fini(void)
2070 {
2071         PyFloatObject *p;
2072         PyFloatBlock *list;
2073         int i;
2074         int u;                  /* total unfreed floats per block */
2075
2076         u = PyFloat_ClearFreeList();
2077
2078         if (!Py_VerboseFlag)
2079                 return;
2080         fprintf(stderr, "# cleanup floats");
2081         if (!u) {
2082                 fprintf(stderr, "\n");
2083         }
2084         else {
2085                 fprintf(stderr,
2086                         ": %d unfreed float%s\n",
2087                         u, u == 1 ? "" : "s");
2088         }
2089         if (Py_VerboseFlag > 1) {
2090                 list = block_list;
2091                 while (list != NULL) {
2092                         for (i = 0, p = &list->objects[0];
2093                              i < N_FLOATOBJECTS;
2094                              i++, p++) {
2095                                 if (PyFloat_CheckExact(p) &&
2096                                     Py_REFCNT(p) != 0) {
2097                                         char buf[100];
2098                                         PyFloat_AsString(buf, p);
2099                                         /* XXX(twouters) cast refcount to
2100                                            long until %zd is universally
2101                                            available
2102                                          */
2103                                         fprintf(stderr,
2104                              "#   <float at %p, refcnt=%ld, val=%s>\n",
2105                                                 p, (long)Py_REFCNT(p), buf);
2106                                 }
2107                         }
2108                         list = list->next;
2109                 }
2110         }
2111 }
2112
2113 /*----------------------------------------------------------------------------
2114  * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
2115  */
2116 int
2117 _PyFloat_Pack4(double x, unsigned char *p, int le)
2118 {
2119         if (float_format == unknown_format) {
2120                 unsigned char sign;
2121                 int e;
2122                 double f;
2123                 unsigned int fbits;
2124                 int incr = 1;
2125
2126                 if (le) {
2127                         p += 3;
2128                         incr = -1;
2129                 }
2130
2131                 if (x < 0) {
2132                         sign = 1;
2133                         x = -x;
2134                 }
2135                 else
2136                         sign = 0;
2137
2138                 f = frexp(x, &e);
2139
2140                 /* Normalize f to be in the range [1.0, 2.0) */
2141                 if (0.5 <= f && f < 1.0) {
2142                         f *= 2.0;
2143                         e--;
2144                 }
2145                 else if (f == 0.0)
2146                         e = 0;
2147                 else {
2148                         PyErr_SetString(PyExc_SystemError,
2149                                         "frexp() result out of range");
2150                         return -1;
2151                 }
2152
2153                 if (e >= 128)
2154                         goto Overflow;
2155                 else if (e < -126) {
2156                         /* Gradual underflow */
2157                         f = ldexp(f, 126 + e);
2158                         e = 0;
2159                 }
2160                 else if (!(e == 0 && f == 0.0)) {
2161                         e += 127;
2162                         f -= 1.0; /* Get rid of leading 1 */
2163                 }
2164
2165                 f *= 8388608.0; /* 2**23 */
2166                 fbits = (unsigned int)(f + 0.5); /* Round */
2167                 assert(fbits <= 8388608);
2168                 if (fbits >> 23) {
2169                         /* The carry propagated out of a string of 23 1 bits. */
2170                         fbits = 0;
2171                         ++e;
2172                         if (e >= 255)
2173                                 goto Overflow;
2174                 }
2175
2176                 /* First byte */
2177                 *p = (sign << 7) | (e >> 1);
2178                 p += incr;
2179
2180                 /* Second byte */
2181                 *p = (char) (((e & 1) << 7) | (fbits >> 16));
2182                 p += incr;
2183
2184                 /* Third byte */
2185                 *p = (fbits >> 8) & 0xFF;
2186                 p += incr;
2187
2188                 /* Fourth byte */
2189                 *p = fbits & 0xFF;
2190
2191                 /* Done */
2192                 return 0;
2193
2194         }
2195         else {
2196                 float y = (float)x;
2197                 const char *s = (char*)&y;
2198                 int i, incr = 1;
2199
2200                 if (Py_IS_INFINITY(y) && !Py_IS_INFINITY(x))
2201                         goto Overflow;
2202
2203                 if ((float_format == ieee_little_endian_format && !le)
2204                     || (float_format == ieee_big_endian_format && le)) {
2205                         p += 3;
2206                         incr = -1;
2207                 }
2208
2209                 for (i = 0; i < 4; i++) {
2210                         *p = *s++;
2211                         p += incr;
2212                 }
2213                 return 0;
2214         }
2215   Overflow:
2216         PyErr_SetString(PyExc_OverflowError,
2217                         "float too large to pack with f format");
2218         return -1;
2219 }
2220
2221 int
2222 _PyFloat_Pack8(double x, unsigned char *p, int le)
2223 {
2224         if (double_format == unknown_format) {
2225                 unsigned char sign;
2226                 int e;
2227                 double f;
2228                 unsigned int fhi, flo;
2229                 int incr = 1;
2230
2231                 if (le) {
2232                         p += 7;
2233                         incr = -1;
2234                 }
2235
2236                 if (x < 0) {
2237                         sign = 1;
2238                         x = -x;
2239                 }
2240                 else
2241                         sign = 0;
2242
2243                 f = frexp(x, &e);
2244
2245                 /* Normalize f to be in the range [1.0, 2.0) */
2246                 if (0.5 <= f && f < 1.0) {
2247                         f *= 2.0;
2248                         e--;
2249                 }
2250                 else if (f == 0.0)
2251                         e = 0;
2252                 else {
2253                         PyErr_SetString(PyExc_SystemError,
2254                                         "frexp() result out of range");
2255                         return -1;
2256                 }
2257
2258                 if (e >= 1024)
2259                         goto Overflow;
2260                 else if (e < -1022) {
2261                         /* Gradual underflow */
2262                         f = ldexp(f, 1022 + e);
2263                         e = 0;
2264                 }
2265                 else if (!(e == 0 && f == 0.0)) {
2266                         e += 1023;
2267                         f -= 1.0; /* Get rid of leading 1 */
2268                 }
2269
2270                 /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
2271                 f *= 268435456.0; /* 2**28 */
2272                 fhi = (unsigned int)f; /* Truncate */
2273                 assert(fhi < 268435456);
2274
2275                 f -= (double)fhi;
2276                 f *= 16777216.0; /* 2**24 */
2277                 flo = (unsigned int)(f + 0.5); /* Round */
2278                 assert(flo <= 16777216);
2279                 if (flo >> 24) {
2280                         /* The carry propagated out of a string of 24 1 bits. */
2281                         flo = 0;
2282                         ++fhi;
2283                         if (fhi >> 28) {
2284                                 /* And it also progagated out of the next 28 bits. */
2285                                 fhi = 0;
2286                                 ++e;
2287                                 if (e >= 2047)
2288                                         goto Overflow;
2289                         }
2290                 }
2291
2292                 /* First byte */
2293                 *p = (sign << 7) | (e >> 4);
2294                 p += incr;
2295
2296                 /* Second byte */
2297                 *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
2298                 p += incr;
2299
2300                 /* Third byte */
2301                 *p = (fhi >> 16) & 0xFF;
2302                 p += incr;
2303
2304                 /* Fourth byte */
2305                 *p = (fhi >> 8) & 0xFF;
2306                 p += incr;
2307
2308                 /* Fifth byte */
2309                 *p = fhi & 0xFF;
2310                 p += incr;
2311
2312                 /* Sixth byte */
2313                 *p = (flo >> 16) & 0xFF;
2314                 p += incr;
2315
2316                 /* Seventh byte */
2317                 *p = (flo >> 8) & 0xFF;
2318                 p += incr;
2319
2320                 /* Eighth byte */
2321                 *p = flo & 0xFF;
2322                 p += incr;
2323
2324                 /* Done */
2325                 return 0;
2326
2327           Overflow:
2328                 PyErr_SetString(PyExc_OverflowError,
2329                                 "float too large to pack with d format");
2330                 return -1;
2331         }
2332         else {
2333                 const char *s = (char*)&x;
2334                 int i, incr = 1;
2335
2336                 if ((double_format == ieee_little_endian_format && !le)
2337                     || (double_format == ieee_big_endian_format && le)) {
2338                         p += 7;
2339                         incr = -1;
2340                 }
2341
2342                 for (i = 0; i < 8; i++) {
2343                         *p = *s++;
2344                         p += incr;
2345                 }
2346                 return 0;
2347         }
2348 }
2349
2350 double
2351 _PyFloat_Unpack4(const unsigned char *p, int le)
2352 {
2353         if (float_format == unknown_format) {
2354                 unsigned char sign;
2355                 int e;
2356                 unsigned int f;
2357                 double x;
2358                 int incr = 1;
2359
2360                 if (le) {
2361                         p += 3;
2362                         incr = -1;
2363                 }
2364
2365                 /* First byte */
2366                 sign = (*p >> 7) & 1;
2367                 e = (*p & 0x7F) << 1;
2368                 p += incr;
2369
2370                 /* Second byte */
2371                 e |= (*p >> 7) & 1;
2372                 f = (*p & 0x7F) << 16;
2373                 p += incr;
2374
2375                 if (e == 255) {
2376                         PyErr_SetString(
2377                                 PyExc_ValueError,
2378                                 "can't unpack IEEE 754 special value "
2379                                 "on non-IEEE platform");
2380                         return -1;
2381                 }
2382
2383                 /* Third byte */
2384                 f |= *p << 8;
2385                 p += incr;
2386
2387                 /* Fourth byte */
2388                 f |= *p;
2389
2390                 x = (double)f / 8388608.0;
2391
2392                 /* XXX This sadly ignores Inf/NaN issues */
2393                 if (e == 0)
2394                         e = -126;
2395                 else {
2396                         x += 1.0;
2397                         e -= 127;
2398                 }
2399                 x = ldexp(x, e);
2400
2401                 if (sign)
2402                         x = -x;
2403
2404                 return x;
2405         }
2406         else {
2407                 float x;
2408
2409                 if ((float_format == ieee_little_endian_format && !le)
2410                     || (float_format == ieee_big_endian_format && le)) {
2411                         char buf[4];
2412                         char *d = &buf[3];
2413                         int i;
2414
2415                         for (i = 0; i < 4; i++) {
2416                                 *d-- = *p++;
2417                         }
2418                         memcpy(&x, buf, 4);
2419                 }
2420                 else {
2421                         memcpy(&x, p, 4);
2422                 }
2423
2424                 return x;
2425         }
2426 }
2427
2428 double
2429 _PyFloat_Unpack8(const unsigned char *p, int le)
2430 {
2431         if (double_format == unknown_format) {
2432                 unsigned char sign;
2433                 int e;
2434                 unsigned int fhi, flo;
2435                 double x;
2436                 int incr = 1;
2437
2438                 if (le) {
2439                         p += 7;
2440                         incr = -1;
2441                 }
2442
2443                 /* First byte */
2444                 sign = (*p >> 7) & 1;
2445                 e = (*p & 0x7F) << 4;
2446
2447                 p += incr;
2448
2449                 /* Second byte */
2450                 e |= (*p >> 4) & 0xF;
2451                 fhi = (*p & 0xF) << 24;
2452                 p += incr;
2453
2454                 if (e == 2047) {
2455                         PyErr_SetString(
2456                                 PyExc_ValueError,
2457                                 "can't unpack IEEE 754 special value "
2458                                 "on non-IEEE platform");
2459                         return -1.0;
2460                 }
2461
2462                 /* Third byte */
2463                 fhi |= *p << 16;
2464                 p += incr;
2465
2466                 /* Fourth byte */
2467                 fhi |= *p  << 8;
2468                 p += incr;
2469
2470                 /* Fifth byte */
2471                 fhi |= *p;
2472                 p += incr;
2473
2474                 /* Sixth byte */
2475                 flo = *p << 16;
2476                 p += incr;
2477
2478                 /* Seventh byte */
2479                 flo |= *p << 8;
2480                 p += incr;
2481
2482                 /* Eighth byte */
2483                 flo |= *p;
2484
2485                 x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
2486                 x /= 268435456.0; /* 2**28 */
2487
2488                 if (e == 0)
2489                         e = -1022;
2490                 else {
2491                         x += 1.0;
2492                         e -= 1023;
2493                 }
2494                 x = ldexp(x, e);
2495
2496                 if (sign)
2497                         x = -x;
2498
2499                 return x;
2500         }
2501         else {
2502                 double x;
2503
2504                 if ((double_format == ieee_little_endian_format && !le)
2505                     || (double_format == ieee_big_endian_format && le)) {
2506                         char buf[8];
2507                         char *d = &buf[7];
2508                         int i;
2509
2510                         for (i = 0; i < 8; i++) {
2511                                 *d-- = *p++;
2512                         }
2513                         memcpy(&x, buf, 8);
2514                 }
2515                 else {
2516                         memcpy(&x, p, 8);
2517                 }
2518
2519                 return x;
2520         }
2521 }