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]
22 class 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
34 static struct _Py_float_state *
35 get_float_state(void)
36 {
37 PyInterpreterState *interp = _PyInterpreterState_GET();
38 return &interp->float_state;
39 }
40 #endif
41
42
43 double
44 PyFloat_GetMax(void)
45 {
46 return DBL_MAX;
47 }
48
49 double
50 PyFloat_GetMin(void)
51 {
52 return DBL_MIN;
53 }
54
55 static PyTypeObject FloatInfoType;
56
57 PyDoc_STRVAR(floatinfo__doc__,
58 "sys.float_info\n\
59 \n\
60 A named tuple holding information about the float type. It contains low level\n\
61 information about the precision and internal representation. Please study\n\
62 your system's :file:`float.h` for more information.");
63
64 static 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
86 static PyStructSequence_Desc floatinfo_desc = {
87 "sys.float_info", /* name */
88 floatinfo__doc__, /* doc */
89 floatinfo_fields, /* fields */
90 11
91 };
92
93 PyObject *
94 PyFloat_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
130 PyObject *
131 PyFloat_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
159 static PyObject *
160 float_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
199 PyObject *
200 PyFloat_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
249 void
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
273 static void
274 float_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
288 double
289 PyFloat_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
361 static int
362 convert_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
381 static PyObject *
382 float_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
413 static PyObject*
414 float_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
592 static Py_hash_t
593 float_hash(PyFloatObject *v)
594 {
595 return _Py_HashDouble((PyObject *)v, v->ob_fval);
596 }
597
598 static PyObject *
599 float_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
608 static PyObject *
609 float_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
618 static PyObject *
619 float_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
628 static PyObject *
629 float_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
643 static PyObject *
644 float_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
671 static 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
709 static PyObject *
710 float_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
724 static PyObject *
725 float_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
743 static PyObject *
744 float_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
865 static PyObject *
866 float_neg(PyFloatObject *v)
867 {
868 return PyFloat_FromDouble(-v->ob_fval);
869 }
870
871 static PyObject *
872 float_abs(PyFloatObject *v)
873 {
874 return PyFloat_FromDouble(fabs(v->ob_fval));
875 }
876
877 static int
878 float_bool(PyFloatObject *v)
879 {
880 return v->ob_fval != 0.0;
881 }
882
883 /*[clinic input]
884 float.is_integer
885
886 Return True if the float is an integer.
887 [clinic start generated code]*/
888
889 static PyObject *
890 float_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]
912 float.__trunc__
913
914 Return the Integral closest to x between 0 and x.
915 [clinic start generated code]*/
916
917 static PyObject *
918 float___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]
925 float.__floor__
926
927 Return the floor as an Integral.
928 [clinic start generated code]*/
929
930 static PyObject *
931 float___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]
939 float.__ceil__
940
941 Return the ceiling as an Integral.
942 [clinic start generated code]*/
943
944 static PyObject *
945 float___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
961 static PyObject *
962 double_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
1019 static PyObject *
1020 double_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]
1069 float.__round__
1070
1071 ndigits as o_ndigits: object = None
1072 /
1073
1074 Return the Integral closest to x, rounding half toward even.
1075
1076 When an argument is passed, work like built-in round(x, ndigits).
1077 [clinic start generated code]*/
1078
1079 static PyObject *
1080 float___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
1124 static PyObject *
1125 float_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]
1135 float.conjugate
1136
1137 Return self, the complex conjugate of any float.
1138 [clinic start generated code]*/
1139
1140 static PyObject *
1141 float_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
1149 static char
1150 char_from_hex(int x)
1151 {
1152 assert(0 <= x && x < 16);
1153 return Py_hexdigits[x];
1154 }
1155
1156 static int
1157 hex_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]
1228 float.hex
1229
1230 Return 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
1238 static PyObject *
1239 float_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
1296 float.fromhex
1297
1298 string: object
1299 /
1300
1301 Create a floating-point number from a hexadecimal string.
1302
1303 >>> float.fromhex('0x1.ffffp10')
1304 2047.984375
1305 >>> float.fromhex('-0x1p-1074')
1306 -5e-324
1307 [clinic start generated code]*/
1308
1309 static PyObject *
1310 float_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]
1551 float.as_integer_ratio
1552
1553 Return integer ratio.
1554
1555 Return a pair of integers, whose ratio is exactly equal to the original float
1556 and with a positive denominator.
1557
1558 Raise 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
1568 static PyObject *
1569 float_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
1632 error:
1633 Py_XDECREF(py_exponent);
1634 Py_XDECREF(denominator);
1635 Py_XDECREF(numerator);
1636 return result_pair;
1637 }
1638
1639 static PyObject *
1640 float_subtype_new(PyTypeObject *type, PyObject *x);
1641
1642 /*[clinic input]
1643 @classmethod
1644 float.__new__ as float_new
1645 x: object(c_default="NULL") = 0
1646 /
1647
1648 Convert a string or number to a floating point number, if possible.
1649 [clinic start generated code]*/
1650
1651 static PyObject *
1652 float_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 */
1677 static PyObject *
1678 float_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
1697 static PyObject *
1698 float_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]
1716 float.__getnewargs__
1717 [clinic start generated code]*/
1718
1719 static PyObject *
1720 float___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
1728 typedef enum {
1729 unknown_format, ieee_big_endian_format, ieee_little_endian_format
1730 } float_format_type;
1731
1732 static float_format_type double_format, float_format;
1733
1734 /*[clinic input]
1735 @classmethod
1736 float.__getformat__
1737
1738 typestr: str
1739 Must be 'double' or 'float'.
1740 /
1741
1742 You probably don't want to use this function.
1743
1744 It exists mainly to be used in Python's test suite.
1745
1746 This function returns whichever of 'unknown', 'IEEE, big-endian' or 'IEEE,
1747 little-endian' best describes the format of floating point numbers used by the
1748 C type named by typestr.
1749 [clinic start generated code]*/
1750
1751 static PyObject *
1752 float___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
1785 static PyObject *
1786 float_getreal(PyObject *v, void *closure)
1787 {
1788 return float_float(v);
1789 }
1790
1791 static PyObject *
1792 float_getimag(PyObject *v, void *closure)
1793 {
1794 return PyFloat_FromDouble(0.0);
1795 }
1796
1797 /*[clinic input]
1798 float.__format__
1799
1800 format_spec: unicode
1801 /
1802
1803 Formats the float according to format_spec.
1804 [clinic start generated code]*/
1805
1806 static PyObject *
1807 float___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
1825 static 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
1841 static 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
1854 static 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
1890 PyTypeObject 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
1933 void
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
1989 PyStatus
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
2010 void
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
2026 void
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
2036 void
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 */
2045 void
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
2068 int
2069 PyFloat_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
2173 int
2174 PyFloat_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
2281 int
2282 PyFloat_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
2411 double
2412 PyFloat_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
2474 double
2475 PyFloat_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
2553 double
2554 PyFloat_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 }