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