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