[Python-checkins] r47095 - sandbox/trunk/decimal-c/_decimal.c
mateusz.rukowicz
python-checkins at python.org
Mon Jun 26 04:19:07 CEST 2006
Author: mateusz.rukowicz
Date: Mon Jun 26 04:19:06 2006
New Revision: 47095
Modified:
sandbox/trunk/decimal-c/_decimal.c
Log:
Square rooting works, some minor bugs in rounding functions fixed, when limb_count != (ob_size + LOG -1)/LOG.
Modified: sandbox/trunk/decimal-c/_decimal.c
==============================================================================
--- sandbox/trunk/decimal-c/_decimal.c (original)
+++ sandbox/trunk/decimal-c/_decimal.c Mon Jun 26 04:19:06 2006
@@ -689,6 +689,11 @@
static PyObject *context_regard_flags(contextobject *, PyObject*);
static decimalobject *_do_decimal_absolute(decimalobject *, contextobject *, int);
static PyObject *decimal_int(decimalobject*);
+static decimalobject *_do_decimal_add(decimalobject *, decimalobject *, contextobject *);
+static PyObject *_do_decimal__divide(decimalobject *, decimalobject *, int, contextobject *);
+static decimalobject *_decimal_fromliteral(PyTypeObject *, char *str, long, contextobject *);
+static decimalobject *_do_decimal_multiply(decimalobject *, decimalobject *, contextobject *);
+static decimalobject *_do_decimal_subtract(decimalobject *, decimalobject *, contextobject *);
/* Exception handlers *********************************************************/
@@ -1014,6 +1019,7 @@
}
if(_limb_get_digit(new->limbs, new->ob_size, 0) == 0)
new->ob_size --;
+ new->limb_count = (new->ob_size + LOG - 1)/ LOG;
return new;
}
@@ -1052,6 +1058,7 @@
if (new2->ob_size > prec) {
_limb_cut_one_digit(new2->limbs,new2->ob_size);
new2->ob_size--;
+ new2->limb_count = (new2->ob_size + LOG -1)/LOG;
new2->exp++;
}
return new2;
@@ -1077,6 +1084,7 @@
if (new->ob_size > prec) {
_limb_cut_one_digit(new->limbs,new->ob_size);
new->ob_size--;
+ new->limb_count = (new->ob_size + LOG - 1)/LOG;
new->exp++;
}
return new;
@@ -1277,6 +1285,7 @@
_limb_cut_one_digit(new->limbs, new->ob_size); /* VERY SLOW */
new->ob_size --;
}
+ new->limb_count = (new->ob_size + LOG - 1)/LOG;
new->exp += expdiff;
if (handle_Rounded(ctx, NULL) != 0) {
Py_DECREF(new);
@@ -2070,10 +2079,356 @@
}
+/* this code is rewritten from python version, this might change
+ * in future */
+
static decimalobject *
_do_decimal_sqrt(decimalobject *self, contextobject *ctx)
{
- Py_RETURN_NONE;
+ decimalobject *ret = 0;
+ decimalobject *ans = 0;
+ decimalobject *tmp = 0, *tmp2 = 0, *tmp3 = 0;
+ contextobject *ctx2 = 0;
+ decimalobject *half = 0;
+ PyObject *flags = 0;
+ long expadd;
+ long firstprec;
+ long i;
+ long Emax;
+ long Emin;
+ long maxp;
+ long rounding;
+ long prevexp;
+
+ if (ISSPECIAL(self)) {
+ decimalobject *nan;
+ if (_check_nans(self, NULL, ctx, &nan))
+ return nan;
+
+ if (ISINF(self) && (self->sign&1) == 0)
+ return _decimal_get_copy(self);
+ }
+
+ if (!decimal_nonzero(self)) {
+ ret = _NEW_decimalobj(1, 0, 0);
+ if (!ret)
+ return NULL;
+
+ ret->limbs[0] = 0;
+ ret->exp = self->exp / 2;
+ if (self->exp < 0 && self->exp%2)
+ ret->exp --;
+ ret->sign = self->sign & 1;
+
+ return ret;
+ }
+
+ if (self->sign&1) {
+ return handle_InvalidOperation(self->ob_type, ctx, "sqrt(-x), x > 0", NULL);
+ }
+
+ tmp = _NEW_decimalobj(self->ob_size + 1, self->sign, self->exp);
+
+ if (!tmp)
+ return NULL;
+
+ expadd = tmp->exp/2;
+ if (tmp->exp < 0 && tmp->exp % 2)
+ expadd --;
+
+ if (tmp->exp &1) {
+ _limb_first_n_digits(self->limbs, self->ob_size, 0, tmp->limbs, tmp->ob_size);
+ }
+ else {
+ tmp->ob_size --;
+ tmp->limb_count = self->limb_count;
+ for(i = 0; i < tmp->limb_count ; i++)
+ tmp->limbs[i] = self->limbs[i];
+ }
+
+ tmp->exp = 0;
+
+ ctx2 = context_copy(ctx);
+
+ if (!ctx2) {
+ Py_DECREF(tmp);
+ return NULL;
+ }
+
+ flags = context_ignore_all_flags(ctx2);
+
+ if (!flags) {
+ Py_DECREF(tmp);
+ Py_DECREF(ctx2);
+ return NULL;
+ }
+
+ ans = _NEW_decimalobj(3, 0, 0);
+
+ if (!ans) {
+ Py_DECREF(tmp);
+ Py_DECREF(ctx2);
+ Py_DECREF(flags);
+ }
+
+ tmp2 = _NEW_decimalobj(3, 0, 0);
+
+ if (!tmp2) {
+ Py_DECREF(tmp);
+ Py_DECREF(ctx2);
+ Py_DECREF(flags);
+ Py_DECREF(ans);
+ }
+
+ if ((ADJUSTED(tmp) & 1) == 0) {
+ ans->limbs[0] = 819;
+ ans->exp = ADJUSTED(tmp) - 2;
+ tmp2->limbs[0] = 259;
+ tmp2->exp = -2;
+ }
+ else {
+ ans->limbs[0] = 259;
+ ans->exp = tmp->exp + tmp->ob_size - 3;
+ tmp2->limbs[0] = 819;
+ tmp2->exp = -3;
+ }
+
+ firstprec = ctx2->prec;
+ ctx2->prec = 3;
+ /* ans += tmp * tmp2 */
+
+ tmp3 = _do_decimal_multiply(tmp, tmp2, ctx2);
+
+ if (!tmp3)
+ goto err;
+
+ Py_DECREF(tmp2);
+ tmp2 = _do_decimal_add(ans, tmp3, ctx2);
+ Py_DECREF(tmp3);
+ tmp3 = 0;
+
+ if (!tmp2)
+ goto err;
+
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+ ans->exp -= 1 + ADJUSTED(tmp)/2;
+ if (1 + ADJUSTED(tmp) < 0 && (1 + ADJUSTED(tmp)) % 2)
+ ans->exp --;
+
+ Emax = ctx2->Emax;
+ Emin = ctx2->Emin;
+
+ ctx2->Emax = PyDecimal_DefaultContext->Emax;
+ ctx2->Emin = PyDecimal_DefaultContext->Emin;
+
+ half = _decimal_fromliteral(self->ob_type, "0.5", 3, ctx2);
+
+ if (!half)
+ goto err;
+
+ maxp = firstprec + 2;
+ rounding = ctx2->rounding;
+ ctx2->rounding = ROUND_HALF_EVEN;
+
+ while (1) {
+ ctx2->prec = 2 * ctx2->prec - 2 < maxp ? 2 * ctx2->prec - 2: maxp;
+ /* ans = half * (ans + tmp/ans) */
+ tmp2 = _do_decimal__divide(tmp, ans, 0, ctx2);
+ if (!tmp2)
+ goto err;
+ /* ans = half * (ans + tmp2) */
+ tmp3 = _do_decimal_add(ans, tmp2, ctx2);
+ if (!tmp3)
+ goto err;
+
+ /* ans = half * tmp3 */
+ Py_DECREF(ans);
+ ans = _do_decimal_multiply(half, tmp3, ctx2);
+
+ if (!ans)
+ goto err;
+
+ Py_DECREF(tmp2);
+ Py_DECREF(tmp3);
+ tmp2 = 0;
+ tmp3 = 0;
+
+ if (ctx2->prec == maxp)
+ break;
+ }
+
+ ctx2->prec = firstprec;
+ prevexp = ADJUSTED(ans);
+ tmp2 = _decimal_round(ans, -1, ctx2, -1);
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = _NEW_decimalobj(tmp2->ob_size + 1, tmp2->sign, tmp2->exp);
+ if (!ans)
+ goto err;
+
+ ctx2->prec = firstprec + 1;
+ if (prevexp != ADJUSTED(tmp2)) {
+ _limb_first_n_digits(tmp2->limbs, tmp2->ob_size, 0, ans->limbs, ans->ob_size);
+ ans->exp --;
+ }
+ else {
+ ans->limb_count = tmp2->limb_count;
+ ans->ob_size = tmp2->ob_size;
+ for (i = 0; i < ans->limb_count; i++)
+ ans->limbs[i] = tmp2->limbs[i];
+ }
+
+ Py_DECREF(tmp2);
+ tmp2 = 0;
+
+ {
+ int cmp;
+ decimalobject *lower;
+ half->exp = ans->exp - 1;
+ half->limbs[0] = 5;
+ lower = _do_decimal_subtract(ans, half, ctx2);
+ if (!lower)
+ goto err;
+
+ ctx2->rounding = ROUND_UP;
+ tmp2 = _do_decimal_multiply(lower, lower, ctx2);
+ Py_DECREF(lower);
+ lower = 0;
+ if (!tmp2) {
+ goto err;
+ }
+
+ cmp = _do_real_decimal_compare(tmp2, tmp, ctx2);
+ if (PyErr_Occurred()) {
+ goto err;
+ }
+
+ Py_DECREF(tmp2);
+ tmp2 = 0;
+
+ if (cmp == 1) {
+ half->exp = ans->exp;
+ half->limbs[0] = 1;
+ tmp2 = _do_decimal_subtract(ans, half, ctx2);
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+ }
+ else {
+ decimalobject *upper;
+ half->exp = ans->exp-1;
+ half->limbs[0] = 5;
+ upper = _do_decimal_add(ans, half, ctx2);
+ if (!upper)
+ goto err;
+ ctx2->rounding = ROUND_DOWN;
+
+ tmp2 = _do_decimal_multiply(upper, upper, ctx2);
+
+ Py_DECREF(upper);
+ upper = 0;
+
+ cmp = _do_real_decimal_compare(tmp2, tmp, ctx2);
+ if (PyErr_Occurred())
+ goto err;
+
+ Py_DECREF(tmp2);
+ tmp2 = 0;
+
+ if (cmp == -1) {
+ half->exp = ans->exp;
+ half->limbs[0] = 1;
+ tmp2 = _do_decimal_add(ans, half, ctx2);
+
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+ }
+ }
+ }
+
+ ans->exp += expadd;
+ ctx2->rounding = rounding;
+
+ tmp2 = _decimal_fix(ans, ctx2);
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+
+ rounding = ctx2->rounding_dec;
+ ctx2->rounding_dec = NEVER_ROUND;
+
+ {
+ int cmp;
+ tmp2 = _do_decimal_multiply(ans, ans, ctx2);
+ if (!tmp2)
+ goto err;
+
+ cmp = _do_real_decimal_compare(tmp2, self, ctx2);
+
+ if (PyErr_Occurred())
+ goto err;
+
+ Py_DECREF(tmp2);
+ tmp2 = 0;
+
+ if (cmp != 0) {
+ if (handle_Rounded(ctx, NULL))
+ goto err;
+
+ if (handle_Inexact(ctx, NULL))
+ goto err;
+ }
+
+ else {
+ long exp = self->exp / 2;
+ if (self->exp < 0 && self->exp % 2)
+ exp --;
+ ctx2->prec += ans->exp - exp;
+ tmp2 = _decimal_rescale(ans, exp, ctx2, -1, 1);
+
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+ ctx2->prec = firstprec;
+ }
+ }
+
+ tmp2 = _decimal_fix(ans, ctx);
+ if (!tmp2)
+ goto err;
+ Py_DECREF(ans);
+ ans = tmp2;
+ tmp2 = 0;
+
+ Py_DECREF(flags);
+ Py_DECREF(ctx2);
+ Py_DECREF(half);
+ Py_DECREF(tmp);
+
+ return ans;
+
+err:
+ Py_XDECREF(tmp);
+ Py_XDECREF(tmp2);
+ Py_XDECREF(tmp3);
+ Py_XDECREF(ans);
+ Py_XDECREF(flags);
+ Py_XDECREF(ctx2);
+ Py_XDECREF(ret);
+ Py_XDECREF(half);
+ return NULL;
}
DECIMAL_UNARY_FUNC(sqrt)
@@ -3510,6 +3865,7 @@
size = _limb_multiply(self->limbs, self->ob_size, other->limbs, other->ob_size, ans->limbs);
ans->ob_size = size;
+ ans->limb_count = (size + LOG -1 )/LOG;
done:
More information about the Python-checkins
mailing list