/* Module implementing high-precision natural number arithmetic in a decimal base. As the name 'deccoeff' suggests, these numbers are intended to be used as the coefficients for Decimal instances. Author: Mark Dickinson. Licensed to PSF under a Contributor Agreement. */ #include "Python.h" #include "longintrepr.h" #include "string.h" /* deccoeffs are represented in base BASE, for BASE a suitable power of 10. I'll use the word 'limb' to refer to a base BASE digit, and 'digit' to refer to a usual decimal digit in the range 0 through 9. Thus a single limb holds some constant number of digits; that constant is called LIMB_DIGITS below. The C type 'limb' should be large enough to hold a single digit in this base: i.e. all numbers in the range [0, BASE). Type 'double_limb' should be large enough to hold all numbers in the range [0, BASE*BASE). */ #define MODULE_NAME "deccoeff" #define CLASS_NAME "Deccoeff" #ifdef HAVE_LONG_LONG #define LIMB_DIGITS 9 /* BASE == 10**LIMB_DIGITS */ #define LIMB_MAX 999999999 /* BASE - 1 */ typedef unsigned long limb; typedef unsigned PY_LONG_LONG double_limb; static limb powers_of_ten[LIMB_DIGITS+1] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000 }; #else #define LIMB_DIGITS 4 #define LIMB_MAX 9999 typedef unsigned short limb; typedef unsigned long double_limb; static limb powers_of_ten[LIMB_DIGITS+1] = { 1, 10, 100, 1000, 10000 }; #endif /* macros to extract the high and low limbs from a double_limb, and to combine high and low limbs into a double_limb. */ #define deccoeff_BASE (LIMB_MAX+1) #define HIGH_LIMB(x) ((limb)(x/deccoeff_BASE)) #define LOW_LIMB(x) ((limb)(x%deccoeff_BASE)) #define HILO(x, y) ((double_limb)x * deccoeff_BASE + y) /*********************** * Arithmetic on limbs * ***********************/ /* add with carry: compute a + b + carry; put result in *result and return new carry. Here, carry should be 0 or 1. */ static int limb_addc(limb *result, limb a, limb b, int carry) { assert(a <= LIMB_MAX && b <= LIMB_MAX); assert(carry == 0 || carry == 1); /* XXX the code below takes care never to do any computation that might give a result outside the range [0, LIMB_MAX]. In fact, unless we're doing something crazy like having a limb be 19 digits then 2*LIMB_MAX+1 still fits in a limb; we could use this to make the code below simpler. */ if (carry) { if (a >= LIMB_MAX - b) { *result = a - (LIMB_MAX-b); return 1; } else { *result = a + b + 1; return 0; } } else { if (a > LIMB_MAX - b) { *result = a - 1 - (LIMB_MAX-b); return 1; } else { *result = a + b; return 0; } } } /* fused multiply-add: compute a * b + c; put less significant limb of result in *low, and return more significant limb. */ static limb limb_fma(limb *low, limb a, limb b, limb c) { double_limb hilo; assert(a <= LIMB_MAX && b <= LIMB_MAX && c <= LIMB_MAX); hilo = (double_limb)a * b + c; *low = LOW_LIMB(hilo); return HIGH_LIMB(hilo); } /* fused multiply-add with two addends: compute a*b + c + d; return results as in limb_fma. This is needed in the multiplication algorithm. */ static limb limb_fmaa(limb *low, limb a, limb b, limb c, limb d) { double_limb hilo; assert(a <= LIMB_MAX && b <= LIMB_MAX && c <= LIMB_MAX && d <= LIMB_MAX); hilo = (double_limb)a * b + c + d; *low = LOW_LIMB(hilo); return HIGH_LIMB(hilo); } /* division, quotient only: divide high*BASE + low by c, and return the integer part of the result. Assumes that high < c, so that the quotient fits into a single limb. */ static limb limb_quot(limb high, limb low, limb c) { double_limb hilo; assert(high < c && low <= LIMB_MAX && c <= LIMB_MAX); hilo = HILO(high, low); return (limb)(hilo/c); } /* division with quotient and remainder. Like limb_quot, except that it also returns the remainder in *rem. */ static limb limb_quotrem(limb *rem, limb high, limb low, limb c) { double_limb hilo; assert(high <= LIMB_MAX && low <= LIMB_MAX && c <= LIMB_MAX); assert(high < c); hilo = HILO(high, low); *rem = (limb)(hilo%c); return (limb)(hilo/c); } /* turn a limb, all of whose digits are 0 or 1, into a binary equivalent. Returns 1 on success, 0 if any of the digits was not 0 or 1. */ static int limb_logical_to_binary(limb *result, limb a) { limb acc; int i; /* shift digits out of a and shift corresponding bits into acc. */ acc = 0; for (i = 0; i < LIMB_DIGITS; i++) { if (a%10 >= 2) return 0; acc = (acc << 1) + (a%10); a /= 10; } *result = acc; return 1; } /* ...and convert back again... */ static limb limb_binary_to_logical(limb a) { limb acc; int i; acc = 0; for (i=0; i < LIMB_DIGITS; i++) { acc = acc*10 + (a&1); a >>= 1; } return acc; } /************************** * deccoeff : definitions * **************************/ typedef struct { PyObject_VAR_HEAD limb ob_limbs[0]; } deccoeff; /* the number of decimal digits (*not* the number of limbs) should fit in a Py_ssize_t, so that length and indexing always make sense; this means that the number of limbs is limited by MAX_LIMBS = ceiling(MAX_DIGITS/LIMB_DIGITS). MAX_LIMBS should be at most PY_SSIZE_T_MAX-1; this makes it possible to allocate (MAX_LIMBS+1) limbs in situations where eventual normalization might result in a value that has only MAX_LIMBS limbs. Without this, the multiplication algorithm would be more complicated. MAX_DIGITS should also be at least 1. */ #if LIMB_DIGITS == 1 #define MAX_DIGITS (PY_SSIZE_T_MAX-1) #else #define MAX_DIGITS PY_SSIZE_T_MAX #endif #define TOP_LIMB_BOUND (powers_of_ten[(MAX_DIGITS-1)%LIMB_DIGITS+1]-1) #define MAX_LIMBS ((MAX_DIGITS-1)/LIMB_DIGITS + 1) #define GETLIMB(v, i) ((i < Py_SIZE(v) && i >= 0) ? v->ob_limbs[i] : 0) static deccoeff *deccoeff_zero; static deccoeff *deccoeff_one; static deccoeff *deccoeff_two; static deccoeff *deccoeff_PyLong_BASE; static PyTypeObject deccoeff_DeccoeffType; /* allocate a new decimal integer with 'size' limbs */ static deccoeff * _deccoeff_new(Py_ssize_t size) { /* should never be allocating more than MAX_LIMBS+1 limbs */ assert(size <= (MAX_LIMBS+1)); return PyObject_NEW_VAR(deccoeff, &deccoeff_DeccoeffType, size); } /* remove leading zero limbs from a possibly unnormalized deccoeff. Actually, this just modifies the recorded size; it doesn't allocate a new object. Returns the same object, without affecting reference counts. */ static deccoeff * deccoeff_normalize(deccoeff *v) { Py_ssize_t v_size, i; i = v_size = Py_SIZE(v); while (i > 0 && v->ob_limbs[i-1] == 0) i--; Py_SIZE(v) = i; return v; } /* return true if v has MAX_DIGITS or fewer digits. */ static int deccoeff_checksize(deccoeff *v) { Py_ssize_t v_size; v_size = Py_SIZE(v); return (v_size < MAX_LIMBS || (v_size == MAX_LIMBS && v->ob_limbs[v_size-1] <= TOP_LIMB_BOUND)); } static deccoeff * deccoeff_copy(deccoeff *v) { deccoeff *w; Py_ssize_t v_size, i; v_size = Py_SIZE(v); w = _deccoeff_new(v_size); if (w == NULL) return w; for (i = 0; i < v_size; i++) w->ob_limbs[i] = v->ob_limbs[i]; return w; } /* return a new reference to the zero deccoeff */ static deccoeff * zero(void) { Py_INCREF(deccoeff_zero); return deccoeff_zero; } /* Create a deccoeff from a C unsigned long. */ #if LIMB_MAX < ULONG_MAX #define deccoeff_BASE (LIMB_MAX+1) static deccoeff * deccoeff_from_ulong(unsigned long x) { unsigned long xcopy; Py_ssize_t z_size; deccoeff *z; limb *z_limb, *z_end; /* wasteful method of figuring out number of limbs needed */ z_size = 0; xcopy = x; while(xcopy > 0) { xcopy /= deccoeff_BASE; z_size++; } z = _deccoeff_new(z_size); if (z==NULL) return NULL; z_limb = z->ob_limbs; z_end = z_limb+z_size; while (z_limb < z_end) { *z_limb++ = x % deccoeff_BASE; x /= deccoeff_BASE; } return z; } #undef deccoeff_BASE #else static deccoeff * deccoeff_from_ulong(unsigned long x) { if (x == 0) return zero(); z = _deccoeff_new(1); if (z==NULL) return NULL; z->ob_limbs[0] = (limb)x; return z; } #endif /*************************** * Arithmetic on deccoeffs * ***************************/ /* logical and. raise ValueError if either operand is not a logical operand */ static deccoeff * deccoeff_and(deccoeff *a, deccoeff *b) { Py_ssize_t a_size, b_size, z_size, size_temp; deccoeff *temp, *z; limb abin, bbin; int i; a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size < b_size) { temp = a; a = b; b = temp; size_temp = a_size; a_size = b_size; b_size = size_temp; } z_size = b_size; z = _deccoeff_new(z_size); if (z == NULL) return NULL; for (i = 0; i < b_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i]) || !limb_logical_to_binary(&bbin, b->ob_limbs[i])) goto Fail; z->ob_limbs[i] = limb_binary_to_logical(abin & bbin); } for (; i < a_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i])) goto Fail; } return deccoeff_normalize(z); Fail: Py_DECREF(z); PyErr_SetString(PyExc_ValueError, "not a logical operand in logical operation"); return NULL; } /* logical xor. raise ValueError if either operand is not a logical operand */ static deccoeff * deccoeff_xor(deccoeff *a, deccoeff *b) { Py_ssize_t a_size, b_size, z_size, size_temp; deccoeff *temp, *z; limb abin, bbin; int i; a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size < b_size) { temp = a; a = b; b = temp; size_temp = a_size; a_size = b_size; b_size = size_temp; } z_size = a_size; z = _deccoeff_new(z_size); if (z == NULL) return NULL; for (i = 0; i < b_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i]) || !limb_logical_to_binary(&bbin, b->ob_limbs[i])) goto Fail; z->ob_limbs[i] = limb_binary_to_logical(abin ^ bbin); } for (; i < a_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i])) goto Fail; z->ob_limbs[i] = a->ob_limbs[i]; } return deccoeff_normalize(z); Fail: Py_DECREF(z); PyErr_SetString(PyExc_ValueError, "not a logical operand in logical operation"); return NULL; } /* logical or. raise ValueError if either operand is not a logical operand */ static deccoeff * deccoeff_or(deccoeff *a, deccoeff *b) { Py_ssize_t a_size, b_size, z_size, size_temp; deccoeff *temp, *z; limb abin, bbin; int i; a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size < b_size) { temp = a; a = b; b = temp; size_temp = a_size; a_size = b_size; b_size = size_temp; } z_size = a_size; z = _deccoeff_new(z_size); if (z == NULL) return NULL; for (i = 0; i < b_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i]) || !limb_logical_to_binary(&bbin, b->ob_limbs[i])) goto Fail; z->ob_limbs[i] = limb_binary_to_logical(abin | bbin); } for (; i < a_size; i++) { if (!limb_logical_to_binary(&abin, a->ob_limbs[i])) goto Fail; z->ob_limbs[i] = a->ob_limbs[i]; } /* no need to normalize here */ return z; Fail: Py_DECREF(z); PyErr_SetString(PyExc_ValueError, "not a logical operand in logical operation"); return NULL; } /* addition: raise OverflowError if result has more than MAX_DIGITS digits */ static deccoeff * deccoeff_add(deccoeff *a, deccoeff *b) { Py_ssize_t z_size, a_size, b_size, size_temp, i; limb al, bl; deccoeff *temp, *z; int extra_limb, carry; /* exchange so that a is at least as long as b */ a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size < b_size) { temp = a; a = b; b = temp; size_temp = a_size; a_size = b_size; b_size = size_temp; } /* the chance of needing an extra limb is fairly small, so it's wasteful to allocate a_size+1 limbs each time. The following code determines whether we need an extra limb. */ extra_limb = 0; for (i = a_size; i > 0; i--) { al = a->ob_limbs[i-1]; bl = i <= b_size ? b->ob_limbs[i-1] : 0; if (bl != LIMB_MAX - al) { if (bl > LIMB_MAX - al) extra_limb = 1; break; } } z_size = a_size + extra_limb; z = _deccoeff_new(z_size); if (z == NULL) return NULL; carry = 0; for (i = 0; i < b_size; i++) carry = limb_addc(&(z->ob_limbs[i]), a->ob_limbs[i], b->ob_limbs[i], carry); for (; i < a_size; i++) carry = limb_addc(&(z->ob_limbs[i]), a->ob_limbs[i], 0, carry); if (carry) z->ob_limbs[a_size] = (limb)carry; /* result should already be normalized */ if (deccoeff_checksize(z)) return z; /* z has more than MAX_DIGITS digits. Raise OverflowError. */ Py_DECREF(z); PyErr_SetString(PyExc_OverflowError, "Too many digits in sum"); return NULL; } /* compare: return -1 if a < b, 0 if a == b, 1 if a > b */ int deccoeff_compare(deccoeff *a, deccoeff *b) { Py_ssize_t a_size, b_size; limb *a_top, *b_top, *a_start; limb a_limb, b_limb; /* compare sizes */ a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size != b_size) return a_size < b_size ? -1 : 1; /* sizes are identical; compare limbs */ a_start = a->ob_limbs; a_top = a_start + a_size; b_top = b->ob_limbs + a_size; while (a_top > a_start) { a_limb = *--a_top; b_limb = *--b_top; if (a_limb != b_limb) return a_limb < b_limb ? -1 : 1; } /* sizes and limbs both match */ return 0; } static PyObject * deccoeff_richcompare(PyObject *self, PyObject *other, int op) { PyObject *result; result = Py_CmpToRich(op, deccoeff_compare((deccoeff*)self, (deccoeff*)other)); return result; } /* subtraction: return a - b; raise OverflowError if the result is negative */ static deccoeff * deccoeff_sub(deccoeff *a, deccoeff *b) { deccoeff *z; Py_ssize_t a_size, b_size, i; int carry; a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size < b_size) goto Overflow; z = _deccoeff_new(a_size); if (z==NULL) return NULL; carry = 1; for (i = 0; i < b_size; i++) carry = limb_addc(&(z->ob_limbs[i]), a->ob_limbs[i], LIMB_MAX-b->ob_limbs[i], carry); for (; i < a_size; i++) carry = limb_addc(&(z->ob_limbs[i]), a->ob_limbs[i], LIMB_MAX, carry); if (carry == 1) return deccoeff_normalize(z); Py_DECREF(z); Overflow: PyErr_SetString(PyExc_OverflowError, "difference is negative and cannot be represented as a Deccoeff"); return NULL; } /* multiplication: return a * b. If the product has more than MAX_DIGITS digits, raise OverflowError and return NULL. */ static deccoeff * deccoeff_mul(deccoeff *a, deccoeff *b) { Py_ssize_t a_size, b_size, z_size, i, j; limb carry, a_limb; deccoeff *z; a_size = Py_SIZE(a); b_size = Py_SIZE(b); if (a_size == 0 || b_size == 0) return zero(); /* number of limbs of result is either a_size + b_size or a_size + b_size - 1. So if a_size + b_size > MAX_LIMBS + 1 then the result is too large to be represented. Otherwise, if a_size + b_size - 1 <= MAX_LIMBS, we're allocating fewer than MAX_LIMBS+1 limbs. */ if (a_size > MAX_LIMBS + 1 - b_size) goto Overflow; /* space for result */ z_size = a_size + b_size; z = _deccoeff_new(z_size); if (z == NULL) return NULL; for (j = 0; j < b_size; j++) z->ob_limbs[j] = 0; for (i = 0; i < a_size; i++) { a_limb = a->ob_limbs[i]; carry = 0; for (j = 0; j < b_size; j++) { carry = limb_fmaa(z->ob_limbs+(i+j), a->ob_limbs[i], b->ob_limbs[j], z->ob_limbs[i+j], carry); } z-> ob_limbs[i+b_size] = carry; } z = deccoeff_normalize(z); if (deccoeff_checksize(z)) return z; Py_DECREF(z); Overflow: PyErr_SetString(PyExc_OverflowError, "product has too many digits"); return NULL; } /* division, returning quotient and remainder. Raises ZeroDivisionError on division by zero. */ /* division; returns quotient and puts remainder in *rem. On failure, returns NULL and puts NULL in *rem. */ static deccoeff * _deccoeff_divmod(deccoeff **rem, deccoeff *a, deccoeff *b) { deccoeff *qq, *aa, *bb; Py_ssize_t a_size, b_size, q_size, i, aa_size; limb q, scale; limb mul_carry, low, top_limb; limb *a_limb, *b_limb, *q_limb, *aa_limb, *bb_limb; int extra_limb, sub_carry; a_size = Py_SIZE(a); b_size = Py_SIZE(b); *rem = NULL; if (b_size == 0) { PyErr_SetString(PyExc_ZeroDivisionError, "division by zero"); return NULL; } /* if top limb of a >= top limb of b, we'll need to pad a with an extra limb */ extra_limb = a_size > 0 && a->ob_limbs[a_size-1] >= b->ob_limbs[b_size-1]; aa_size = a_size + extra_limb; /* if a is clearly less than b, then quotient is zero, remainder is a */ if (aa_size <= b_size) { Py_INCREF(a); *rem = a; return zero(); } assert(a_size > 0 && b_size > 0); /* scale factor for normalization */ top_limb = b->ob_limbs[b_size-1]; scale = 1 + (LIMB_MAX-top_limb)/(top_limb+1); /* scale b */ bb = _deccoeff_new(b_size); if (bb == NULL) return NULL; mul_carry = 0; bb_limb = bb->ob_limbs; b_limb = b->ob_limbs; for (i = 0; i < b_size; i++) mul_carry = limb_fma(bb->ob_limbs + i, b->ob_limbs[i], scale, mul_carry); assert(mul_carry == 0); b = bb; /* scale a */ aa = _deccoeff_new(aa_size); if (aa == NULL) { Py_DECREF(b); return NULL; } mul_carry = 0; aa_limb = aa->ob_limbs; a_limb = a->ob_limbs; for (i=0; i < a_size; i++) mul_carry = limb_fma(aa_limb++, *a_limb++, scale, mul_carry); if (a_size != aa_size) *aa_limb++ = mul_carry; else assert(mul_carry == 0); a = aa; a_size = aa_size; /* note that we now own the references to a, b */ /* space for quotient */ q_size = a_size - b_size; qq = _deccoeff_new(q_size); if (qq == NULL) { Py_DECREF(a); Py_DECREF(b); return NULL; } q_limb = qq->ob_limbs + q_size; while (q_limb > qq->ob_limbs) { assert(a_size > b_size); assert(a->ob_limbs[a_size-1] <= b->ob_limbs[b_size-1]); /* get estimate for quotient */ if (a->ob_limbs[a_size-1] == b->ob_limbs[b_size-1]) { q = LIMB_MAX; } else { /* divide top two limbs of a by top limb of b; this should only ever produce an overestimate for the true quotient */ assert(a_size >= 2); q = limb_quot(a->ob_limbs[a_size-1], a->ob_limbs[a_size-2], b->ob_limbs[b_size-1]); assert(q <= LIMB_MAX); } /* now we want to replace a with a - q*b */ mul_carry = 0; /* carry from multiplication */ sub_carry = 1; a_limb = a->ob_limbs + (a_size - b_size - 1); b_limb = b->ob_limbs; for (i = 0; i < b_size; i++) { mul_carry = limb_fma(&low, *b_limb++, q, mul_carry); sub_carry = limb_addc(a_limb, *a_limb, LIMB_MAX-low, sub_carry); a_limb++; } sub_carry = limb_addc(a_limb, *a_limb, LIMB_MAX-mul_carry, sub_carry); a_limb++; assert(a_limb == a->ob_limbs + a_size); /* while result is negative, re-add b */ while (sub_carry == 0) { assert(q > 0); --q; a_limb = a->ob_limbs + (a_size - b_size - 1); b_limb = b->ob_limbs; /* add back */ for (i=0; iob_limbs + a_size); } /* at this point, the top limb of a should be 0 */ assert(a->ob_limbs[a_size-1] == 0); a_size--; *--q_limb = q; } assert(a_size == b_size); /* remainder is normalization of a; need to divide through by scale */ /* short division */ mul_carry = 0; for (i = a_size; i > 0; i-- ) a->ob_limbs[i-1] = limb_quotrem(&mul_carry, mul_carry, a->ob_limbs[i-1], scale); /* division should have been exact */ assert(mul_carry == 0); Py_DECREF(b); *rem = deccoeff_normalize(a); return deccoeff_normalize(qq); } static PyObject * deccoeff_divmod(deccoeff *a, deccoeff *b) { deccoeff *quot, *rem; quot = _deccoeff_divmod(&rem, a, b); if (quot == NULL) return NULL; return Py_BuildValue("OO", (PyObject *)quot, (PyObject *)rem); } static deccoeff * deccoeff_floor_divide(deccoeff *a, deccoeff *b) { deccoeff *quot, *rem; quot = _deccoeff_divmod(&rem, a, b); if (quot == NULL) return NULL; Py_DECREF(rem); return quot; } static deccoeff * deccoeff_remainder(deccoeff *a, deccoeff *b) { deccoeff *quot, *rem; quot = _deccoeff_divmod(&rem, a, b); if (quot == NULL) return NULL; Py_DECREF(quot); return rem; } /* halve b in place, returning the lowest bit. Be careful only to use this on 'private' deccoeffs, since any deccoeff the Python interpreter sees is supposed to be immutable. */ #define HALF_BASE (LIMB_MAX/2 + 1) int halve_in_place(deccoeff *b) { Py_ssize_t i; int carry; limb limb; if (Py_SIZE(b) == 0) return 0; limb = b->ob_limbs[Py_SIZE(b)-1]; if (limb < 2) { carry = (int)limb; Py_SIZE(b)--; } else { carry = 0; } for (i=Py_SIZE(b); i>0; i--) { limb = b->ob_limbs[i-1]; b->ob_limbs[i-1] = (limb >> 1) + (carry ? HALF_BASE : 0); carry = limb & 1; } return carry; } /* power: a**b. Slow algorithm involving many unnecessary allocations and deallocations of memory. */ /* To do: speed this up, and consider special-casing 2**n, since it's needed for shifts. */ static deccoeff * deccoeff_pow(deccoeff *a, deccoeff *b, PyObject *c) { deccoeff *acc, *temp; int lowbit; if (c != Py_None) PyErr_SetString(PyExc_TypeError, "deccoeff does not support 3rd argument to pow"); if (Py_SIZE(b) == 0) { Py_INCREF(deccoeff_one); return deccoeff_one; } /* make a copy of b, which we'll modify in-place */ b = deccoeff_copy(b); if (b == NULL) return NULL; /* we now own a reference to b */ Py_INCREF(a); /* result accumulates in acc */ acc = deccoeff_one; Py_INCREF(acc); while (1) { /* invariant quantity: a**b*acc. */ /* b > 0. */ lowbit = halve_in_place(b); if (lowbit) { /* acc *= a */ temp = deccoeff_mul(a, acc); if (temp == NULL) { Py_DECREF(a); Py_DECREF(b); Py_DECREF(acc); return NULL; } Py_DECREF(acc); acc = temp; } if (Py_SIZE(b) == 0) break; /* a *= a */ temp = deccoeff_mul(a, a); Py_DECREF(a); if (temp == NULL) { Py_DECREF(b); Py_DECREF(acc); return NULL; } a = temp; } Py_DECREF(a); Py_DECREF(b); return acc; } /* compute a << b */ static deccoeff * deccoeff_lshift(deccoeff *a, deccoeff *b) { deccoeff *c, *temp; /* 2**b */ c = deccoeff_pow(deccoeff_two, b, Py_None); if (c == NULL) return c; /* c = a*c */ temp = deccoeff_mul(a, c); Py_DECREF(c); c = temp; return c; } /* compute a >> b */ static deccoeff * deccoeff_rshift(deccoeff *a, deccoeff *b) { deccoeff *c, *temp; /* 2**b */ c = deccoeff_pow(deccoeff_two, b, Py_None); if (c == NULL) return c; /* c = a // c */ temp = deccoeff_floor_divide(a, c); Py_DECREF(c); c = temp; return c; } /* Create a deccoeff from a Python long integer */ /* Naive algorithm; needs optimization. But not before testing with a range of limb sizes. */ static deccoeff * deccoeff_from_PyLong(deccoeff *self, PyObject *o) { Py_ssize_t a_size; PyLongObject *a; digit *a_start, *a_end; deccoeff *result, *result2, *addend; if (!PyLong_Check(o)) { PyErr_SetString(PyExc_TypeError, "argument must be a long"); return NULL; } a = (PyLongObject *)o; a_size = Py_SIZE(a); if (a_size < 0) { PyErr_SetString(PyExc_OverflowError, "Can't convert negative integer to a deccoeff"); return NULL; } result = deccoeff_zero; Py_INCREF(result); a_start = a->ob_digit; a_end = a_start+a_size; while (a_end > a_start) { /* multiply result by PyLong_BASE, freeing old reference */ result2 = deccoeff_mul(deccoeff_PyLong_BASE, result); Py_DECREF(result); result = result2; if (result == NULL) return NULL; /* add PyLong digit to result */ addend = deccoeff_from_ulong((unsigned long)(*--a_end)); if (addend == NULL) { Py_DECREF(result); return NULL; } result2 = deccoeff_add(result, addend); Py_DECREF(result); Py_DECREF(addend); result = result2; if (result == NULL) return NULL; } return result; } /* we also need (eventually) a routine to convert a deccoeff back into a Python long. I'm resisting writing this for the moment, since the lack of it provides incentive to keep all computations in decimal. */ static PyObject * deccoeff_new(PyTypeObject *type, PyObject *args, PyObject *kwds) { const char* s; char *start, *end; int lens; int nlimbs, digits_in_limb; limb accum, *limb_pointer; char c; deccoeff *new; if (!PyArg_ParseTuple(args, "s#:" CLASS_NAME, &s, &lens)) return NULL; if (lens > MAX_DIGITS) { PyErr_SetString(PyExc_OverflowError, "too many digits"); return NULL; } nlimbs = (lens+LIMB_DIGITS-1)/LIMB_DIGITS; digits_in_limb = (lens+LIMB_DIGITS-1) % LIMB_DIGITS + 1; new = _deccoeff_new(nlimbs); if (new == NULL) return NULL; /* iterate through digits, from most significant to least */ start = (char *)s; end = (char *)s + lens; accum = 0; limb_pointer = new -> ob_limbs + nlimbs - 1; while (start < end) { c = *start++; if (!(c >= '0') || !(c <= '9')) goto Error; accum = accum*10 + (c-'0'); digits_in_limb--; if (digits_in_limb == 0) { digits_in_limb = LIMB_DIGITS; *limb_pointer-- = accum; accum = 0; } } return (PyObject *)(deccoeff_normalize(new)); Error: PyErr_Format(PyExc_ValueError, "invalid literal for deccoeff(): %s", s); Py_DECREF(new); return NULL; } /* number of digits of a deccoeff. Raises ValueError if the deccoeff is zero. */ static Py_ssize_t deccoeff_length(deccoeff *v) { Py_ssize_t sz, v_size; limb top_limb; v_size = Py_SIZE(v); /* for safety, we raise a ValueError for 0; there are too many situations where returning 0 might turn out to be wrong; 0 almost always needs special treatment anyway. Another way to think of this function is as 1 more than the floor of log_10(v); this again is undefined for v = 0, justifying the ValueError. */ if (v_size == 0) { PyErr_SetString(PyExc_ValueError, "refusing to find length of 0"); return -1; } sz = 1; top_limb = v->ob_limbs[v_size-1]; while (top_limb >= powers_of_ten[sz]) sz++; /* since number of digits is limited to MAX_DIGITS, there should be no danger of overflow here. */ return sz + (v_size-1) * LIMB_DIGITS; } static int deccoeff_nonzero(deccoeff *v) { return Py_SIZE(v) != 0; } static PyObject * deccoeff_str(deccoeff *v) { Py_ssize_t sz, nlimbs; limb *limb_pointer, *last_limb, limb_value; PyObject *str; int i; Py_UNICODE *p; nlimbs = Py_SIZE(v); if (nlimbs == 0) { /* return empty string */ str = PyUnicode_FromUnicode(NULL, 0); if (str == NULL) return NULL; p = PyUnicode_AS_UNICODE(str); *p = '\0'; return str; } sz = deccoeff_length(v); str = PyUnicode_FromUnicode(NULL, sz); if (str == NULL) return NULL; p = PyUnicode_AS_UNICODE(str) + sz; *p = '\0'; /* fill in digits from right to left; start with the least significant limb */ limb_pointer = v -> ob_limbs; last_limb = limb_pointer + nlimbs - 1; while (limb_pointer < last_limb) { limb_value = *limb_pointer++; for (i=0; i < LIMB_DIGITS; i++) { *--p = '0' + limb_value % 10; limb_value /= 10; } } /* most significant limb */ limb_value = *limb_pointer; assert(limb_value != 0); while (limb_value != 0) { *--p = '0' + limb_value % 10; limb_value /= 10; } return str; } static PyObject * deccoeff_repr(deccoeff *v) { PyObject *strv; strv = deccoeff_str(v); if (strv == NULL) return NULL; return PyUnicode_FromFormat(CLASS_NAME "('%U')", deccoeff_str(v)); } static void deccoeff_dealloc(PyObject *v) { Py_TYPE(v)->tp_free(v); } /* increment or decrement a deccoeff; returns a new reference */ static deccoeff* decrement(deccoeff *v) { deccoeff *z; Py_ssize_t v_size, z_size, i; int lose_one; limb *v_start, *v_end; limb *v_limb, *z_limb; v_size = Py_SIZE(v); if (v_size == 0) { PyErr_SetString(PyExc_ValueError, "Cannot decrement 0"); return NULL; } /* skip zero limbs */ v_start = v->ob_limbs; v_end = v_start + v_size; v_limb = v_start; while (*v_limb == 0) v_limb++; lose_one = (v_limb == v_end-1) && (*v_limb == 1); z_size = v_size-lose_one; z = _deccoeff_new(z_size); if (z==NULL) return NULL; z_limb = z->ob_limbs; for (i=0; iob_limbs == z_size); assert(z_size == 0 || z->ob_limbs[z_size-1] != 0); return z; } static deccoeff* increment(deccoeff *v) { deccoeff *z; limb *v_start, *v_end, *v_limb, *z_limb; Py_ssize_t v_size, i; v_size = Py_SIZE(v); v_start = v->ob_limbs; v_end = v->ob_limbs + v_size; /* skip limbs equal to LIMB_MAX */ for (v_limb = v_start; v_limb < v_end; v_limb++) if (*v_limb != LIMB_MAX) break; /* correct number of limbs; no need to normalize XXX deal with possible overflow of v_size + 1 */ z = _deccoeff_new(v_size + (v_limb == v_end ? 1 : 0)); if (z==NULL) return NULL; /* fill in zeros */ z_limb = z->ob_limbs; for (i=0; i < v_limb-v_start; i++) *z_limb++ = 0; /* next limb = corresponding limb of v (if it exists) + 1 */ *z_limb++ = (v_limb < v_end ? *v_limb++ : 0) + 1; /* now copy over any remaining limbs */ while (v_limb < v_end) *z_limb++ = *v_limb++; return z; } static deccoeff* fix05(deccoeff *v) { deccoeff *z; Py_ssize_t v_size; limb *v_end, *v_limb, *z_limb; limb low_limb; v_size = Py_SIZE(v); if (v_size == 0) { Py_INCREF(deccoeff_one); return deccoeff_one; } /* if last digit of v is *not* 0 or 5, return v unchanged */ v_limb = v->ob_limbs; v_end = v_limb + v_size; low_limb = *v_limb++; if (low_limb % 5 != 0) { Py_INCREF(v); return v; } /* otherwise, increment it. */ z = _deccoeff_new(v_size); if (z==NULL) return NULL; z_limb = z->ob_limbs; *z_limb++ = ++low_limb; while (v_limb < v_end) { *z_limb++ = *v_limb++; } return z; } /* Compute v[a:b], for any integers a and b. The result can be interpreted as v % 10**b // 10**a. */ static deccoeff * get_slice(deccoeff *v, Py_ssize_t a, Py_ssize_t b) { Py_ssize_t alimbs, adigits, blimbs, bdigits, z_size, v_digits; limb *v_limb, *v_lastlimb, *z_limb, *v_start; limb divisor, multiplier, low, next; deccoeff *z; /* special cases, and reductions */ if ((Py_SIZE(v) == 0) || (b <= a) || (b <= 0)) return zero(); v_digits = deccoeff_length(v); if (v_digits == -1) return NULL; if (a >= v_digits) return zero(); if (b > v_digits) b = v_digits; /* a = alimbs * LIMB_DIGITS + adigits, 0 <= adigits < LIMB_DIGITS b = blimbs * LIMB_DIGITS + bdigits, 0 < bdigits <= LIMB_DIGITS limbs 'alimbs' through 'blimbs' (inclusive) contribute to result */ alimbs = a >= 0 ? a/LIMB_DIGITS : ~(~a/LIMB_DIGITS); adigits = a >= 0 ? a%LIMB_DIGITS : LIMB_DIGITS-1 - (~a%LIMB_DIGITS); blimbs = (b-1) / LIMB_DIGITS; bdigits = (b-1) % LIMB_DIGITS + 1; /* allocate new deccoeff with sufficient space for result */ z_size = (blimbs - alimbs) + (adigits < bdigits); z = _deccoeff_new(z_size); if (z==NULL) return NULL; z_limb = z->ob_limbs; /* fill in limbs of z */ v_start = v->ob_limbs; v_limb = v_start+alimbs; v_lastlimb = v_start + blimbs; divisor = powers_of_ten[adigits]; if (v_limb == v_lastlimb) { *z_limb++ = *v_limb++ % powers_of_ten[bdigits] / divisor; } else { if (v_limb < v_start) { v_limb++; while (v_limb < v_start) { v_limb++; *z_limb++ = 0; } low = 0; } else { low = *v_limb++ / divisor; } multiplier = powers_of_ten[LIMB_DIGITS-adigits]; while (v_limb < v_lastlimb) { next = *v_limb++; *z_limb++ = next % divisor * multiplier + low; low = next / divisor; } next = *v_limb++ % powers_of_ten[bdigits]; *z_limb++ = next % divisor * multiplier + low; if (adigits < bdigits) *z_limb++ = next/divisor; } /* check we've filled in the number of limbs that we expected to */ assert(z_limb - z->ob_limbs == z_size); return deccoeff_normalize(z); } enum ROUNDING_MODE { ROUND_DOWN, ROUND_UP, ROUND_HALF_DOWN, ROUND_HALF_EVEN, ROUND_HALF_UP, ROUND_05UP, ROUND_INVALID }; typedef enum ROUNDING_MODE rounding_mode; /* decide whether to round a nonnegative decimal number x up or not, based on three pieces of information about x: (1) half is true iff the fractional part of x is >= .5 (2) exact is true iff the fractional part of x is exactly 0 or 0.5 (3) units_digit is the units digit of x; i.e., floor(x) % 10. returns: -1 if the result should be floor(x) and floor(x) != x returns: 0 if x is an integer, and result is x exactly returns: 1 if result should be ceiling(x) */ static int rounding_decision(int half, int inexact, int units_digit, rounding_mode mode) { int increment; assert(half == 0 || half == 1); assert(inexact == 0 || inexact == 1); assert(0 <= units_digit && units_digit < 10); assert(0 <= mode && mode < ROUND_INVALID); if (mode == ROUND_DOWN) increment = 0; else if (mode == ROUND_UP) increment = half || inexact; else if (mode == ROUND_HALF_DOWN) increment = half && inexact; else if (mode == ROUND_HALF_EVEN) increment = half && (inexact || ((units_digit & 1) != 0)); else if (mode == ROUND_HALF_UP) increment = half; else if (mode == ROUND_05UP) increment = (half || inexact) && (units_digit == 0 || units_digit == 5); else assert(0); /* should never get here */ if (!increment && (inexact || half)) increment = -1; return increment; } static PyObject * rounds(deccoeff *v, PyObject *args) { PyObject *ans, *oshift; Py_ssize_t limbs, digits, i, count, current, z_size, v_size, shift; int inexact, half, increment; rounding_mode mode; limb divisor, quot, rem, top, next, *z_limb, multiplier; deccoeff *z; if (!PyArg_ParseTuple(args, "Oi", &oshift, &mode)) return NULL; if (mode < 0 || mode >= ROUND_INVALID) { PyErr_SetString(PyExc_ValueError, "invalid rounding mode"); return NULL; } /* oshift should be an integer */ if (!PyLong_Check(oshift)) { PyErr_SetString(PyExc_TypeError, "first argument of round operation " "should be an integer"); return NULL; } /* convert oshift to a Py_ssize_t; catch overflow error */ shift = PyLong_AsSsize_t(oshift); if (shift == -1 && PyErr_Occurred() != NULL) { if (PyErr_ExceptionMatches(PyExc_OverflowError)) { PyErr_Clear(); if (PyLong_Check(oshift) && Py_SIZE(oshift) >= 0) { /* shift by more than MAX_DIGITS */ increment = rounding_decision(0, Py_SIZE(v) > 0, 0, mode); return Py_BuildValue("OO", increment == 1 ? deccoeff_one : deccoeff_zero, increment == 0 ? Py_False : Py_True); } } else { /* propagate any unexpected exception */ return NULL; } } /* shift should be positive */ if (shift <= 0) { PyErr_SetString(PyExc_ValueError, "second operand to round operation must be positive"); return NULL; } v_size = Py_SIZE(v); /* write shift = limbs*LIMB_DIGITS + digits, with 1 <= digits <= LIMB_DIGITS. To do: optimize for the case when digits == LIMB_DIGITS. */ limbs = (shift-1)/LIMB_DIGITS; digits = (shift-1)%LIMB_DIGITS + 1; divisor = powers_of_ten[digits]; multiplier = powers_of_ten[LIMB_DIGITS-digits]; /* pull apart v->ob_limbs[limbs] */ top = GETLIMB(v, limbs); quot = top / divisor; rem = top % divisor; /* remainder gives us final value of inexact and half */ half = (rem >= divisor/2); inexact = (rem != (half ? divisor/2 : 0)); /* check all other limbs of v, starting with most significant */ if (!inexact) { for (i=(v_size < limbs ? v_size : limbs)-1; i>=0; i--) { if (GETLIMB(v, i) != 0) { inexact = 1; break; } } } /* get first limb of result; we need this to be able to decide whether to increment or not. */ top = GETLIMB(v, limbs+1); next = quot + multiplier*(top%divisor); quot = top/divisor; increment = rounding_decision(half, inexact, (int)(next % 10), mode); /* size of result, assuming that incrementing doesn't change it */ z_size = v_size - limbs; if (GETLIMB(v, v_size-1) < divisor) z_size--; if (z_size < 0) z_size = 0; if (increment==1) { /* first we just count the number of result limbs equal to LIMB_MAX. */ count = 0; current = limbs+2; while (next == LIMB_MAX) { count += 1; top = GETLIMB(v, current); current++; next = quot + multiplier*(top%divisor); quot = top/divisor; } /* Under what circumstances does incrementing change the size? Only if *all* the result limbs are equal to LIMB_MAX. */ if (count == z_size) z_size++; assert(z_size > 0); /* now we know how many limbs the result is: allocate space */ z = _deccoeff_new(z_size); if (z==NULL) return NULL; /* put zeros in */ z_limb = z->ob_limbs; for (i=0; iob_limbs == z_size); } else { /* not incrementing; merely copying. How many limbs are in the result?. There were v_size to begin with; we've lost limbs of them, and if the top limb is less than divisor we'll lose one more. */ if (z_size > 0) { /* put in first limb of result */ z = _deccoeff_new(z_size); if (z==NULL) return NULL; z_limb = z->ob_limbs; *z_limb++ = next; current = limbs + 2; while (current < v_size) { top = GETLIMB(v, current); current++; next = quot + multiplier*(top%divisor); quot = top/divisor; *z_limb++ = next; } /* last limb of v; quot is currently the top limb of v, divided by divisor */ if (quot != 0) { *z_limb++ = quot; } assert(z_limb - z->ob_limbs == z_size); } else { assert(next == 0); z = deccoeff_zero; Py_INCREF(z); } } assert(Py_SIZE(z) == 0 || z->ob_limbs[Py_SIZE(z)-1] != 0); ans = Py_BuildValue("OO", z, increment==0 ? Py_False : Py_True); Py_DECREF(z); return ans; } /* count the number of trailing zeros on a nonzero deccoeff */ static PyObject * count_zeros(deccoeff *v) { Py_ssize_t ans, v_size; limb *v_limb; limb l; v_size = Py_SIZE(v); if (v_size == 0) { PyErr_SetString(PyExc_ValueError, "argument to count_zeros should be positive"); return NULL; } /* count zero_limbs */ v_limb = v->ob_limbs; while (*v_limb == 0) v_limb++; ans = (v_limb - v->ob_limbs)*LIMB_DIGITS; l = *v_limb; assert(l != 0); while (l % 10 == 0) { ans++; l /= 10; } return Py_BuildValue("n", ans); } #define HASH_SHIFT 13 #define HASH_MASK ((1<ob_limbs; v_top = v_start + Py_SIZE(v); x = HASH_START; while (v_top > v_start) { x = ((x & ~HASH_MASK) >> HASH_SHIFT) | ((x & HASH_MASK) << (HASH_BITS-HASH_SHIFT)); x ^= *--v_top; } y = (long)x; return (y == -1) ? -2 : y; } static deccoeff* deccoeff_subscript(deccoeff *v, PyObject* sslice) { Py_ssize_t start, stop; PySliceObject *slice; if (PyIndex_Check(sslice)) { start = PyNumber_AsSsize_t(sslice, PyExc_IndexError); if (start==-1 && PyErr_Occurred()) return NULL; stop = start+1; } else if (PySlice_Check(sslice)) { slice = (PySliceObject*)sslice; if (slice->step != Py_None) { PyErr_SetString(PyExc_ValueError, "Please don't specify a step"); return NULL; } if (slice->start == Py_None) start = 0; else if (!_PyEval_SliceIndex(slice->start, &start)) { return NULL; } if (slice->stop == Py_None) { stop = Py_SIZE(v) ? deccoeff_length(v) : 0; } else if (!_PyEval_SliceIndex(slice->stop, &stop)) { return NULL; } } else { PyErr_SetString(PyExc_TypeError, "Expected an integer or a slice"); return NULL; } return get_slice(v, start, stop); } static PyMethodDef deccoeff_methods[] = { {"rounds", (PyCFunction)rounds, METH_VARARGS, "rounds(shift, rounding_mode)"}, {"increment", (PyCFunction)increment, METH_NOARGS, "increment"}, {"decrement", (PyCFunction)decrement, METH_NOARGS, "increment"}, {"from_int", (PyCFunction)deccoeff_from_PyLong, METH_O | METH_STATIC, "create from an integer"}, {"count_zeros", (PyCFunction)count_zeros, METH_NOARGS, "count trailing zeros on a nonzero deccoeff"}, {"fix05", (PyCFunction)fix05, METH_NOARGS, "if last digit is 0 or 5, increment it"}, {NULL, NULL} }; static PyMappingMethods deccoeff_as_mapping = { (lenfunc)deccoeff_length, /*mp_length*/ (binaryfunc)deccoeff_subscript, /*mp_subscript*/ 0,/*mp_ass_subscript*/ }; static PyNumberMethods deccoeff_as_number = { (binaryfunc) deccoeff_add, /*nb_add*/ (binaryfunc) deccoeff_sub, /*nb_subtract*/ (binaryfunc) deccoeff_mul, /*nb_multiply*/ (binaryfunc) deccoeff_remainder, /*nb_remainder*/ (binaryfunc) deccoeff_divmod, /*nb_divmod*/ (ternaryfunc) deccoeff_pow, /*nb_power*/ 0,/*nb_negative*/ 0,/*tp_positive*/ 0,/*tp_absolute*/ (inquiry)deccoeff_nonzero, /*tp_bool*/ 0,/*nb_invert*/ (binaryfunc)deccoeff_lshift, /*nb_lshift*/ (binaryfunc)deccoeff_rshift, /*nb_rshift*/ (binaryfunc)deccoeff_and, /*nb_and*/ (binaryfunc)deccoeff_xor, /*nb_xor*/ (binaryfunc)deccoeff_or, /*nb_or*/ 0,/*nb_reserved*/ 0,/*nb_int*/ 0,/*nb_long*/ 0,/*nb_float*/ 0,/*nb_oct*/ 0,/*nb_hex*/ 0, /* nb_inplace_add */ 0, /* nb_inplace_subtract */ 0, /* nb_inplace_multiply */ 0, /* nb_inplace_remainder */ 0, /* nb_inplace_power */ 0, /* nb_inplace_lshift */ 0, /* nb_inplace_rshift */ 0, /* nb_inplace_and */ 0, /* nb_inplace_xor */ 0, /* nb_inplace_or */ (binaryfunc)deccoeff_floor_divide, /* nb_floor_divide */ 0, /* nb_true_divide */ 0, /* nb_inplace_floor_divide */ 0, /* nb_inplace_true_divide */ 0, /* nb_index */ }; static PyTypeObject deccoeff_DeccoeffType = { PyVarObject_HEAD_INIT(&PyType_Type, 0) MODULE_NAME "." CLASS_NAME, /* tp_name */ sizeof(deccoeff), /* tp_basicsize */ sizeof(limb), /* tp_itemsize */ deccoeff_dealloc, /* tp_dealloc */ 0, /* tp_print */ 0, /* tp_getattr */ 0, /* tp_setattr */ 0, /* tp_compare */ (reprfunc)deccoeff_repr, /* tp_repr */ &deccoeff_as_number, /* tp_as_number */ 0, /* tp_as_sequence */ &deccoeff_as_mapping, /* tp_as_mapping */ (hashfunc)deccoeff_hash, /* tp_hash */ 0, /* tp_call */ (reprfunc)deccoeff_str, /* tp_str */ 0, /* tp_getattro */ 0, /* tp_setattro */ 0, /* tp_as_buffer */ Py_TPFLAGS_DEFAULT, /* tp_flags */ "Decimal integers", /* tp_doc */ 0, /* tp_traverse */ 0, /* tp_clear */ (richcmpfunc)deccoeff_richcompare, /* tp_richcompare */ 0, /* tp_weaklistoffset */ 0, /* tp_iter */ 0, /* tp_iternext */ deccoeff_methods, /* tp_methods */ 0, /* tp_members */ 0, /* tp_getset */ 0, /* tp_base */ 0, /* tp_dict */ 0, /* tp_descr_get */ 0, /* tp_descr_set */ 0, /* tp_dictoffset */ 0, /* tp_init */ 0, /* tp_alloc */ deccoeff_new, /* tp_new */ PyObject_Del, /* tp_free */ }; static PyMethodDef deccoeff_module_methods[] = { {NULL, NULL} }; #define ADD_CONST(m, name) \ if (PyModule_AddIntConstant(m, #name, name) < 0) return PyMODINIT_FUNC initdeccoeff(void) { PyObject *m; /* a module object */ if (PyType_Ready(&deccoeff_DeccoeffType) < 0) return; m = Py_InitModule3(MODULE_NAME, deccoeff_module_methods, "Decimal integer implementation."); if (m == NULL) return; Py_INCREF(&deccoeff_DeccoeffType); PyModule_AddObject(m, CLASS_NAME, (PyObject *) &deccoeff_DeccoeffType); ADD_CONST(m, ROUND_DOWN); ADD_CONST(m, ROUND_UP); ADD_CONST(m, ROUND_HALF_DOWN); ADD_CONST(m, ROUND_HALF_UP); ADD_CONST(m, ROUND_HALF_EVEN); ADD_CONST(m, ROUND_05UP); deccoeff_zero = deccoeff_from_ulong((unsigned long)0); if (deccoeff_zero == NULL) return; /* XXX what happens if there's an error here? How does the caller know? */ deccoeff_one = deccoeff_from_ulong((unsigned long)1); if (deccoeff_one == NULL) return; deccoeff_two = deccoeff_from_ulong((unsigned long)2); if (deccoeff_two == NULL) return; deccoeff_PyLong_BASE = deccoeff_from_ulong((unsigned long)PyLong_BASE); if (deccoeff_PyLong_BASE == NULL) return; /* XXX do we need to deallocate deccoeff_zero and deccoeff_PyLong_BASE somewhere? Do multiple imports leak references? */ }