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