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