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