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