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