[Python-checkins] r50628 - sandbox/trunk/decimal-c/_decimal.c
mateusz.rukowicz
python-checkins at python.org
Fri Jul 14 03:27:07 CEST 2006
Author: mateusz.rukowicz
Date: Fri Jul 14 03:27:05 2006
New Revision: 50628
Modified:
sandbox/trunk/decimal-c/_decimal.c
Log:
Added exponention function. Some updates to newest standard.
Modified: sandbox/trunk/decimal-c/_decimal.c
==============================================================================
--- sandbox/trunk/decimal-c/_decimal.c (original)
+++ sandbox/trunk/decimal-c/_decimal.c Fri Jul 14 03:27:05 2006
@@ -1203,6 +1203,7 @@
static decimalobject *_do_decimal_multiply(decimalobject *, decimalobject *, contextobject *);
static decimalobject *_do_decimal_subtract(decimalobject *, decimalobject *, contextobject *);
static PyObject *context_ignore_flags(contextobject *self, PyObject *args);
+static decimalobject *_do_decimal_power(decimalobject *, decimalobject *, decimalobject *, contextobject *);
/* Exception handlers *********************************************************/
@@ -1992,6 +1993,10 @@
ans = _decimal_rescale(self, Etiny, ctx, -1, 1);
if (!ans)
return NULL;
+ if (!decimal_nonzero(ans)) {
+ if (handle_Clamped(ctx, NULL) != 0)
+ goto err;
+ }
if (handle_Subnormal(ctx, NULL) != 0)
goto err;
if (_is_flag_set(ctx->flags, S_INEXACT)) {
@@ -2502,6 +2507,8 @@
}
dup = _decimal_fix(self, ctx);
+ if (!dup)
+ return NULL;
if (ISINF(dup))
return dup;
@@ -4228,7 +4235,310 @@
return (PyObject*) _do_decimal__divide(self, other, divmod, ctx);
}
+
+/* The idea is as follows, we need to calculate e^x (self = x),
+ * we use series e^x = x^0/0! + x^1/1! + x^2/2! ...
+ * But first, we divide x, by integer, so that x < 1 - this results
+ * in much faster calculation. So we have e^(x/a), where a = 10^h,
+ * when we got e^(x/a) we raise it to power a. H is at most 8. */
+static PyObject *
+_do_decimal_exponent(decimalobject *self, contextobject *ctx) {
+ decimalobject *ans;
+ long h;
+ int clamped;
+ long precision, prec;
+ int rounding;
+ exp_t new_exp;
+
+ if (!ctx)
+ ctx = getcontext();
+
+ if (!ctx)
+ return NULL;
+
+ if (ISSPECIAL(self)) {
+ decimalobject *nan;
+
+ if (_check_nans(self, NULL, ctx, &nan))
+ return nan;
+
+ /* -Infinity -> +0 */
+ /* Infinity -> self */
+ if (ISINF(self)) {
+ if (self->sign == SIGN_NEGINF) {
+ ans = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ans)
+ return NULL;
+
+ ans->limbs[0] = 0;
+ return ans;
+ }
+ else {
+ ans = _decimal_get_copy(self);
+ return ans;
+ }
+ }
+ }
+
+ if (!decimal_nonzero(self)) {
+ ans = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ans)
+ return NULL;
+
+ ans->limbs[0] = 1;
+ return ans;
+ }
+
+ h = exp_to_i(ADJUSTED(self)) + 1;
+
+ precision = ctx->prec;
+ clamped = ctx->clamp;
+ rounding = ctx->rounding;
+ ctx->rounding = ROUND_HALF_EVEN;
+ /* this will be changed */
+ if (h > 8) {
+ h = 8;
+ prec = 9;
+ ans = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ans)
+ return NULL;
+ ans->limbs[0] = 2;
+
+ if (self->sign&1) {
+ ans->exp = exp_from_i(-2);
+ }
+ }
+ else {
+ decimalobject *t; /* that's what we add */
+ decimalobject *d; /* divisor */
+ decimalobject *rescaled;
+ contextobject *ctx2;
+ d = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!d)
+ return NULL;
+ d->limbs[0] = 2;
+
+ if (h < 0)
+ h = 0;
+
+ /* we increase h, so self is smaller and series go faster */
+ if (self->ob_size > 8 && h < 8)
+ h ++;
+
+ new_exp = exp_sub_i(self->exp, h);
+
+ if (h) {
+ rescaled = _decimal_get_copy(self);
+ if (!rescaled) {
+ Py_DECREF(d);
+ return NULL;
+ }
+
+ rescaled->exp = new_exp;
+ }
+ else {
+ Py_INCREF(self);
+ rescaled = self;
+ }
+
+ prec = (ctx->prec > self->ob_size ? ctx->prec : self->ob_size)
+ + h + 2;
+
+ ctx2 = context_copy(ctx);
+ if (!ctx2) {
+ Py_DECREF(rescaled);
+ Py_DECREF(d);
+ return NULL;
+ }
+ t = _decimal_get_copy(rescaled);
+ if (!t) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ return NULL;
+ }
+
+ ctx2->Emin = exp_from_i(-999999999);
+ ctx->clamp = 0;
+ ctx->prec = 2*prec;
+ ctx2->prec = prec;
+ ctx2->rounding = ROUND_HALF_EVEN;
+
+
+ ans = _NEW_decimalobj(1, 0, exp_from_i(0));
+ if (!ans) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ Py_DECREF(t);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+ ans->limbs[0] = 1;
+ while (1) {
+ decimalobject *tmp;
+ tmp = _do_decimal_add(ans, t, ctx);
+ Py_DECREF(ans);
+ if (!tmp) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ Py_DECREF(t);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+ ans = tmp;
+
+ tmp = _do_decimal_multiply(t, rescaled, ctx2);
+ Py_DECREF(t);
+
+ if (!tmp) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ Py_DECREF(ans);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+
+ t = tmp;
+
+ tmp = _do_decimal__divide(t, d, 0, ctx2);
+ Py_DECREF(t);
+
+ if (!tmp) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ Py_DECREF(ans);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+
+ /* we end, when |x^i/i!| cannot affect result
+ * we now, that when S(i) = x^0/0! + x^1/1! ... x^i/i!, then
+ * S(i) - x^i/i! <= e^x <= S(i) + x^i/i! */
+ t = tmp;
+ if (ans->ob_size >= prec) {
+ exp_t tmp_exp = exp_add_i(ADJUSTED(t), prec + 1);
+ if (exp_ge(ADJUSTED(ans), tmp_exp))
+ break;
+ }
+
+ ctx2->prec = 16;
+ ctx2->Emax = exp_from_i(1000);
+ tmp = _decimal_increment(d, 1, ctx2);
+ Py_DECREF(d);
+ ctx2->prec = prec;
+ ctx2->Emax = ctx->Emax;
+ if (!tmp) {
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ Py_DECREF(d);
+ Py_DECREF(ans);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ }
+
+ d = tmp;
+ }
+
+ Py_DECREF(t);
+ Py_DECREF(d);
+ Py_DECREF(rescaled);
+ Py_DECREF(ctx2);
+ }
+
+
+ ctx->prec = prec + 2;
+ if (h) {
+ long pow;
+ decimalobject *y;
+ decimalobject *tmp;
+ int i;
+ pow = 1;
+ for (i = 0; i < h; i++)
+ pow *= 10;
+
+ y = decimal_from_long(self->ob_type, pow);
+
+ if (!y) {
+ Py_DECREF(ans);
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+
+ Py_INCREF(Py_None);
+ tmp = _do_decimal_power(ans, y, Py_None, ctx);
+ Py_DECREF(Py_None);
+ Py_DECREF(y);
+ Py_DECREF(ans);
+ if (!tmp) {
+ ctx->prec = precision;
+ ctx->clamp = clamped;
+ return NULL;
+ }
+ ans = tmp;
+ }
+
+ ctx->clamp = clamped;
+ ctx->prec = precision;
+
+ /* we'll add 1 at the end, to make proper rounding */
+ if (!ISSPECIAL(ans) && decimal_nonzero(ans)) {
+ long i;
+ decimalobject *tmp = _NEW_decimalobj(ans->ob_size + LOG, ans->sign, exp_sub_i(ans->exp, LOG));
+ if (!tmp) {
+ Py_DECREF(ans);
+ return NULL;
+ }
+
+
+ for (i = 0; i<ans->limb_count;i++) {
+ tmp->limbs[i+1] = ans->limbs[i];
+ }
+ tmp->limbs[0] = 1;
+ Py_DECREF(ans);
+ ans = tmp;
+
+ }
+ {
+ decimalobject *fixed = _decimal_fix(ans, ctx);
+ Py_DECREF(ans);
+ ans = fixed;
+ }
+ ctx->rounding = rounding;
+ return ans;
+}
+
+DECIMAL_UNARY_FUNC(exponent);
+
+static PyObject *
+_do_decimal_exp(decimalobject *self, contextobject *ctx) {
+ if (!ctx)
+ ctx = getcontext();
+ if (!ctx)
+ return NULL;
+
+ if (ctx->prec > 999999 ||
+ exp_g_i(ctx->Emax, 999999) ||
+ exp_l_i(ctx->Emin,-999999))
+ handle_InvalidContext(self->ob_type, ctx, NULL);
+ return _do_decimal_exponent(self, ctx);
+}
+
static PyMethodDef decimal_methods[] = {
+ {"exp", (PyCFunction)decimal_exponent,
+ METH_VARARGS | METH_KEYWORDS,
+ PyDoc_STR("TODO")},
+ {"_exponent", (PyCFunction)decimal_exponent,
+ METH_VARARGS | METH_KEYWORDS,
+ PyDoc_STR("TODO")},
{"_divide", (PyCFunction)decimal__divide,
METH_VARARGS | METH_KEYWORDS,
PyDoc_STR("TODO")},
@@ -4399,7 +4709,12 @@
if (!res) return NULL;
res->digits[0] = 0;
res->limbs[0] = 0;
- return res;
+ ret = _decimal_fix(res, ctx);
+ Py_DECREF(res);
+ if (!ret)
+ return NULL;
+
+ return ret;
}
if (!decimal_nonzero(self)) {
oexp = exp_sub_i(other->exp, ctx->prec + 1);
@@ -6398,6 +6713,7 @@
CONTEXT_BINARY_FUNC(remainder_near, remainder_near)
CONTEXT_UNARY_FUNC(sqrt, sqrt)
CONTEXT_UNARY_FUNC(to_eng_string, to_eng_string)
+CONTEXT_UNARY_FUNC(exp, exp)
/* Unfortunately, the following methods are non-standard and can't
@@ -6920,6 +7236,8 @@
METH_VARARGS},
{"sqrt", (PyCFunction)context_sqrt,
METH_O},
+ {"exp", (PyCFunction)context_exp,
+ METH_O},
{"subtract", (PyCFunction)context_subtract,
METH_VARARGS},
{"to_eng_string", (PyCFunction)context_to_eng_string,
More information about the Python-checkins
mailing list