diff options
Diffstat (limited to 'ic-reals-6.3/base/realLib.c')
| -rw-r--r-- | ic-reals-6.3/base/realLib.c | 486 |
1 files changed, 486 insertions, 0 deletions
diff --git a/ic-reals-6.3/base/realLib.c b/ic-reals-6.3/base/realLib.c new file mode 100644 index 0000000..15da22a --- /dev/null +++ b/ic-reals-6.3/base/realLib.c @@ -0,0 +1,486 @@ +/* + * 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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "real.h" +#include "real-impl.h" + +void initEmitDigit(void); +void initEmitSign(void); +void initEpsDel(); +void initTmp(); +void initStack(); +void initForceMethodLookupTable(); +void initStrategy(); + +int defaultForceCount = DEFAULT_FORCE_COUNT; +int forceDecUpperBound = FORCE_DEC_UPPER_BOUND; +int stackSize = STACK_SIZE * 1024; + +/* + * Initialization of the real library + */ +void +initRealBase() +{ + char *p; + + if ((p = getenv("ICR_DEFAULT_FORCE_COUNT")) != NULL) + defaultForceCount = atoi(p); + + if ((p = getenv("ICR_FORCE_DEC_UPPER_BOUND")) != NULL) + forceDecUpperBound = atoi(p); + + if ((p = getenv("ICR_STACK_SIZE")) != NULL) + stackSize = atoi(p) * 1024; + +#ifdef DAVINCI + initDaVinci(); +#endif + + initStack(); + initForceMethodLookupTable(); + initEmitDigit(); + initEmitSign(); + initTmp(); + initEpsDel(); + initStrategy(); +} + +/* + * For a real x, this returns a new real which denotes the same value as + * x, but where the new real is a stream (ie prefixed by a SignX or DigsX). + * The only interesting cases are when the expression rooted at x is an + * lft. In this case we create a new real object but also record in x + * a link to the new real to avoid doing it again. + */ +Real +makeStream(Real x) +{ + Real r; + DigsX *digsX; + void force_To_DigsX_From_Vec(); + void force_To_DigsX_From_MatX_Entry(); + void force_To_DigsX_From_TenXY_Entry(); + void force_To_DigsX_From_Cls_Entry(); + void force_To_DigsX_From_Alt_Entry(); + + /* + * When we want to make a stream from a vector or matrix, we + * must make a copy of the lft first since the LFT might be shared + * and emitting a sign or digit changes the lft. For tensors, + * this is unnecessary as, while the tensor might be shared, + * any consumer of the tensor can consume only the stream + * and changes to the original tensor are acceptable. + */ + switch (x->gen.tag.type) { + case VECTOR : + /* + * Check that we haven't already made a stream for this lft + */ + if (x->vec.strm == NULL) { + + /* create a copy of the lft */ + r = vector_Z(x->vec.vec[0], x->vec.vec[1]); + + /* to make a signed lft into a stream, we prefix it by a SignX */ + if (r->gen.tag.isSigned) + x->vec.strm = (Real) allocSignX(r, SIGN_UNKN); + + /* for an unsigned lft, we prefix it by a DigsX */ + else { + digsX = allocDigsX(); + digsX->force = force_To_DigsX_From_Vec; + digsX->x = r; +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, r); + endGraphUpdate(); +#endif + x->vec.strm = (Real) digsX; + } + } + + /* + * If daVinci is running then we draw a double line from the original + * real and its stream version to indicate that there are two + * expressions in the heap which point denote the same real value. + */ +#ifdef DAVINCI + beginGraphUpdate(); + drawEqEdge(x, x->vec.strm); + endGraphUpdate(); +#endif + return x->vec.strm; + break; + + /* this case is exactly the same as the above */ + case MATX : + if (x->matX.strm == NULL) { + r = matrix_Z(x->matX.x, x->matX.mat[0][0], x->matX.mat[0][1], + x->matX.mat[1][0], x->matX.mat[1][1]); + if (r->gen.tag.isSigned) + x->matX.strm = (Real) allocSignX(r, SIGN_UNKN); + else { + digsX = allocDigsX(); + digsX->force = force_To_DigsX_From_MatX_Entry; + digsX->x = r; +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, r); + endGraphUpdate(); +#endif + x->matX.strm = (Real) digsX; + } + } +#ifdef DAVINCI + beginGraphUpdate(); + drawEqEdge(x, x->matX.strm); + endGraphUpdate(); +#endif + return x->matX.strm; + break; + + /* same as for vectors */ + case TENXY : + if (x->tenXY.strm == NULL) { + r = tensor_Z(x->tenXY.x, x->tenXY.y, + x->tenXY.ten[0][0], x->tenXY.ten[0][1], + x->tenXY.ten[1][0], x->tenXY.ten[1][1], + x->tenXY.ten[2][0], x->tenXY.ten[2][1], + x->tenXY.ten[3][0], x->tenXY.ten[3][1]); + if (r->gen.tag.isSigned) + x->tenXY.strm = (Real) allocSignX(r, SIGN_UNKN); + else { + digsX = allocDigsX(); + digsX->force = force_To_DigsX_From_TenXY_Entry; + digsX->x = r; +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, r); + endGraphUpdate(); +#endif + x->tenXY.strm = (Real) digsX; + } +#ifdef DAVINCI + beginGraphUpdate(); + drawEqEdge(x, x->tenXY.strm); + endGraphUpdate(); +#endif + } + return x->tenXY.strm; + break; + + case ALT : + if (x->alt.tag.isSigned) + r = (Real) allocSignX(x, SIGN_UNKN); + else { + r = (Real) allocDigsX(); + r->digsX.force = force_To_DigsX_From_Alt_Entry; + r->digsX.x = x; +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(r, x); + endGraphUpdate(); +#endif + } + return r; + break; + + case CLOSURE : + if (x->cls.tag.isSigned) + r = (Real) allocSignX(x, SIGN_UNKN); + else { + r = (Real) allocDigsX(); + r->digsX.force = force_To_DigsX_From_Cls_Entry; + r->digsX.x = x; +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(r, x); + endGraphUpdate(); +#endif + } + return r; + break; + + case DIGSX : + case SIGNX : + return x; + + default : + Error(FATAL, E_INT, "makeStream", + "trying to make a stream from a non-real"); + return NULL; + } +} + +/* + * This creates a closure (typically a function which unfolds a tensor + * chain). + */ +Cls * +allocCls(void (*force)(), void *userData) +{ + Cls *cls; + + if ((cls = (Cls *) malloc (sizeof(Cls))) == NULL) + Error(FATAL, E_INT, "allocCls", "malloc failed"); + +#ifdef DAVINCI + newNodeId(cls); +#else +#ifdef TRACE + newNodeId(cls); +#endif +#endif + + cls->tag.type = CLOSURE; + cls->tag.dumped = FALSE; + cls->tag.isSigned = FALSE; + cls->force = force; + cls->userData = userData; + cls->redirect = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + newNode(cls, CLOSURE); + endGraphUpdate(); +#endif + + return cls; +} + +/* + * The next family of functions are simple arithmetic operations. + */ +Real +add_R_R(Real x, Real y) +{ + return tensor_Int(x, y, 0, 0, 1, 0, 1, 0, 0, 1); +} + +Real +add_R_Int(Real x, int y) +{ + return matrix_Int(x, 1, 0, y, 1); +} + +Real +add_R_QInt(Real x, int a, int b) +{ + return matrix_Int(x, b, 0, a, b); +} + +Real +sub_R_R(Real x, Real y) +{ + return tensor_Int(x, y, 0, 0, 1, 0, -1, 0, 0, 1); +} + +Real +sub_R_Int(Real x, int y) +{ + return matrix_Int(x, 1, 0, -y, 1); +} + +Real +sub_R_QInt(Real x, int a, int b) +{ + return matrix_Int(x, b, 0, -a, b); +} + +Real +sub_Int_R(int x, Real y) +{ + return matrix_Int(y, -1, 0, x, 1); +} + +Real +sub_QInt_R(int a, int b, Real y) +{ + return matrix_Int(y, -b, 0, a, b); +} + +Real +mul_R_R(Real x, Real y) +{ + return tensor_Int(x, y, 1, 0, 0, 0, 0, 0, 0, 1); +} + +Real +mul_R_Int(Real x, int y) +{ + return matrix_Int(x, y, 0, 0, 1); +} + +Real +mul_R_QInt(Real x, int a, int b) +{ + return matrix_Int(x, a, 0, 0, b); +} + +Real +div_R_R(Real x, Real y) +{ + return tensor_Int(x, y, 0, 0, 1, 0, 0, 1, 0, 0); +} + +Real +div_R_Int(Real x, int y) +{ + return matrix_Int(x, 1, 0, 0, y); +} + +Real +div_Int_R(int x, Real y) +{ + return matrix_Int(y, 0, 1, x, 0); +} + +Real +div_R_QInt(Real x, int a, int b) +{ + return matrix_Int(x, b, 0, 0, a); +} + +Real +div_QInt_R(int a, int b, Real x) +{ + return matrix_Int(x, 0, b, a, 0); +} + + +/* + * This creates a vector (a,b) and then prepends it by the characteristic + * pair (c,n) and finally sets the sign to sign. This is used only for + * debugging .. to set up a real argument for a function. + */ +Real +makeRealSignCNQInt(Sign sign, char *c, int n, int a, int b) +{ + DigsX *digsX; + Real y; + void force_To_DigsX_From_Vec(); + void force_To_DigsX_From_DigsX_Entry(); + + y = (Real) vector_Int(a, b); + + if (n > 0) { + digsX = allocDigsX(); + digsX->x = y; + digsX->force = force_To_DigsX_From_Vec; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, y); + endGraphUpdate(); +#endif + + y = (Real) digsX; + + digsX = allocDigsX(); + digsX->x = y; + digsX->force = force_To_DigsX_From_DigsX_Entry; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, y); + endGraphUpdate(); +#endif + + digsX->count = n; +#ifdef PACK_DIGITS + if (n <= DIGITS_PER_WORD) + digsX->word.small = atoi(c); + else + mpz_init_set_str(digsX->word.big, c, 10); +#else + mpz_set_str(digsX->word.big, c, 10); +#endif + return (Real) allocSignX((Real) digsX, sign); + } + else { + switch (sign) { + case SZERO : + case SPOS : + case SNEG : + case SINF : + digsX = allocDigsX(); + digsX->x = y; + digsX->force = force_To_DigsX_From_Vec; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, y); + endGraphUpdate(); +#endif + + digsX->count = 0; + return (Real) allocSignX((Real) digsX, sign); + + case SIGN_UNKN : + return (Real) allocSignX(y, sign); + + default : + Error(FATAL, E_INT, "makeRealSignCNQInt", "invalid sign"); + return NULL; + } + } +} + +/* + * These are the new preferred names for the basic lft functions. + * The other names will disappear in a future release. + */ +Real +real_QInt(int a, int b) +{ + return vector_Int(a, b); +} + +Real +real_QZ(mpz_t a, mpz_t b) +{ + return vector_Z(a, b); +} + +/* + * These are the same as matrix_ functions but take their arguments + * in a different order (row order rather than column order). + */ +Real +lft_R_Int(Real x, int a, int b, int c, int d) +{ + return matrix_Int(x, a, c, b, d); +} + +Real +lft_R_Z(Real x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) +{ + return matrix_Z(x, a, c, b, d); +} + +/* + * These are the same as tensor_ functions but take their arguments + * in a different order (row order rather than column order). + */ +Real +lft_R_R_Int(Real x, Real y, int a, int b, int c, int d, + int e, int f, int g, int h) +{ + return tensor_Int(x, y, a, e, b, f, c, g, d, h); +} + +Real +lft_R_R_Z(Real x, Real y, mpz_t a, mpz_t b, mpz_t c, mpz_t d, + mpz_t e, mpz_t f, mpz_t g, mpz_t h) +{ + return tensor_Z(x, y, a, e, b, f, c, g, d, h); +} |