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