/* * 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 #include #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); }