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