[Python-checkins] r50663 - sandbox/trunk/decimal-c/_decimal.c
mateusz.rukowicz
python-checkins at python.org
Sat Jul 15 08:20:21 CEST 2006
Author: mateusz.rukowicz
Date: Sat Jul 15 08:20:19 2006
New Revision: 50663
Modified:
sandbox/trunk/decimal-c/_decimal.c
Log:
Added ln function and x^y with non-integer y (this still need some fixing). Improved efficency of dividing a little bit.
Modified: sandbox/trunk/decimal-c/_decimal.c
==============================================================================
--- sandbox/trunk/decimal-c/_decimal.c (original)
+++ sandbox/trunk/decimal-c/_decimal.c Sat Jul 15 08:20:19 2006
@@ -241,7 +241,7 @@
static long
_limb_get_digit(long *self, long ndigits, long i)
{
- if(i >= ndigits)
+ if(i >= ndigits || i < 0)
return 0;
long pos = ndigits - i - 1;
long limb = pos / LOG;
@@ -397,6 +397,27 @@
}
+
+static long
+_limb_multiply_int(long *first, long flimbs, long second, long *out) {
+ long i;
+ for (i=0;i<flimbs;i++)
+ out[i] = first[i] * second;
+ out[flimbs] = 0;
+
+ for (i=0;i<=flimbs;i++) {
+ if (out[i] >= BASE) {
+ assert (i+1 <=flimbs);
+ out[i+1] += out[i]/BASE;
+ out[i] %= BASE;
+ }
+ }
+
+ if (out[flimbs])
+ return flimbs+1;
+ return flimbs;
+
+}
/* temporary solution, will be speeded up */
static long
_limb_multiply_core(long *first, long flimbs, long *second, long slimbs, long *out)
@@ -473,12 +494,45 @@
while(1) {
long candidate = 0;
long cmp;
+ long up, down;
+ up = BASE;
+ down = 0;
+
+ while(1) {
+ long tmp[slimbs+1];
+ long tmp_limbs;
+ long diff;
+ long mid;
+ int bla;
+ diff = up - down;
+ mid = down + diff/2;
+
+ tmp_limbs = _limb_multiply_int(second, slimbs, mid, tmp);
+ cmp = _limb_compare_un(rest, rlimbs, tmp, tmp_limbs);
+ if (cmp == 0) {
+ up = mid + 1;
+ down = mid;
+ }
+ if (cmp == -1) {
+ up = mid;
+ }
+ if (cmp == 1) {
+ down = mid;
+ }
+
+ if (down == up-1) {
+ tmp_limbs = _limb_multiply_int(second, slimbs, down, tmp);
+ rlimbs = _limb_sub_sl(rest, rlimbs, tmp, tmp_limbs);
+ candidate = down;
+ break;
+ }
+ }
- while ((cmp = _limb_compare_un(rest, rlimbs, second, slimbs)) >= 0) {
+ /* while ((cmp = _limb_compare_un(rest, rlimbs, second, slimbs)) >= 0) {
candidate ++;
rlimbs = _limb_sub_sl(rest, rlimbs, second, slimbs);
}
-
+*/
if (candidate)
is_significant = 1;
@@ -1074,6 +1128,8 @@
#define SIGN_POSSNAN 6
#define SIGN_NEGSNAN 7
+/* Bigest exponent in transcendetial functions */
+#define MAX_MATH 999999
/* Rounding constants */
/* for context->rounding */
@@ -1900,7 +1956,7 @@
ans->digits[1] = 1;
ans->limbs[0] = 1;
digits = 1;
- } else {
+ } else { /* SLOW */
ans = _NEW_decimalobj(self->ob_size+1, self->sign, self->exp);
if (!ans)
return NULL;
@@ -4507,12 +4563,12 @@
ans = tmp;
}
+ ctx->rounding = rounding;
{
decimalobject *fixed = _decimal_fix(ans, ctx);
Py_DECREF(ans);
ans = fixed;
}
- ctx->rounding = rounding;
return ans;
}
@@ -4525,14 +4581,261 @@
if (!ctx)
return NULL;
- if (ctx->prec > 999999 ||
- exp_g_i(ctx->Emax, 999999) ||
- exp_l_i(ctx->Emin,-999999))
+ if (ctx->prec > MAX_MATH ||
+ exp_g_i(ctx->Emax, MAX_MATH) ||
+ exp_l_i(ctx->Emin,-1 * MAX_MATH))
handle_InvalidContext(self->ob_type, ctx, NULL);
return _do_decimal_exponent(self, ctx);
}
+int ln_lookup[] = {9016, 8652, 8316, 8008, 7724, 7456, 7208,
+ 6972, 6748, 6540, 6340, 6148, 5968, 5792, 5628, 5464, 5312,
+ 5164, 5020, 4884, 4748, 4620, 4496, 4376, 4256, 4144, 4032,
+ 39233, 38181, 37157, 36157, 35181, 34229, 33297, 32389, 31501, 30629,
+ 29777, 28945, 28129, 27329, 26545, 25777, 25021, 24281, 23553, 22837,
+ 22137, 21445, 20769, 20101, 19445, 18801, 18165, 17541, 16925, 16321,
+ 15721, 15133, 14553, 13985, 13421, 12865, 12317, 11777, 11241, 10717,
+ 10197, 9685, 9177, 8677, 8185, 7697, 7213, 6737, 6269, 5801,
+ 5341, 4889, 4437, 39930, 35534, 31186, 26886, 22630, 18418, 14254,
+ 10130, 6046, 20055};
+
+static PyObject *
+_do_decimal__ln(decimalobject *self, contextobject *ctx) {
+
+ decimalobject *ans = 0, *b = 0, *tmp, *one = 0;
+ contextobject *ctx1 = 0, *ctx2 = 0;
+ PyObject *flags;
+ long precision;
+ int rounding;
+ int clamped;
+ long prec, cp;
+ int t;
+ if (!ctx)
+ ctx = getcontext();
+ if (!ctx)
+ return NULL;
+
+ if (ISSPECIAL(self)) {
+ decimalobject *nan;
+ if (_check_nans(self, NULL, ctx, &nan) != 0)
+ return nan;
+
+ /*-inf -> error
+ * inf -> self */
+ if (ISINF(self)) {
+ if (self->sign &1)
+ return handle_InvalidOperation(self->ob_type, ctx, "ln(-inf)", NULL);
+ else
+ return _decimal_get_copy(self);
+ }
+ }
+
+ if (!decimal_nonzero(self)) {
+ ans = _NEW_decimalobj(1, SIGN_NEGINF, exp_from_i(0));
+ return ans;
+ }
+
+ if (self->sign &1) {
+ return handle_InvalidOperation(self->ob_type, ctx, "ln(-x)", NULL);
+ }
+
+ one = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!one)
+ return NULL;
+ one->limbs[0] = 1;
+
+ ctx1 = context_copy(ctx);
+ if (!ctx1)
+ goto err;
+
+ flags = context_ignore_all_flags(ctx1);
+ if (!flags)
+ goto err;
+ Py_DECREF(flags);
+
+
+
+ ctx1->prec = 16;
+
+ prec = (self->ob_size > ctx->prec ? self->ob_size : ctx->prec) + 2;
+ prec = prec > 9 ? prec : 9;
+
+ ans = decimal_from_long(self->ob_type, exp_to_i(ADJUSTED(self)) + 1);
+ if (!ans)
+ goto err;
+ b = decimal_from_long(self->ob_type, 2302585);
+ if (!b)
+ goto err;
+ b->exp = exp_from_i(-6);
+
+
+ tmp = _do_decimal_multiply(ans, b, ctx1);
+ if (!tmp)
+ goto err;
+
+ Py_DECREF(ans);
+ ans = tmp;
+
+ t = _limb_get_digit(self->limbs, self->ob_size, 0);
+ if (self->ob_size > 1) {
+ t *= 10;
+ t += _limb_get_digit(self->limbs, self->ob_size, 1);
+ }
+ if (t<10) t *= 10;
+
+ t = ln_lookup[t-10];
+
+ Py_DECREF(b);
+ b = decimal_from_long(self->ob_type, t >> 2);
+ if (!b)
+ goto err;
+ b->exp = exp_from_i(-(t&3) - 3);
+ b->sign = 1;
+
+ rounding = ctx1->rounding;
+ ctx1->rounding = ROUND_HALF_EVEN;
+ tmp = _do_decimal_add(ans, b, ctx1);
+ if (!tmp)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp;
+
+ ctx1->prec = ctx->prec;
+ ctx1->clamp = 0;
+ ctx2 = context_shallow_copy(ctx1);
+ if (!ctx2)
+ goto err;
+ ctx2->Emax = exp_from_i(2 * MAX_MATH);
+ ctx2->Emin = exp_from_i(-2 * MAX_MATH);
+
+ cp = 9; /* precision we'll get after every iteration */
+ ctx1->prec = cp;
+ ctx2->prec = cp + self->ob_size;
+
+
+ while (1) {
+ ans->sign ^= 1;
+ Py_DECREF(b);
+ b = _do_decimal_exponent(ans, ctx2);
+ if (!b)
+ goto err;
+ ans->sign ^= 1;
+
+ tmp = _do_decimal_multiply(b, self, ctx2);
+ if (!tmp)
+ goto err;
+ Py_DECREF(b);
+ b = tmp;
+
+ tmp = _do_decimal_subtract(b, one, ctx2);
+ if (!tmp)
+ goto err;
+ Py_DECREF(b);
+ b = tmp;
+
+ /* ADJUSTED(ans) >= ADJUSTED(b) + ctx->prec + 1 */
+ if (!decimal_nonzero(b) || exp_ge(ADJUSTED(ans),
+ exp_add_i(ADJUSTED(b), ctx->prec + 1))) {
+ if (ans->ob_size == prec)
+ break;
+
+ if (!decimal_nonzero(ans)) {
+ if (!_do_real_decimal_compare(self, one, ctx1)) {
+ ans->exp = exp_from_i(0);
+ break;
+ }
+ else {
+ if (handle_Rounded(ctx, NULL))
+ goto err;
+ if (handle_Inexact(ctx, NULL))
+ goto err;
+
+ }
+ }
+
+ /* to make addition easier */
+ if (!decimal_nonzero(b))
+ b->exp = exp_sub_i(ans->exp, prec);
+ }
+
+ tmp = _do_decimal_add(ans, b, ctx1);
+ if (!tmp)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp;
+
+ /* we're done */
+ if (cp == prec)
+ continue;
+
+ cp *= 2;
+ if (cp > prec)
+ cp = prec;
+ ctx1->prec = cp;
+ ctx2->prec = cp + self->ob_size;
+ }
+
+ rounding = ctx->rounding;
+ ctx->rounding = ROUND_HALF_EVEN;
+ /* we add extra 1 at the end to make proper rounding */
+ if (decimal_nonzero(ans) && !ISSPECIAL(ans)) {
+ int i;
+
+
+ tmp = _NEW_decimalobj(ans->ob_size + LOG, ans->sign, exp_sub_i(ans->exp, LOG));
+ if (!tmp) {
+ ctx->rounding = rounding;
+ goto err;
+ }
+
+ for (i = 0; i < ans->limb_count ; i++){
+ tmp->limbs[i+1] = ans->limbs[i];
+ }
+ tmp->limbs[0] = 1;
+
+ Py_DECREF(ans);
+ ans = tmp;
+ }
+
+ ctx->rounding = rounding;
+ tmp = _decimal_fix(ans, ctx);
+
+ if (!tmp)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp;
+
+ Py_DECREF(b);
+ Py_DECREF(ctx1);
+ Py_DECREF(ctx2);
+ Py_DECREF(one);
+ return ans;
+err:
+ Py_DECREF(one);
+ Py_XDECREF(ans);
+ Py_XDECREF(b);
+ Py_XDECREF(ctx1);
+ Py_XDECREF(ctx2);
+
+ return NULL;
+}
+
+DECIMAL_UNARY_FUNC(_ln);
+
+static PyObject *
+_do_decimal_ln(decimalobject *self, contextobject *ctx) {
+ return _do_decimal__ln(self, ctx);
+}
+
+DECIMAL_UNARY_FUNC(ln);
+
static PyMethodDef decimal_methods[] = {
+ {"ln", (PyCFunction)decimal_ln,
+ METH_VARARGS | METH_KEYWORDS,
+ PyDoc_STR("TODO")},
+ {"_ln", (PyCFunction)decimal__ln,
+ METH_VARARGS | METH_KEYWORDS,
+ PyDoc_STR("TODO")},
{"exp", (PyCFunction)decimal_exponent,
METH_VARARGS | METH_KEYWORDS,
PyDoc_STR("TODO")},
@@ -5121,20 +5424,126 @@
int mod = modulo != Py_None;
long firstprec = ctx->prec;
int cmp;
+ int use_exp = 0; /* when we should use exp/ln method */
+
+ if (!ctx)
+ ctx = getcontext();
+ if (!ctx)
+ return NULL;
+
+// if ((ISINF(other) || exp_g_i(ADJUSTED(other), 8)) && decimal_nonzero(other) ) {
+// return handle_InvalidOperation(self->ob_type, ctx, "x ** INF", NULL);
+// }
+
+ if ((exp_g_i(ADJUSTED(other), 8)) && decimal_nonzero(other)) {
+ use_exp = 1;
+ }
+ if(0)
+ if (exp_eq_i(ADJUSTED(other), 9) && _limb_get_digit(other->limbs, other->ob_size, 0) >=2)
+ use_exp = 1;
- if (ISINF(other) || exp_g_i(ADJUSTED(other), 8)) {
- return handle_InvalidOperation(self->ob_type, ctx, "x ** INF", NULL);
- }
-
if (ISSPECIAL(self) || ISSPECIAL(other)) {
decimalobject *nan;
if (_check_nans(self, other, ctx, &nan))
return nan;
}
+ if (ISINF(other)) {
+ int cmp;
+
+ if (!decimal_nonzero(self)) {
+ if (other->sign &1) {
+ return _NEW_decimalobj(1, SIGN_POSINF, exp_from_i(0));
+ }
+ else {
+ ret = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ret)
+ return NULL;
+ ret->limbs[0] = 0;
+ return ret;
+ }
+ }
+
+ if (self->sign &1)
+ return handle_InvalidOperation(self->ob_type, ctx, "-x ** [+-]INF", NULL);
+
+ {
+ contextobject *ctx2;
+ decimalobject *one;
+ ctx2 = context_copy(ctx);
+ if (!ctx2)
+ return NULL;
+
+
+ PyObject *flags = context_ignore_all_flags(ctx2);
+ if (!flags){
+ Py_DECREF(ctx2);
+ return NULL;
+ }
+ Py_DECREF(flags);
+
+ one = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!one) {
+ Py_DECREF(ctx2);
+ return NULL;
+ }
+ one->limbs[0] = 1;
+
+ cmp = _do_real_decimal_compare(self, one, ctx2);
+ Py_DECREF(ctx2);
+ Py_DECREF(one);
+ if (PyErr_Occurred())
+ return NULL;
+ }
+
+ if (cmp == 0) {
+ int i;
+ if (handle_Inexact(ctx, NULL))
+ return NULL;
+ if (handle_Rounded(ctx, NULL))
+ return NULL;
+
+ ret = _NEW_decimalobj(ctx->prec, 0, exp_from_i(-ctx->prec+1));
+ if (!ret)
+ return NULL;
+ for (i=0;i<ret->limb_count;i++)
+ ret->limbs[i] = 0;
+
+ {
+ long mult = 1;
+ long where;
+ long limb = (ctx->prec-1) / LOG;
+ where = (ctx->prec-1) % LOG;
+ while(where --) mult *= 10;
+ ret->limbs[limb] = mult;
+ }
+
+ return ret;
+ }
+ /* cmp == 1 for self > 1
+ * cmp == 0 for self < 1*/
+ cmp += 1;
+ cmp /= 2;
+
+ /* if self > 1 self^inf = inf, self^-inf = 0
+ * if self < 1 self^inf = 0, self^-inf = inf */
+
+ /* inf */
+ if (cmp ^ (other->sign &1)) {
+ return _NEW_decimalobj(1, SIGN_POSINF, exp_from_i(0));
+ }
+ /* 0 */
+ else {
+ ret = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ret)
+ return NULL;
+ ret->limbs[0] = 0;
+ return ret;
+ }
+ }
+
if (!_decimal_isint(other)) {
- return handle_InvalidOperation(self->ob_type, ctx,
- "x ** (non-integer)", NULL);
+ use_exp = 1;
}
if (!decimal_nonzero(self) && !decimal_nonzero(other)) {
@@ -5147,6 +5556,7 @@
return ret;
}
+ if (!use_exp)
{
PyObject *tmp = decimal_int(other);
if (!tmp)
@@ -5160,24 +5570,141 @@
}
}
+ /* all we really need here is n%2 - for special cases */
+ else {
+ /* where = other->ob_size - other->exp - 1 */
+ exp_t where = exp_from_i(other->ob_size);
+ exp_inp_sub(&where, exp_add_i(other->exp, 1));
+
+ n = _limb_get_digit(other->limbs, other->ob_size, exp_to_i(where));
+ }
sign = (self->sign&1) && n&1;
+ if (!decimal_nonzero(self)) {
+ if (other->sign &1) {
+ ret = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ret)
+ return NULL;
+ ret->sign = sign ? SIGN_NEGINF : SIGN_POSINF;
+ return ret;
+ }
+ else {
+ ret = _NEW_decimalobj(1, sign, exp_from_i(0));
+ if (!ret)
+ return NULL;
+ ret->limbs[0] = 0;
+ return ret;
+ }
+ }
if (ISINF(self)) {
if (mod) {
return handle_InvalidOperation(self->ob_type, ctx,
"INF % x", NULL);
}
+ if (!_decimal_isint(other) && self->sign&1)
+ return handle_InvalidOperation(self->ob_type, ctx, "INF ** -(non-int)", NULL);
ret = _NEW_decimalobj(1, sign, exp_from_i(0));
ret->limbs[0] = 0;
- if (n > 0)
+ if (decimal_nonzero(other) && !(other->sign &1))
ret->sign = sign ? SIGN_NEGINF : SIGN_POSINF;
return ret;
}
+ /* non-integer case */
+ /* we calculate it using exp(ln(self) * other) */
+ if (use_exp) {
+ decimalobject *tmp;
+ contextobject *ctx2;
+ ctx2 = context_shallow_copy(ctx);
+ if (!ctx2)
+ return NULL;
+
+ ctx2->Emax = exp_from_i(MAX_MATH);
+ ctx2->Emin = exp_from_i(-MAX_MATH);
+ ctx2->clamp = 0;
+ ctx2->rounding = ctx->rounding;
+ ctx2->rounding_dec = ALWAYS_ROUND;
+ /* we take context precision or size of self if greater
+ * and add some constant */
+ ctx2->prec = ctx->prec > self->ob_size ? ctx->prec : self->ob_size;
+ ctx2->prec += 10;
+
+ ret = _do_decimal__ln(self, ctx2);
+ if (!ret) {
+ Py_DECREF(ctx2);
+ return NULL;
+ }
+
+ if (!decimal_nonzero(ret)) {
+ int i;
+ Py_DECREF(ret);
+ Py_DECREF(ctx2);
+ if (!_decimal_isint(other))
+ {
+ ret = _NEW_decimalobj(ctx->prec, 0, exp_from_i(-ctx->prec + 1));
+ if (!ret)
+ return NULL;
+
+ if (handle_Inexact(ctx, NULL))
+ return NULL;
+ if (handle_Rounded(ctx, NULL))
+ return NULL;
+
+ for (i = 0; i < ret -> limb_count; i++) {
+ ret->limbs[i] = 0;
+ }
+ {
+ long mult = 1;
+ long limb;
+ long where = (ctx->prec - 1)% LOG;
+ limb = (ctx->prec - 1)/LOG;
+ while(where--) mult *= 10;
+ ret->limbs[limb] = mult;
+ }
+
+ return ret;
+ }
+ else {
+ ret = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ret)
+ return NULL;
+ ret->limbs[0] = 1;
+ return ret;
+ }
+ }
+
+ tmp = _do_decimal_multiply(ret, other, ctx2);
+ Py_DECREF(ret);
+ if (!tmp) {
+ Py_DECREF(ctx2);
+ return NULL;
+ }
+ ret = tmp;
+
+ tmp = _do_decimal_exponent(ret, ctx2);
+ Py_DECREF(ret);
+ Py_DECREF(ctx2);
+
+ if (!tmp)
+ return NULL;
+
+ ret = tmp;
+
+ if (ctx->rounding_dec == ALWAYS_ROUND) {
+ tmp = _decimal_fix(ret, ctx);
+ Py_DECREF(ret);
+ if (!tmp)
+ return NULL;
+ ret = tmp;
+ }
+
+ return ret;
+ }
+
/* XXX temporary solution */
{
#ifdef BIG_EXP
@@ -5222,7 +5749,11 @@
t/=10;
}
}
+ /* when n < 0 we need to extend precision - then every tests pass =] */
+ if (n<0)
+ ctx->prec += 1;
+ /* TODO shouldn't it be Invalid_Context ? */
if (!mod && exp_l_i(PyDecimal_DefaultContext->Emax, ctx->prec)) {
ctx->prec = firstprec;
return handle_Overflow(self->ob_type, ctx, "Too much precision", sign);
@@ -5234,6 +5765,7 @@
if(!val || !mul) {
Py_XDECREF(mul);
Py_XDECREF(val);
+ ctx->prec = firstprec;
return NULL;
}
@@ -5244,6 +5776,7 @@
if (!tmp) {
Py_DECREF(mul);
Py_DECREF(val);
+ ctx->prec = firstprec;
return NULL;
}
@@ -5615,6 +6148,7 @@
* whether the fractional part consists only of zeroes. */
i = exp_from_i(d->ob_size-1);
cmp = exp_add_i(d->exp, d->ob_size);
+ /* TODO bail out when we go out of bound */
for (; exp_ge(i, cmp); exp_dec(&i)) {
/* if (d->digits[i] > 0)*/
if(_limb_get_digit(d->limbs, d->ob_size, exp_to_i(i)) > 0)
@@ -6714,6 +7248,7 @@
CONTEXT_UNARY_FUNC(sqrt, sqrt)
CONTEXT_UNARY_FUNC(to_eng_string, to_eng_string)
CONTEXT_UNARY_FUNC(exp, exp)
+CONTEXT_UNARY_FUNC(ln, ln);
/* Unfortunately, the following methods are non-standard and can't
@@ -7238,6 +7773,8 @@
METH_O},
{"exp", (PyCFunction)context_exp,
METH_O},
+ {"ln", (PyCFunction)context_ln,
+ METH_O},
{"subtract", (PyCFunction)context_subtract,
METH_VARARGS},
{"to_eng_string", (PyCFunction)context_to_eng_string,
More information about the Python-checkins
mailing list