History

add hg and python

cinap_lenrek committed 15 years ago

Commit 4bb9e221 parent 987265c3706 files · +1228772 −-233Browse files

Showing the changes to sys/src/cmd/python/Objects/floatobject.c (1 of 3706 files). Show all

Unified · Split

sys/src/cmd/python/Objects/floatobject.c+1748−0
@@ -0,0 +1,1748 @@
1+
2+/* Float object implementation */
3+
4+/* XXX There should be overflow checks here, but it's hard to check
5+ for any kind of float exception without losing portability. */
6+
7+#include "Python.h"
8+
9+#include <ctype.h>
10+
11+#if !defined(__STDC__)
12+extern double fmod(double, double);
13+extern double pow(double, double);
14+#endif
15+
16+/* Special free list -- see comments for same code in intobject.c. */
17+#define BLOCK_SIZE 1000 /* 1K less typical malloc overhead */
18+#define BHEAD_SIZE 8 /* Enough for a 64-bit pointer */
19+#define N_FLOATOBJECTS ((BLOCK_SIZE - BHEAD_SIZE) / sizeof(PyFloatObject))
20+
21+struct _floatblock {
22+ struct _floatblock *next;
23+ PyFloatObject objects[N_FLOATOBJECTS];
24+};
25+
26+typedef struct _floatblock PyFloatBlock;
27+
28+static PyFloatBlock *block_list = NULL;
29+static PyFloatObject *free_list = NULL;
30+
31+static PyFloatObject *
32+fill_free_list(void)
33+{
34+ PyFloatObject *p, *q;
35+ /* XXX Float blocks escape the object heap. Use PyObject_MALLOC ??? */
36+ p = (PyFloatObject *) PyMem_MALLOC(sizeof(PyFloatBlock));
37+ if (p == NULL)
38+ return (PyFloatObject *) PyErr_NoMemory();
39+ ((PyFloatBlock *)p)->next = block_list;
40+ block_list = (PyFloatBlock *)p;
41+ p = &((PyFloatBlock *)p)->objects[0];
42+ q = p + N_FLOATOBJECTS;
43+ while (--q > p)
44+ q->ob_type = (struct _typeobject *)(q-1);
45+ q->ob_type = NULL;
46+ return p + N_FLOATOBJECTS - 1;
47+}
48+
49+PyObject *
50+PyFloat_FromDouble(double fval)
51+{
52+ register PyFloatObject *op;
53+ if (free_list == NULL) {
54+ if ((free_list = fill_free_list()) == NULL)
55+ return NULL;
56+ }
57+ /* Inline PyObject_New */
58+ op = free_list;
59+ free_list = (PyFloatObject *)op->ob_type;
60+ PyObject_INIT(op, &PyFloat_Type);
61+ op->ob_fval = fval;
62+ return (PyObject *) op;
63+}
64+
65+/**************************************************************************
66+RED_FLAG 22-Sep-2000 tim
67+PyFloat_FromString's pend argument is braindead. Prior to this RED_FLAG,
68+
69+1. If v was a regular string, *pend was set to point to its terminating
70+ null byte. That's useless (the caller can find that without any
71+ help from this function!).
72+
73+2. If v was a Unicode string, or an object convertible to a character
74+ buffer, *pend was set to point into stack trash (the auto temp
75+ vector holding the character buffer). That was downright dangerous.
76+
77+Since we can't change the interface of a public API function, pend is
78+still supported but now *officially* useless: if pend is not NULL,
79+*pend is set to NULL.
80+**************************************************************************/
81+PyObject *
82+PyFloat_FromString(PyObject *v, char **pend)
83+{
84+ const char *s, *last, *end;
85+ double x;
86+ char buffer[256]; /* for errors */
87+#ifdef Py_USING_UNICODE
88+ char s_buffer[256]; /* for objects convertible to a char buffer */
89+#endif
90+ Py_ssize_t len;
91+
92+ if (pend)
93+ *pend = NULL;
94+ if (PyString_Check(v)) {
95+ s = PyString_AS_STRING(v);
96+ len = PyString_GET_SIZE(v);
97+ }
98+#ifdef Py_USING_UNICODE
99+ else if (PyUnicode_Check(v)) {
100+ if (PyUnicode_GET_SIZE(v) >= (Py_ssize_t)sizeof(s_buffer)) {
101+ PyErr_SetString(PyExc_ValueError,
102+ "Unicode float() literal too long to convert");
103+ return NULL;
104+ }
105+ if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
106+ PyUnicode_GET_SIZE(v),
107+ s_buffer,
108+ NULL))
109+ return NULL;
110+ s = s_buffer;
111+ len = strlen(s);
112+ }
113+#endif
114+ else if (PyObject_AsCharBuffer(v, &s, &len)) {
115+ PyErr_SetString(PyExc_TypeError,
116+ "float() argument must be a string or a number");
117+ return NULL;
118+ }
119+
120+ last = s + len;
121+ while (*s && isspace(Py_CHARMASK(*s)))
122+ s++;
123+ if (*s == '\0') {
124+ PyErr_SetString(PyExc_ValueError, "empty string for float()");
125+ return NULL;
126+ }
127+ /* We don't care about overflow or underflow. If the platform supports
128+ * them, infinities and signed zeroes (on underflow) are fine.
129+ * However, strtod can return 0 for denormalized numbers, where atof
130+ * does not. So (alas!) we special-case a zero result. Note that
131+ * whether strtod sets errno on underflow is not defined, so we can't
132+ * key off errno.
133+ */
134+ PyFPE_START_PROTECT("strtod", return NULL)
135+ x = PyOS_ascii_strtod(s, (char **)&end);
136+ PyFPE_END_PROTECT(x)
137+ errno = 0;
138+ /* Believe it or not, Solaris 2.6 can move end *beyond* the null
139+ byte at the end of the string, when the input is inf(inity). */
140+ if (end > last)
141+ end = last;
142+ if (end == s) {
143+ PyOS_snprintf(buffer, sizeof(buffer),
144+ "invalid literal for float(): %.200s", s);
145+ PyErr_SetString(PyExc_ValueError, buffer);
146+ return NULL;
147+ }
148+ /* Since end != s, the platform made *some* kind of sense out
149+ of the input. Trust it. */
150+ while (*end && isspace(Py_CHARMASK(*end)))
151+ end++;
152+ if (*end != '\0') {
153+ PyOS_snprintf(buffer, sizeof(buffer),
154+ "invalid literal for float(): %.200s", s);
155+ PyErr_SetString(PyExc_ValueError, buffer);
156+ return NULL;
157+ }
158+ else if (end != last) {
159+ PyErr_SetString(PyExc_ValueError,
160+ "null byte in argument for float()");
161+ return NULL;
162+ }
163+ if (x == 0.0) {
164+ /* See above -- may have been strtod being anal
165+ about denorms. */
166+ PyFPE_START_PROTECT("atof", return NULL)
167+ x = PyOS_ascii_atof(s);
168+ PyFPE_END_PROTECT(x)
169+ errno = 0; /* whether atof ever set errno is undefined */
170+ }
171+ return PyFloat_FromDouble(x);
172+}
173+
174+static void
175+float_dealloc(PyFloatObject *op)
176+{
177+ if (PyFloat_CheckExact(op)) {
178+ op->ob_type = (struct _typeobject *)free_list;
179+ free_list = op;
180+ }
181+ else
182+ op->ob_type->tp_free((PyObject *)op);
183+}
184+
185+double
186+PyFloat_AsDouble(PyObject *op)
187+{
188+ PyNumberMethods *nb;
189+ PyFloatObject *fo;
190+ double val;
191+
192+ if (op && PyFloat_Check(op))
193+ return PyFloat_AS_DOUBLE((PyFloatObject*) op);
194+
195+ if (op == NULL) {
196+ PyErr_BadArgument();
197+ return -1;
198+ }
199+
200+ if ((nb = op->ob_type->tp_as_number) == NULL || nb->nb_float == NULL) {
201+ PyErr_SetString(PyExc_TypeError, "a float is required");
202+ return -1;
203+ }
204+
205+ fo = (PyFloatObject*) (*nb->nb_float) (op);
206+ if (fo == NULL)
207+ return -1;
208+ if (!PyFloat_Check(fo)) {
209+ PyErr_SetString(PyExc_TypeError,
210+ "nb_float should return float object");
211+ return -1;
212+ }
213+
214+ val = PyFloat_AS_DOUBLE(fo);
215+ Py_DECREF(fo);
216+
217+ return val;
218+}
219+
220+/* Methods */
221+
222+static void
223+format_float(char *buf, size_t buflen, PyFloatObject *v, int precision)
224+{
225+ register char *cp;
226+ char format[32];
227+ /* Subroutine for float_repr and float_print.
228+ We want float numbers to be recognizable as such,
229+ i.e., they should contain a decimal point or an exponent.
230+ However, %g may print the number as an integer;
231+ in such cases, we append ".0" to the string. */
232+
233+ assert(PyFloat_Check(v));
234+ PyOS_snprintf(format, 32, "%%.%ig", precision);
235+ PyOS_ascii_formatd(buf, buflen, format, v->ob_fval);
236+ cp = buf;
237+ if (*cp == '-')
238+ cp++;
239+ for (; *cp != '\0'; cp++) {
240+ /* Any non-digit means it's not an integer;
241+ this takes care of NAN and INF as well. */
242+ if (!isdigit(Py_CHARMASK(*cp)))
243+ break;
244+ }
245+ if (*cp == '\0') {
246+ *cp++ = '.';
247+ *cp++ = '0';
248+ *cp++ = '\0';
249+ }
250+}
251+
252+/* XXX PyFloat_AsStringEx should not be a public API function (for one
253+ XXX thing, its signature passes a buffer without a length; for another,
254+ XXX it isn't useful outside this file).
255+*/
256+void
257+PyFloat_AsStringEx(char *buf, PyFloatObject *v, int precision)
258+{
259+ format_float(buf, 100, v, precision);
260+}
261+
262+/* Macro and helper that convert PyObject obj to a C double and store
263+ the value in dbl; this replaces the functionality of the coercion
264+ slot function. If conversion to double raises an exception, obj is
265+ set to NULL, and the function invoking this macro returns NULL. If
266+ obj is not of float, int or long type, Py_NotImplemented is incref'ed,
267+ stored in obj, and returned from the function invoking this macro.
268+*/
269+#define CONVERT_TO_DOUBLE(obj, dbl) \
270+ if (PyFloat_Check(obj)) \
271+ dbl = PyFloat_AS_DOUBLE(obj); \
272+ else if (convert_to_double(&(obj), &(dbl)) < 0) \
273+ return obj;
274+
275+static int
276+convert_to_double(PyObject **v, double *dbl)
277+{
278+ register PyObject *obj = *v;
279+
280+ if (PyInt_Check(obj)) {
281+ *dbl = (double)PyInt_AS_LONG(obj);
282+ }
283+ else if (PyLong_Check(obj)) {
284+ *dbl = PyLong_AsDouble(obj);
285+ if (*dbl == -1.0 && PyErr_Occurred()) {
286+ *v = NULL;
287+ return -1;
288+ }
289+ }
290+ else {
291+ Py_INCREF(Py_NotImplemented);
292+ *v = Py_NotImplemented;
293+ return -1;
294+ }
295+ return 0;
296+}
297+
298+/* Precisions used by repr() and str(), respectively.
299+
300+ The repr() precision (17 significant decimal digits) is the minimal number
301+ that is guaranteed to have enough precision so that if the number is read
302+ back in the exact same binary value is recreated. This is true for IEEE
303+ floating point by design, and also happens to work for all other modern
304+ hardware.
305+
306+ The str() precision is chosen so that in most cases, the rounding noise
307+ created by various operations is suppressed, while giving plenty of
308+ precision for practical use.
309+
310+*/
311+
312+#define PREC_REPR 17
313+#define PREC_STR 12
314+
315+/* XXX PyFloat_AsString and PyFloat_AsReprString should be deprecated:
316+ XXX they pass a char buffer without passing a length.
317+*/
318+void
319+PyFloat_AsString(char *buf, PyFloatObject *v)
320+{
321+ format_float(buf, 100, v, PREC_STR);
322+}
323+
324+void
325+PyFloat_AsReprString(char *buf, PyFloatObject *v)
326+{
327+ format_float(buf, 100, v, PREC_REPR);
328+}
329+
330+/* ARGSUSED */
331+static int
332+float_print(PyFloatObject *v, FILE *fp, int flags)
333+{
334+ char buf[100];
335+ format_float(buf, sizeof(buf), v,
336+ (flags & Py_PRINT_RAW) ? PREC_STR : PREC_REPR);
337+ fputs(buf, fp);
338+ return 0;
339+}
340+
341+static PyObject *
342+float_repr(PyFloatObject *v)
343+{
344+ char buf[100];
345+ format_float(buf, sizeof(buf), v, PREC_REPR);
346+ return PyString_FromString(buf);
347+}
348+
349+static PyObject *
350+float_str(PyFloatObject *v)
351+{
352+ char buf[100];
353+ format_float(buf, sizeof(buf), v, PREC_STR);
354+ return PyString_FromString(buf);
355+}
356+
357+/* Comparison is pretty much a nightmare. When comparing float to float,
358+ * we do it as straightforwardly (and long-windedly) as conceivable, so
359+ * that, e.g., Python x == y delivers the same result as the platform
360+ * C x == y when x and/or y is a NaN.
361+ * When mixing float with an integer type, there's no good *uniform* approach.
362+ * Converting the double to an integer obviously doesn't work, since we
363+ * may lose info from fractional bits. Converting the integer to a double
364+ * also has two failure modes: (1) a long int may trigger overflow (too
365+ * large to fit in the dynamic range of a C double); (2) even a C long may have
366+ * more bits than fit in a C double (e.g., on a a 64-bit box long may have
367+ * 63 bits of precision, but a C double probably has only 53), and then
368+ * we can falsely claim equality when low-order integer bits are lost by
369+ * coercion to double. So this part is painful too.
370+ */
371+
372+static PyObject*
373+float_richcompare(PyObject *v, PyObject *w, int op)
374+{
375+ double i, j;
376+ int r = 0;
377+
378+ assert(PyFloat_Check(v));
379+ i = PyFloat_AS_DOUBLE(v);
380+
381+ /* Switch on the type of w. Set i and j to doubles to be compared,
382+ * and op to the richcomp to use.
383+ */
384+ if (PyFloat_Check(w))
385+ j = PyFloat_AS_DOUBLE(w);
386+
387+ else if (!Py_IS_FINITE(i)) {
388+ if (PyInt_Check(w) || PyLong_Check(w))
389+ /* If i is an infinity, its magnitude exceeds any
390+ * finite integer, so it doesn't matter which int we
391+ * compare i with. If i is a NaN, similarly.
392+ */
393+ j = 0.0;
394+ else
395+ goto Unimplemented;
396+ }
397+
398+ else if (PyInt_Check(w)) {
399+ long jj = PyInt_AS_LONG(w);
400+ /* In the worst realistic case I can imagine, C double is a
401+ * Cray single with 48 bits of precision, and long has 64
402+ * bits.
403+ */
404+#if SIZEOF_LONG > 6
405+ unsigned long abs = (unsigned long)(jj < 0 ? -jj : jj);
406+ if (abs >> 48) {
407+ /* Needs more than 48 bits. Make it take the
408+ * PyLong path.
409+ */
410+ PyObject *result;
411+ PyObject *ww = PyLong_FromLong(jj);
412+
413+ if (ww == NULL)
414+ return NULL;
415+ result = float_richcompare(v, ww, op);
416+ Py_DECREF(ww);
417+ return result;
418+ }
419+#endif
420+ j = (double)jj;
421+ assert((long)j == jj);
422+ }
423+
424+ else if (PyLong_Check(w)) {
425+ int vsign = i == 0.0 ? 0 : i < 0.0 ? -1 : 1;
426+ int wsign = _PyLong_Sign(w);
427+ size_t nbits;
428+ int exponent;
429+
430+ if (vsign != wsign) {
431+ /* Magnitudes are irrelevant -- the signs alone
432+ * determine the outcome.
433+ */
434+ i = (double)vsign;
435+ j = (double)wsign;
436+ goto Compare;
437+ }
438+ /* The signs are the same. */
439+ /* Convert w to a double if it fits. In particular, 0 fits. */
440+ nbits = _PyLong_NumBits(w);
441+ if (nbits == (size_t)-1 && PyErr_Occurred()) {
442+ /* This long is so large that size_t isn't big enough
443+ * to hold the # of bits. Replace with little doubles
444+ * that give the same outcome -- w is so large that
445+ * its magnitude must exceed the magnitude of any
446+ * finite float.
447+ */
448+ PyErr_Clear();
449+ i = (double)vsign;
450+ assert(wsign != 0);
451+ j = wsign * 2.0;
452+ goto Compare;
453+ }
454+ if (nbits <= 48) {
455+ j = PyLong_AsDouble(w);
456+ /* It's impossible that <= 48 bits overflowed. */
457+ assert(j != -1.0 || ! PyErr_Occurred());
458+ goto Compare;
459+ }
460+ assert(wsign != 0); /* else nbits was 0 */
461+ assert(vsign != 0); /* if vsign were 0, then since wsign is
462+ * not 0, we would have taken the
463+ * vsign != wsign branch at the start */
464+ /* We want to work with non-negative numbers. */
465+ if (vsign < 0) {
466+ /* "Multiply both sides" by -1; this also swaps the
467+ * comparator.
468+ */
469+ i = -i;
470+ op = _Py_SwappedOp[op];
471+ }
472+ assert(i > 0.0);
473+ (void) frexp(i, &exponent);
474+ /* exponent is the # of bits in v before the radix point;
475+ * we know that nbits (the # of bits in w) > 48 at this point
476+ */
477+ if (exponent < 0 || (size_t)exponent < nbits) {
478+ i = 1.0;
479+ j = 2.0;
480+ goto Compare;
481+ }
482+ if ((size_t)exponent > nbits) {
483+ i = 2.0;
484+ j = 1.0;
485+ goto Compare;
486+ }
487+ /* v and w have the same number of bits before the radix
488+ * point. Construct two longs that have the same comparison
489+ * outcome.
490+ */
491+ {
492+ double fracpart;
493+ double intpart;
494+ PyObject *result = NULL;
495+ PyObject *one = NULL;
496+ PyObject *vv = NULL;
497+ PyObject *ww = w;
498+
499+ if (wsign < 0) {
500+ ww = PyNumber_Negative(w);
501+ if (ww == NULL)
502+ goto Error;
503+ }
504+ else
505+ Py_INCREF(ww);
506+
507+ fracpart = modf(i, &intpart);
508+ vv = PyLong_FromDouble(intpart);
509+ if (vv == NULL)
510+ goto Error;
511+
512+ if (fracpart != 0.0) {
513+ /* Shift left, and or a 1 bit into vv
514+ * to represent the lost fraction.
515+ */
516+ PyObject *temp;
517+
518+ one = PyInt_FromLong(1);
519+ if (one == NULL)
520+ goto Error;
521+
522+ temp = PyNumber_Lshift(ww, one);
523+ if (temp == NULL)
524+ goto Error;
525+ Py_DECREF(ww);
526+ ww = temp;
527+
528+ temp = PyNumber_Lshift(vv, one);
529+ if (temp == NULL)
530+ goto Error;
531+ Py_DECREF(vv);
532+ vv = temp;
533+
534+ temp = PyNumber_Or(vv, one);
535+ if (temp == NULL)
536+ goto Error;
537+ Py_DECREF(vv);
538+ vv = temp;
539+ }
540+
541+ r = PyObject_RichCompareBool(vv, ww, op);
542+ if (r < 0)
543+ goto Error;
544+ result = PyBool_FromLong(r);
545+ Error:
546+ Py_XDECREF(vv);
547+ Py_XDECREF(ww);
548+ Py_XDECREF(one);
549+ return result;
550+ }
551+ } /* else if (PyLong_Check(w)) */
552+
553+ else /* w isn't float, int, or long */
554+ goto Unimplemented;
555+
556+ Compare:
557+ PyFPE_START_PROTECT("richcompare", return NULL)
558+ switch (op) {
559+ case Py_EQ:
560+ r = i == j;
561+ break;
562+ case Py_NE:
563+ r = i != j;
564+ break;
565+ case Py_LE:
566+ r = i <= j;
567+ break;
568+ case Py_GE:
569+ r = i >= j;
570+ break;
571+ case Py_LT:
572+ r = i < j;
573+ break;
574+ case Py_GT:
575+ r = i > j;
576+ break;
577+ }
578+ PyFPE_END_PROTECT(r)
579+ return PyBool_FromLong(r);
580+
581+ Unimplemented:
582+ Py_INCREF(Py_NotImplemented);
583+ return Py_NotImplemented;
584+}
585+
586+static long
587+float_hash(PyFloatObject *v)
588+{
589+ return _Py_HashDouble(v->ob_fval);
590+}
591+
592+static PyObject *
593+float_add(PyObject *v, PyObject *w)
594+{
595+ double a,b;
596+ CONVERT_TO_DOUBLE(v, a);
597+ CONVERT_TO_DOUBLE(w, b);
598+ PyFPE_START_PROTECT("add", return 0)
599+ a = a + b;
600+ PyFPE_END_PROTECT(a)
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+ PyFPE_START_PROTECT("subtract", return 0)
611+ a = a - b;
612+ PyFPE_END_PROTECT(a)
613+ return PyFloat_FromDouble(a);
614+}
615+
616+static PyObject *
617+float_mul(PyObject *v, PyObject *w)
618+{
619+ double a,b;
620+ CONVERT_TO_DOUBLE(v, a);
621+ CONVERT_TO_DOUBLE(w, b);
622+ PyFPE_START_PROTECT("multiply", return 0)
623+ a = a * b;
624+ PyFPE_END_PROTECT(a)
625+ return PyFloat_FromDouble(a);
626+}
627+
628+static PyObject *
629+float_div(PyObject *v, PyObject *w)
630+{
631+ double a,b;
632+ CONVERT_TO_DOUBLE(v, a);
633+ CONVERT_TO_DOUBLE(w, b);
634+ if (b == 0.0) {
635+ PyErr_SetString(PyExc_ZeroDivisionError, "float division");
636+ return NULL;
637+ }
638+ PyFPE_START_PROTECT("divide", return 0)
639+ a = a / b;
640+ PyFPE_END_PROTECT(a)
641+ return PyFloat_FromDouble(a);
642+}
643+
644+static PyObject *
645+float_classic_div(PyObject *v, PyObject *w)
646+{
647+ double a,b;
648+ CONVERT_TO_DOUBLE(v, a);
649+ CONVERT_TO_DOUBLE(w, b);
650+ if (Py_DivisionWarningFlag >= 2 &&
651+ PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)
652+ return NULL;
653+ if (b == 0.0) {
654+ PyErr_SetString(PyExc_ZeroDivisionError, "float division");
655+ return NULL;
656+ }
657+ PyFPE_START_PROTECT("divide", return 0)
658+ a = a / b;
659+ PyFPE_END_PROTECT(a)
660+ return PyFloat_FromDouble(a);
661+}
662+
663+static PyObject *
664+float_rem(PyObject *v, PyObject *w)
665+{
666+ double vx, wx;
667+ double mod;
668+ CONVERT_TO_DOUBLE(v, vx);
669+ CONVERT_TO_DOUBLE(w, wx);
670+ if (wx == 0.0) {
671+ PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
672+ return NULL;
673+ }
674+ PyFPE_START_PROTECT("modulo", return 0)
675+ mod = fmod(vx, wx);
676+ /* note: checking mod*wx < 0 is incorrect -- underflows to
677+ 0 if wx < sqrt(smallest nonzero double) */
678+ if (mod && ((wx < 0) != (mod < 0))) {
679+ mod += wx;
680+ }
681+ PyFPE_END_PROTECT(mod)
682+ return PyFloat_FromDouble(mod);
683+}
684+
685+static PyObject *
686+float_divmod(PyObject *v, PyObject *w)
687+{
688+ double vx, wx;
689+ double div, mod, floordiv;
690+ CONVERT_TO_DOUBLE(v, vx);
691+ CONVERT_TO_DOUBLE(w, wx);
692+ if (wx == 0.0) {
693+ PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");
694+ return NULL;
695+ }
696+ PyFPE_START_PROTECT("divmod", return 0)
697+ mod = fmod(vx, wx);
698+ /* fmod is typically exact, so vx-mod is *mathematically* an
699+ exact multiple of wx. But this is fp arithmetic, and fp
700+ vx - mod is an approximation; the result is that div may
701+ not be an exact integral value after the division, although
702+ it will always be very close to one.
703+ */
704+ div = (vx - mod) / wx;
705+ if (mod) {
706+ /* ensure the remainder has the same sign as the denominator */
707+ if ((wx < 0) != (mod < 0)) {
708+ mod += wx;
709+ div -= 1.0;
710+ }
711+ }
712+ else {
713+ /* the remainder is zero, and in the presence of signed zeroes
714+ fmod returns different results across platforms; ensure
715+ it has the same sign as the denominator; we'd like to do
716+ "mod = wx * 0.0", but that may get optimized away */
717+ mod *= mod; /* hide "mod = +0" from optimizer */
718+ if (wx < 0.0)
719+ mod = -mod;
720+ }
721+ /* snap quotient to nearest integral value */
722+ if (div) {
723+ floordiv = floor(div);
724+ if (div - floordiv > 0.5)
725+ floordiv += 1.0;
726+ }
727+ else {
728+ /* div is zero - get the same sign as the true quotient */
729+ div *= div; /* hide "div = +0" from optimizers */
730+ floordiv = div * vx / wx; /* zero w/ sign of vx/wx */
731+ }
732+ PyFPE_END_PROTECT(floordiv)
733+ return Py_BuildValue("(dd)", floordiv, mod);
734+}
735+
736+static PyObject *
737+float_floor_div(PyObject *v, PyObject *w)
738+{
739+ PyObject *t, *r;
740+
741+ t = float_divmod(v, w);
742+ if (t == NULL || t == Py_NotImplemented)
743+ return t;
744+ assert(PyTuple_CheckExact(t));
745+ r = PyTuple_GET_ITEM(t, 0);
746+ Py_INCREF(r);
747+ Py_DECREF(t);
748+ return r;
749+}
750+
751+static PyObject *
752+float_pow(PyObject *v, PyObject *w, PyObject *z)
753+{
754+ double iv, iw, ix;
755+
756+ if ((PyObject *)z != Py_None) {
757+ PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "
758+ "allowed unless all arguments are integers");
759+ return NULL;
760+ }
761+
762+ CONVERT_TO_DOUBLE(v, iv);
763+ CONVERT_TO_DOUBLE(w, iw);
764+
765+ /* Sort out special cases here instead of relying on pow() */
766+ if (iw == 0) { /* v**0 is 1, even 0**0 */
767+ PyFPE_START_PROTECT("pow", return NULL)
768+ if ((PyObject *)z != Py_None) {
769+ double iz;
770+ CONVERT_TO_DOUBLE(z, iz);
771+ ix = fmod(1.0, iz);
772+ if (ix != 0 && iz < 0)
773+ ix += iz;
774+ }
775+ else
776+ ix = 1.0;
777+ PyFPE_END_PROTECT(ix)
778+ return PyFloat_FromDouble(ix);
779+ }
780+ if (iv == 0.0) { /* 0**w is error if w<0, else 1 */
781+ if (iw < 0.0) {
782+ PyErr_SetString(PyExc_ZeroDivisionError,
783+ "0.0 cannot be raised to a negative power");
784+ return NULL;
785+ }
786+ return PyFloat_FromDouble(0.0);
787+ }
788+ if (iv < 0.0) {
789+ /* Whether this is an error is a mess, and bumps into libm
790+ * bugs so we have to figure it out ourselves.
791+ */
792+ if (iw != floor(iw)) {
793+ PyErr_SetString(PyExc_ValueError, "negative number "
794+ "cannot be raised to a fractional power");
795+ return NULL;
796+ }
797+ /* iw is an exact integer, albeit perhaps a very large one.
798+ * -1 raised to an exact integer should never be exceptional.
799+ * Alas, some libms (chiefly glibc as of early 2003) return
800+ * NaN and set EDOM on pow(-1, large_int) if the int doesn't
801+ * happen to be representable in a *C* integer. That's a
802+ * bug; we let that slide in math.pow() (which currently
803+ * reflects all platform accidents), but not for Python's **.
804+ */
805+ if (iv == -1.0 && Py_IS_FINITE(iw)) {
806+ /* Return 1 if iw is even, -1 if iw is odd; there's
807+ * no guarantee that any C integral type is big
808+ * enough to hold iw, so we have to check this
809+ * indirectly.
810+ */
811+ ix = floor(iw * 0.5) * 2.0;
812+ return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0);
813+ }
814+ /* Else iv != -1.0, and overflow or underflow are possible.
815+ * Unless we're to write pow() ourselves, we have to trust
816+ * the platform to do this correctly.
817+ */
818+ }
819+ errno = 0;
820+ PyFPE_START_PROTECT("pow", return NULL)
821+ ix = pow(iv, iw);
822+ PyFPE_END_PROTECT(ix)
823+ Py_ADJUST_ERANGE1(ix);
824+ if (errno != 0) {
825+ /* We don't expect any errno value other than ERANGE, but
826+ * the range of libm bugs appears unbounded.
827+ */
828+ PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
829+ PyExc_ValueError);
830+ return NULL;
831+ }
832+ return PyFloat_FromDouble(ix);
833+}
834+
835+static PyObject *
836+float_neg(PyFloatObject *v)
837+{
838+ return PyFloat_FromDouble(-v->ob_fval);
839+}
840+
841+static PyObject *
842+float_pos(PyFloatObject *v)
843+{
844+ if (PyFloat_CheckExact(v)) {
845+ Py_INCREF(v);
846+ return (PyObject *)v;
847+ }
848+ else
849+ return PyFloat_FromDouble(v->ob_fval);
850+}
851+
852+static PyObject *
853+float_abs(PyFloatObject *v)
854+{
855+ return PyFloat_FromDouble(fabs(v->ob_fval));
856+}
857+
858+static int
859+float_nonzero(PyFloatObject *v)
860+{
861+ return v->ob_fval != 0.0;
862+}
863+
864+static int
865+float_coerce(PyObject **pv, PyObject **pw)
866+{
867+ if (PyInt_Check(*pw)) {
868+ long x = PyInt_AsLong(*pw);
869+ *pw = PyFloat_FromDouble((double)x);
870+ Py_INCREF(*pv);
871+ return 0;
872+ }
873+ else if (PyLong_Check(*pw)) {
874+ double x = PyLong_AsDouble(*pw);
875+ if (x == -1.0 && PyErr_Occurred())
876+ return -1;
877+ *pw = PyFloat_FromDouble(x);
878+ Py_INCREF(*pv);
879+ return 0;
880+ }
881+ else if (PyFloat_Check(*pw)) {
882+ Py_INCREF(*pv);
883+ Py_INCREF(*pw);
884+ return 0;
885+ }
886+ return 1; /* Can't do it */
887+}
888+
889+static PyObject *
890+float_long(PyObject *v)
891+{
892+ double x = PyFloat_AsDouble(v);
893+ return PyLong_FromDouble(x);
894+}
895+
896+static PyObject *
897+float_int(PyObject *v)
898+{
899+ double x = PyFloat_AsDouble(v);
900+ double wholepart; /* integral portion of x, rounded toward 0 */
901+
902+ (void)modf(x, &wholepart);
903+ /* Try to get out cheap if this fits in a Python int. The attempt
904+ * to cast to long must be protected, as C doesn't define what
905+ * happens if the double is too big to fit in a long. Some rare
906+ * systems raise an exception then (RISCOS was mentioned as one,
907+ * and someone using a non-default option on Sun also bumped into
908+ * that). Note that checking for >= and <= LONG_{MIN,MAX} would
909+ * still be vulnerable: if a long has more bits of precision than
910+ * a double, casting MIN/MAX to double may yield an approximation,
911+ * and if that's rounded up, then, e.g., wholepart=LONG_MAX+1 would
912+ * yield true from the C expression wholepart<=LONG_MAX, despite
913+ * that wholepart is actually greater than LONG_MAX.
914+ */
915+ if (LONG_MIN < wholepart && wholepart < LONG_MAX) {
916+ const long aslong = (long)wholepart;
917+ return PyInt_FromLong(aslong);
918+ }
919+ return PyLong_FromDouble(wholepart);
920+}
921+
922+static PyObject *
923+float_float(PyObject *v)
924+{
925+ if (PyFloat_CheckExact(v))
926+ Py_INCREF(v);
927+ else
928+ v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
929+ return v;
930+}
931+
932+
933+static PyObject *
934+float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
935+
936+static PyObject *
937+float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
938+{
939+ PyObject *x = Py_False; /* Integer zero */
940+ static char *kwlist[] = {"x", 0};
941+
942+ if (type != &PyFloat_Type)
943+ return float_subtype_new(type, args, kwds); /* Wimp out */
944+ if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
945+ return NULL;
946+ if (PyString_Check(x))
947+ return PyFloat_FromString(x, NULL);
948+ return PyNumber_Float(x);
949+}
950+
951+/* Wimpy, slow approach to tp_new calls for subtypes of float:
952+ first create a regular float from whatever arguments we got,
953+ then allocate a subtype instance and initialize its ob_fval
954+ from the regular float. The regular float is then thrown away.
955+*/
956+static PyObject *
957+float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
958+{
959+ PyObject *tmp, *newobj;
960+
961+ assert(PyType_IsSubtype(type, &PyFloat_Type));
962+ tmp = float_new(&PyFloat_Type, args, kwds);
963+ if (tmp == NULL)
964+ return NULL;
965+ assert(PyFloat_CheckExact(tmp));
966+ newobj = type->tp_alloc(type, 0);
967+ if (newobj == NULL) {
968+ Py_DECREF(tmp);
969+ return NULL;
970+ }
971+ ((PyFloatObject *)newobj)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
972+ Py_DECREF(tmp);
973+ return newobj;
974+}
975+
976+static PyObject *
977+float_getnewargs(PyFloatObject *v)
978+{
979+ return Py_BuildValue("(d)", v->ob_fval);
980+}
981+
982+/* this is for the benefit of the pack/unpack routines below */
983+
984+typedef enum {
985+ unknown_format, ieee_big_endian_format, ieee_little_endian_format
986+} float_format_type;
987+
988+static float_format_type double_format, float_format;
989+static float_format_type detected_double_format, detected_float_format;
990+
991+static PyObject *
992+float_getformat(PyTypeObject *v, PyObject* arg)
993+{
994+ char* s;
995+ float_format_type r;
996+
997+ if (!PyString_Check(arg)) {
998+ PyErr_Format(PyExc_TypeError,
999+ "__getformat__() argument must be string, not %.500s",
1000+ arg->ob_type->tp_name);
1001+ return NULL;
1002+ }
1003+ s = PyString_AS_STRING(arg);
1004+ if (strcmp(s, "double") == 0) {
1005+ r = double_format;
1006+ }
1007+ else if (strcmp(s, "float") == 0) {
1008+ r = float_format;
1009+ }
1010+ else {
1011+ PyErr_SetString(PyExc_ValueError,
1012+ "__getformat__() argument 1 must be "
1013+ "'double' or 'float'");
1014+ return NULL;
1015+ }
1016+
1017+ switch (r) {
1018+ case unknown_format:
1019+ return PyString_FromString("unknown");
1020+ case ieee_little_endian_format:
1021+ return PyString_FromString("IEEE, little-endian");
1022+ case ieee_big_endian_format:
1023+ return PyString_FromString("IEEE, big-endian");
1024+ default:
1025+ Py_FatalError("insane float_format or double_format");
1026+ return NULL;
1027+ }
1028+}
1029+
1030+PyDoc_STRVAR(float_getformat_doc,
1031+"float.__getformat__(tstr) -> string\n"
1032+"\n"
1033+"You probably don't want to use this function. It exists mainly to be\n"
1034+"used in Python's test suite.\n"
1035+"\n"
1036+"tstr must be 'double' or 'float'. This function returns whichever of\n"
1037+"'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1038+"format of floating point numbers used by the C type named by tstr.");
1039+
1040+static PyObject *
1041+float_setformat(PyTypeObject *v, PyObject* args)
1042+{
1043+ char* tstr;
1044+ char* format;
1045+ float_format_type f;
1046+ float_format_type detected;
1047+ float_format_type *p;
1048+
1049+ if (!PyArg_ParseTuple(args, "ss:__setformat__", &tstr, &format))
1050+ return NULL;
1051+
1052+ if (strcmp(tstr, "double") == 0) {
1053+ p = &double_format;
1054+ detected = detected_double_format;
1055+ }
1056+ else if (strcmp(tstr, "float") == 0) {
1057+ p = &float_format;
1058+ detected = detected_float_format;
1059+ }
1060+ else {
1061+ PyErr_SetString(PyExc_ValueError,
1062+ "__setformat__() argument 1 must "
1063+ "be 'double' or 'float'");
1064+ return NULL;
1065+ }
1066+
1067+ if (strcmp(format, "unknown") == 0) {
1068+ f = unknown_format;
1069+ }
1070+ else if (strcmp(format, "IEEE, little-endian") == 0) {
1071+ f = ieee_little_endian_format;
1072+ }
1073+ else if (strcmp(format, "IEEE, big-endian") == 0) {
1074+ f = ieee_big_endian_format;
1075+ }
1076+ else {
1077+ PyErr_SetString(PyExc_ValueError,
1078+ "__setformat__() argument 2 must be "
1079+ "'unknown', 'IEEE, little-endian' or "
1080+ "'IEEE, big-endian'");
1081+ return NULL;
1082+
1083+ }
1084+
1085+ if (f != unknown_format && f != detected) {
1086+ PyErr_Format(PyExc_ValueError,
1087+ "can only set %s format to 'unknown' or the "
1088+ "detected platform value", tstr);
1089+ return NULL;
1090+ }
1091+
1092+ *p = f;
1093+ Py_RETURN_NONE;
1094+}
1095+
1096+PyDoc_STRVAR(float_setformat_doc,
1097+"float.__setformat__(tstr, fmt) -> None\n"
1098+"\n"
1099+"You probably don't want to use this function. It exists mainly to be\n"
1100+"used in Python's test suite.\n"
1101+"\n"
1102+"tstr must be 'double' or 'float'. fmt must be one of 'unknown',\n"
1103+"'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1104+"one of the latter two if it appears to match the underlying C reality.\n"
1105+"\n"
1106+"Overrides the automatic determination of C-level floating point type.\n"
1107+"This affects how floats are converted to and from binary strings.");
1108+
1109+static PyMethodDef float_methods[] = {
1110+ {"__getnewargs__", (PyCFunction)float_getnewargs, METH_NOARGS},
1111+ {"__getformat__", (PyCFunction)float_getformat,
1112+ METH_O|METH_CLASS, float_getformat_doc},
1113+ {"__setformat__", (PyCFunction)float_setformat,
1114+ METH_VARARGS|METH_CLASS, float_setformat_doc},
1115+ {NULL, NULL} /* sentinel */
1116+};
1117+
1118+PyDoc_STRVAR(float_doc,
1119+"float(x) -> floating point number\n\
1120+\n\
1121+Convert a string or number to a floating point number, if possible.");
1122+
1123+
1124+static PyNumberMethods float_as_number = {
1125+ float_add, /*nb_add*/
1126+ float_sub, /*nb_subtract*/
1127+ float_mul, /*nb_multiply*/
1128+ float_classic_div, /*nb_divide*/
1129+ float_rem, /*nb_remainder*/
1130+ float_divmod, /*nb_divmod*/
1131+ float_pow, /*nb_power*/
1132+ (unaryfunc)float_neg, /*nb_negative*/
1133+ (unaryfunc)float_pos, /*nb_positive*/
1134+ (unaryfunc)float_abs, /*nb_absolute*/
1135+ (inquiry)float_nonzero, /*nb_nonzero*/
1136+ 0, /*nb_invert*/
1137+ 0, /*nb_lshift*/
1138+ 0, /*nb_rshift*/
1139+ 0, /*nb_and*/
1140+ 0, /*nb_xor*/
1141+ 0, /*nb_or*/
1142+ float_coerce, /*nb_coerce*/
1143+ float_int, /*nb_int*/
1144+ float_long, /*nb_long*/
1145+ float_float, /*nb_float*/
1146+ 0, /* nb_oct */
1147+ 0, /* nb_hex */
1148+ 0, /* nb_inplace_add */
1149+ 0, /* nb_inplace_subtract */
1150+ 0, /* nb_inplace_multiply */
1151+ 0, /* nb_inplace_divide */
1152+ 0, /* nb_inplace_remainder */
1153+ 0, /* nb_inplace_power */
1154+ 0, /* nb_inplace_lshift */
1155+ 0, /* nb_inplace_rshift */
1156+ 0, /* nb_inplace_and */
1157+ 0, /* nb_inplace_xor */
1158+ 0, /* nb_inplace_or */
1159+ float_floor_div, /* nb_floor_divide */
1160+ float_div, /* nb_true_divide */
1161+ 0, /* nb_inplace_floor_divide */
1162+ 0, /* nb_inplace_true_divide */
1163+};
1164+
1165+PyTypeObject PyFloat_Type = {
1166+ PyObject_HEAD_INIT(&PyType_Type)
1167+ 0,
1168+ "float",
1169+ sizeof(PyFloatObject),
1170+ 0,
1171+ (destructor)float_dealloc, /* tp_dealloc */
1172+ (printfunc)float_print, /* tp_print */
1173+ 0, /* tp_getattr */
1174+ 0, /* tp_setattr */
1175+ 0, /* tp_compare */
1176+ (reprfunc)float_repr, /* tp_repr */
1177+ &float_as_number, /* tp_as_number */
1178+ 0, /* tp_as_sequence */
1179+ 0, /* tp_as_mapping */
1180+ (hashfunc)float_hash, /* tp_hash */
1181+ 0, /* tp_call */
1182+ (reprfunc)float_str, /* tp_str */
1183+ PyObject_GenericGetAttr, /* tp_getattro */
1184+ 0, /* tp_setattro */
1185+ 0, /* tp_as_buffer */
1186+ Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
1187+ Py_TPFLAGS_BASETYPE, /* tp_flags */
1188+ float_doc, /* tp_doc */
1189+ 0, /* tp_traverse */
1190+ 0, /* tp_clear */
1191+ float_richcompare, /* tp_richcompare */
1192+ 0, /* tp_weaklistoffset */
1193+ 0, /* tp_iter */
1194+ 0, /* tp_iternext */
1195+ float_methods, /* tp_methods */
1196+ 0, /* tp_members */
1197+ 0, /* tp_getset */
1198+ 0, /* tp_base */
1199+ 0, /* tp_dict */
1200+ 0, /* tp_descr_get */
1201+ 0, /* tp_descr_set */
1202+ 0, /* tp_dictoffset */
1203+ 0, /* tp_init */
1204+ 0, /* tp_alloc */
1205+ float_new, /* tp_new */
1206+};
1207+
1208+void
1209+_PyFloat_Init(void)
1210+{
1211+ /* We attempt to determine if this machine is using IEEE
1212+ floating point formats by peering at the bits of some
1213+ carefully chosen values. If it looks like we are on an
1214+ IEEE platform, the float packing/unpacking routines can
1215+ just copy bits, if not they resort to arithmetic & shifts
1216+ and masks. The shifts & masks approach works on all finite
1217+ values, but what happens to infinities, NaNs and signed
1218+ zeroes on packing is an accident, and attempting to unpack
1219+ a NaN or an infinity will raise an exception.
1220+
1221+ Note that if we're on some whacked-out platform which uses
1222+ IEEE formats but isn't strictly little-endian or big-
1223+ endian, we will fall back to the portable shifts & masks
1224+ method. */
1225+
1226+#if SIZEOF_DOUBLE == 8
1227+ {
1228+ double x = 9006104071832581.0;
1229+ if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
1230+ detected_double_format = ieee_big_endian_format;
1231+ else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
1232+ detected_double_format = ieee_little_endian_format;
1233+ else
1234+ detected_double_format = unknown_format;
1235+ }
1236+#else
1237+ detected_double_format = unknown_format;
1238+#endif
1239+
1240+#if SIZEOF_FLOAT == 4
1241+ {
1242+ float y = 16711938.0;
1243+ if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
1244+ detected_float_format = ieee_big_endian_format;
1245+ else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
1246+ detected_float_format = ieee_little_endian_format;
1247+ else
1248+ detected_float_format = unknown_format;
1249+ }
1250+#else
1251+ detected_float_format = unknown_format;
1252+#endif
1253+
1254+ double_format = detected_double_format;
1255+ float_format = detected_float_format;
1256+}
1257+
1258+void
1259+PyFloat_Fini(void)
1260+{
1261+ PyFloatObject *p;
1262+ PyFloatBlock *list, *next;
1263+ unsigned i;
1264+ int bc, bf; /* block count, number of freed blocks */
1265+ int frem, fsum; /* remaining unfreed floats per block, total */
1266+
1267+ bc = 0;
1268+ bf = 0;
1269+ fsum = 0;
1270+ list = block_list;
1271+ block_list = NULL;
1272+ free_list = NULL;
1273+ while (list != NULL) {
1274+ bc++;
1275+ frem = 0;
1276+ for (i = 0, p = &list->objects[0];
1277+ i < N_FLOATOBJECTS;
1278+ i++, p++) {
1279+ if (PyFloat_CheckExact(p) && p->ob_refcnt != 0)
1280+ frem++;
1281+ }
1282+ next = list->next;
1283+ if (frem) {
1284+ list->next = block_list;
1285+ block_list = list;
1286+ for (i = 0, p = &list->objects[0];
1287+ i < N_FLOATOBJECTS;
1288+ i++, p++) {
1289+ if (!PyFloat_CheckExact(p) ||
1290+ p->ob_refcnt == 0) {
1291+ p->ob_type = (struct _typeobject *)
1292+ free_list;
1293+ free_list = p;
1294+ }
1295+ }
1296+ }
1297+ else {
1298+ PyMem_FREE(list); /* XXX PyObject_FREE ??? */
1299+ bf++;
1300+ }
1301+ fsum += frem;
1302+ list = next;
1303+ }
1304+ if (!Py_VerboseFlag)
1305+ return;
1306+ fprintf(stderr, "# cleanup floats");
1307+ if (!fsum) {
1308+ fprintf(stderr, "\n");
1309+ }
1310+ else {
1311+ fprintf(stderr,
1312+ ": %d unfreed float%s in %d out of %d block%s\n",
1313+ fsum, fsum == 1 ? "" : "s",
1314+ bc - bf, bc, bc == 1 ? "" : "s");
1315+ }
1316+ if (Py_VerboseFlag > 1) {
1317+ list = block_list;
1318+ while (list != NULL) {
1319+ for (i = 0, p = &list->objects[0];
1320+ i < N_FLOATOBJECTS;
1321+ i++, p++) {
1322+ if (PyFloat_CheckExact(p) &&
1323+ p->ob_refcnt != 0) {
1324+ char buf[100];
1325+ PyFloat_AsString(buf, p);
1326+ /* XXX(twouters) cast refcount to
1327+ long until %zd is universally
1328+ available
1329+ */
1330+ fprintf(stderr,
1331+ "# <float at %p, refcnt=%ld, val=%s>\n",
1332+ p, (long)p->ob_refcnt, buf);
1333+ }
1334+ }
1335+ list = list->next;
1336+ }
1337+ }
1338+}
1339+
1340+/*----------------------------------------------------------------------------
1341+ * _PyFloat_{Pack,Unpack}{4,8}. See floatobject.h.
1342+ *
1343+ * TODO: On platforms that use the standard IEEE-754 single and double
1344+ * formats natively, these routines could simply copy the bytes.
1345+ */
1346+int
1347+_PyFloat_Pack4(double x, unsigned char *p, int le)
1348+{
1349+ if (float_format == unknown_format) {
1350+ unsigned char sign;
1351+ int e;
1352+ double f;
1353+ unsigned int fbits;
1354+ int incr = 1;
1355+
1356+ if (le) {
1357+ p += 3;
1358+ incr = -1;
1359+ }
1360+
1361+ if (x < 0) {
1362+ sign = 1;
1363+ x = -x;
1364+ }
1365+ else
1366+ sign = 0;
1367+
1368+ f = frexp(x, &e);
1369+
1370+ /* Normalize f to be in the range [1.0, 2.0) */
1371+ if (0.5 <= f && f < 1.0) {
1372+ f *= 2.0;
1373+ e--;
1374+ }
1375+ else if (f == 0.0)
1376+ e = 0;
1377+ else {
1378+ PyErr_SetString(PyExc_SystemError,
1379+ "frexp() result out of range");
1380+ return -1;
1381+ }
1382+
1383+ if (e >= 128)
1384+ goto Overflow;
1385+ else if (e < -126) {
1386+ /* Gradual underflow */
1387+ f = ldexp(f, 126 + e);
1388+ e = 0;
1389+ }
1390+ else if (!(e == 0 && f == 0.0)) {
1391+ e += 127;
1392+ f -= 1.0; /* Get rid of leading 1 */
1393+ }
1394+
1395+ f *= 8388608.0; /* 2**23 */
1396+ fbits = (unsigned int)(f + 0.5); /* Round */
1397+ assert(fbits <= 8388608);
1398+ if (fbits >> 23) {
1399+ /* The carry propagated out of a string of 23 1 bits. */
1400+ fbits = 0;
1401+ ++e;
1402+ if (e >= 255)
1403+ goto Overflow;
1404+ }
1405+
1406+ /* First byte */
1407+ *p = (sign << 7) | (e >> 1);
1408+ p += incr;
1409+
1410+ /* Second byte */
1411+ *p = (char) (((e & 1) << 7) | (fbits >> 16));
1412+ p += incr;
1413+
1414+ /* Third byte */
1415+ *p = (fbits >> 8) & 0xFF;
1416+ p += incr;
1417+
1418+ /* Fourth byte */
1419+ *p = fbits & 0xFF;
1420+
1421+ /* Done */
1422+ return 0;
1423+
1424+ Overflow:
1425+ PyErr_SetString(PyExc_OverflowError,
1426+ "float too large to pack with f format");
1427+ return -1;
1428+ }
1429+ else {
1430+ float y = (float)x;
1431+ const char *s = (char*)&y;
1432+ int i, incr = 1;
1433+
1434+ if ((float_format == ieee_little_endian_format && !le)
1435+ || (float_format == ieee_big_endian_format && le)) {
1436+ p += 3;
1437+ incr = -1;
1438+ }
1439+
1440+ for (i = 0; i < 4; i++) {
1441+ *p = *s++;
1442+ p += incr;
1443+ }
1444+ return 0;
1445+ }
1446+}
1447+
1448+int
1449+_PyFloat_Pack8(double x, unsigned char *p, int le)
1450+{
1451+ if (double_format == unknown_format) {
1452+ unsigned char sign;
1453+ int e;
1454+ double f;
1455+ unsigned int fhi, flo;
1456+ int incr = 1;
1457+
1458+ if (le) {
1459+ p += 7;
1460+ incr = -1;
1461+ }
1462+
1463+ if (x < 0) {
1464+ sign = 1;
1465+ x = -x;
1466+ }
1467+ else
1468+ sign = 0;
1469+
1470+ f = frexp(x, &e);
1471+
1472+ /* Normalize f to be in the range [1.0, 2.0) */
1473+ if (0.5 <= f && f < 1.0) {
1474+ f *= 2.0;
1475+ e--;
1476+ }
1477+ else if (f == 0.0)
1478+ e = 0;
1479+ else {
1480+ PyErr_SetString(PyExc_SystemError,
1481+ "frexp() result out of range");
1482+ return -1;
1483+ }
1484+
1485+ if (e >= 1024)
1486+ goto Overflow;
1487+ else if (e < -1022) {
1488+ /* Gradual underflow */
1489+ f = ldexp(f, 1022 + e);
1490+ e = 0;
1491+ }
1492+ else if (!(e == 0 && f == 0.0)) {
1493+ e += 1023;
1494+ f -= 1.0; /* Get rid of leading 1 */
1495+ }
1496+
1497+ /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
1498+ f *= 268435456.0; /* 2**28 */
1499+ fhi = (unsigned int)f; /* Truncate */
1500+ assert(fhi < 268435456);
1501+
1502+ f -= (double)fhi;
1503+ f *= 16777216.0; /* 2**24 */
1504+ flo = (unsigned int)(f + 0.5); /* Round */
1505+ assert(flo <= 16777216);
1506+ if (flo >> 24) {
1507+ /* The carry propagated out of a string of 24 1 bits. */
1508+ flo = 0;
1509+ ++fhi;
1510+ if (fhi >> 28) {
1511+ /* And it also progagated out of the next 28 bits. */
1512+ fhi = 0;
1513+ ++e;
1514+ if (e >= 2047)
1515+ goto Overflow;
1516+ }
1517+ }
1518+
1519+ /* First byte */
1520+ *p = (sign << 7) | (e >> 4);
1521+ p += incr;
1522+
1523+ /* Second byte */
1524+ *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
1525+ p += incr;
1526+
1527+ /* Third byte */
1528+ *p = (fhi >> 16) & 0xFF;
1529+ p += incr;
1530+
1531+ /* Fourth byte */
1532+ *p = (fhi >> 8) & 0xFF;
1533+ p += incr;
1534+
1535+ /* Fifth byte */
1536+ *p = fhi & 0xFF;
1537+ p += incr;
1538+
1539+ /* Sixth byte */
1540+ *p = (flo >> 16) & 0xFF;
1541+ p += incr;
1542+
1543+ /* Seventh byte */
1544+ *p = (flo >> 8) & 0xFF;
1545+ p += incr;
1546+
1547+ /* Eighth byte */
1548+ *p = flo & 0xFF;
1549+ p += incr;
1550+
1551+ /* Done */
1552+ return 0;
1553+
1554+ Overflow:
1555+ PyErr_SetString(PyExc_OverflowError,
1556+ "float too large to pack with d format");
1557+ return -1;
1558+ }
1559+ else {
1560+ const char *s = (char*)&x;
1561+ int i, incr = 1;
1562+
1563+ if ((double_format == ieee_little_endian_format && !le)
1564+ || (double_format == ieee_big_endian_format && le)) {
1565+ p += 7;
1566+ incr = -1;
1567+ }
1568+
1569+ for (i = 0; i < 8; i++) {
1570+ *p = *s++;
1571+ p += incr;
1572+ }
1573+ return 0;
1574+ }
1575+}
1576+
1577+double
1578+_PyFloat_Unpack4(const unsigned char *p, int le)
1579+{
1580+ if (float_format == unknown_format) {
1581+ unsigned char sign;
1582+ int e;
1583+ unsigned int f;
1584+ double x;
1585+ int incr = 1;
1586+
1587+ if (le) {
1588+ p += 3;
1589+ incr = -1;
1590+ }
1591+
1592+ /* First byte */
1593+ sign = (*p >> 7) & 1;
1594+ e = (*p & 0x7F) << 1;
1595+ p += incr;
1596+
1597+ /* Second byte */
1598+ e |= (*p >> 7) & 1;
1599+ f = (*p & 0x7F) << 16;
1600+ p += incr;
1601+
1602+ if (e == 255) {
1603+ PyErr_SetString(
1604+ PyExc_ValueError,
1605+ "can't unpack IEEE 754 special value "
1606+ "on non-IEEE platform");
1607+ return -1;
1608+ }
1609+
1610+ /* Third byte */
1611+ f |= *p << 8;
1612+ p += incr;
1613+
1614+ /* Fourth byte */
1615+ f |= *p;
1616+
1617+ x = (double)f / 8388608.0;
1618+
1619+ /* XXX This sadly ignores Inf/NaN issues */
1620+ if (e == 0)
1621+ e = -126;
1622+ else {
1623+ x += 1.0;
1624+ e -= 127;
1625+ }
1626+ x = ldexp(x, e);
1627+
1628+ if (sign)
1629+ x = -x;
1630+
1631+ return x;
1632+ }
1633+ else {
1634+ float x;
1635+
1636+ if ((float_format == ieee_little_endian_format && !le)
1637+ || (float_format == ieee_big_endian_format && le)) {
1638+ char buf[4];
1639+ char *d = &buf[3];
1640+ int i;
1641+
1642+ for (i = 0; i < 4; i++) {
1643+ *d-- = *p++;
1644+ }
1645+ memcpy(&x, buf, 4);
1646+ }
1647+ else {
1648+ memcpy(&x, p, 4);
1649+ }
1650+
1651+ return x;
1652+ }
1653+}
1654+
1655+double
1656+_PyFloat_Unpack8(const unsigned char *p, int le)
1657+{
1658+ if (double_format == unknown_format) {
1659+ unsigned char sign;
1660+ int e;
1661+ unsigned int fhi, flo;
1662+ double x;
1663+ int incr = 1;
1664+
1665+ if (le) {
1666+ p += 7;
1667+ incr = -1;
1668+ }
1669+
1670+ /* First byte */
1671+ sign = (*p >> 7) & 1;
1672+ e = (*p & 0x7F) << 4;
1673+
1674+ p += incr;
1675+
1676+ /* Second byte */
1677+ e |= (*p >> 4) & 0xF;
1678+ fhi = (*p & 0xF) << 24;
1679+ p += incr;
1680+
1681+ if (e == 2047) {
1682+ PyErr_SetString(
1683+ PyExc_ValueError,
1684+ "can't unpack IEEE 754 special value "
1685+ "on non-IEEE platform");
1686+ return -1.0;
1687+ }
1688+
1689+ /* Third byte */
1690+ fhi |= *p << 16;
1691+ p += incr;
1692+
1693+ /* Fourth byte */
1694+ fhi |= *p << 8;
1695+ p += incr;
1696+
1697+ /* Fifth byte */
1698+ fhi |= *p;
1699+ p += incr;
1700+
1701+ /* Sixth byte */
1702+ flo = *p << 16;
1703+ p += incr;
1704+
1705+ /* Seventh byte */
1706+ flo |= *p << 8;
1707+ p += incr;
1708+
1709+ /* Eighth byte */
1710+ flo |= *p;
1711+
1712+ x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
1713+ x /= 268435456.0; /* 2**28 */
1714+
1715+ if (e == 0)
1716+ e = -1022;
1717+ else {
1718+ x += 1.0;
1719+ e -= 1023;
1720+ }
1721+ x = ldexp(x, e);
1722+
1723+ if (sign)
1724+ x = -x;
1725+
1726+ return x;
1727+ }
1728+ else {
1729+ double x;
1730+
1731+ if ((double_format == ieee_little_endian_format && !le)
1732+ || (double_format == ieee_big_endian_format && le)) {
1733+ char buf[8];
1734+ char *d = &buf[7];
1735+ int i;
1736+
1737+ for (i = 0; i < 8; i++) {
1738+ *d-- = *p++;
1739+ }
1740+ memcpy(&x, buf, 8);
1741+ }
1742+ else {
1743+ memcpy(&x, p, 8);
1744+ }
1745+
1746+ return x;
1747+ }
1748+}