From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- ic-reals-6.3/math-lib/log_R.c~ | 227 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 ic-reals-6.3/math-lib/log_R.c~ (limited to 'ic-reals-6.3/math-lib/log_R.c~') diff --git a/ic-reals-6.3/math-lib/log_R.c~ b/ic-reals-6.3/math-lib/log_R.c~ new file mode 100644 index 0000000..bc66045 --- /dev/null +++ b/ic-reals-6.3/math-lib/log_R.c~ @@ -0,0 +1,227 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +void log2Cont(); +void logInside(); +void logLow(); +void logHigh(); + +/* + * This diverges when the argument is zero. + */ +Real +log_R(Real x) +{ + Bool xGtEq0, xLtEq2, xGtEqOneHalf, xLtEq2_and_GtEqOneHalf; + Bool xLtEqOverOneHalf_and_GtEq0; + Real in, low, high, ltZero; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(logInside, "logInside", 2); + registerForceFunc(logLow, "logLow", 2); + registerForceFunc(logHigh, "logHigh", 2); + registerForceFunc(log2Cont, "log2Cont", 3); + doneInit++; + } + + xGtEq0 = gtEq_R_0(x); + xLtEq2 = ltEq_R_QInt(x, 2, 1); + xGtEqOneHalf = gtEq_R_QInt(x, 1, 2); + xLtEq2_and_GtEqOneHalf = and_B_B(xLtEq2, xGtEqOneHalf); + xLtEqOverOneHalf_and_GtEq0 = and_B_B(ltEq_R_QInt(x, 1001, 2000), xGtEq0); + + in = (Real) allocCls(logInside, (void *) x); + in->cls.tag.isSigned = TRUE; + low = (Real) allocCls(logLow, (void *) x); + low->cls.tag.isSigned = TRUE; + high = (Real) allocCls(logHigh, (void *) x); + high->cls.tag.isSigned = TRUE; + ltZero = realError("(log_R x) and x < 0"); + + return realIf(4, + xLtEq2_and_GtEqOneHalf, in, + gtEq_R_QInt(x, 1999, 1000), high, + xLtEqOverOneHalf_and_GtEq0, low, + not_B(xGtEq0), ltZero); +} + +static TenXY *nextTensor(Real, Real, int); + +void +logInside() +{ + Cls *cls, *newCls; + ClsData *data; + Real x; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "logInside", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 0, 1, 1, -1, 1, -1, 0); + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +void +logHigh() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = div_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = add_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +void +logLow() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = mul_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = sub_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +/* + * This closure works a little different from most. Most simply create + * an lft and set the redirect field of the closure to point to the new + * heap object. In this case, we can do a bit better. This continuation + * represents (log 2). One option would be to follow the usual strategy + * and to arrange for (log 2) to be shared whenever possible. However, + * we know that this is used internally to log_R only and only in + * those circumstances when the argument falls outside the working range + * of the log tensor chain. Moreover, we know that the consumer of this + * real is either an addition or substraction tensor. Whatever it is, + * it is certain to be a tensor. The other argument of the tensor will + * reduce at most to a matrix and hence the tensor itself will never + * reduce. + * + * This property allows us to adopt a different strategy. Rather than create + * the next matrix in the (log 2) chain, we put the matrix directly into + * the consuming tensor. This avoids creating garbage in the stack and + * avoids a separate reduction step. + * + * There is a slighlt better scheme than this. Rather than plunk in + * each successive + * matrix into the tensor (each providing 4 digits), it would be better + * to accumulate a sequence of matrices (up to the capacity of a small matrix) + * and then plunk it into the tensor. This would avoid using bignum stuff until + * it becomes necessary. Perhaps I'll do this later. + * + * For this it is useful to bear in mind that when a sequence of matrices + * starts from n = 0, then the largest entry is always d. When the sequence + * starts with n > 0, then the largest entry is c. To decide if there + * is going to be an overflow, it suffices to check if c (for example) + * will overflow. + */ +static void nextMatrix(Tensor, int); + +void +log2Cont() +{ + TenXY *tenXY; + Cls *cls; + int digitsNeeded; + int n; + + tenXY = (TenXY *) POP; + digitsNeeded = POP; + cls = (Cls *) tenXY->y; + + n = (int) cls->userData; + + while (digitsNeeded > 0) { + nextMatrix(tenXY->ten, n); + if (n == 0) + digitsNeeded -= 2; + else + digitsNeeded -= 4; + n += 1; + } + cls->userData = (void *) n; +} + +/* + * We set a rather arbitrary (but large) limit on the value of n + * for log 2. On a 32 bit machine it is 536,870,911. + * Larger than this and the matrix entries are no longer small. + */ +static void +nextMatrix(Tensor ten, int n) +{ + SmallMatrix smallMat; + + if (n == 0) { + smallMat[0][0] = 1; + smallMat[0][1] = 1; + smallMat[1][0] = 1; + smallMat[1][1] = 2; + } + else { + if (n <= (MAXINT - 2) / 4) { + smallMat[0][0] = 3 * n + 1; + smallMat[0][1] = 2 * n + 1; + smallMat[1][0] = 4 * n + 2; + smallMat[1][1] = 3 * n + 2; + } + else + Error(FATAL, E_INT, "nextMatrix (log2)", + "n > %d", (MAXINT - 2) / 4); + } + multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat); + multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat); + normalizeTensor(ten); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, n, 0, 2*n+1, n+1, n+1, 2*n+1, 0, n); +} -- cgit v1.2.3