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/base/emitSign.c | 618 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 618 insertions(+) create mode 100644 ic-reals-6.3/base/emitSign.c (limited to 'ic-reals-6.3/base/emitSign.c') diff --git a/ic-reals-6.3/base/emitSign.c b/ic-reals-6.3/base/emitSign.c new file mode 100644 index 0000000..d25aa03 --- /dev/null +++ b/ic-reals-6.3/base/emitSign.c @@ -0,0 +1,618 @@ +/* + * 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" + +/* + * This file deals with emitting signs from LFTs. The algorithm is similar + * to digit emission. See ``Notes on sign and digit emission'' for details. + */ + +/* + * This version of emitSign favours SZERO over all other signs. This + * is because most of the analytic functions prefer SZERO. + */ +typedef struct TE { + Sign sign; + Comparison compare; + void (*residual)(Vector); /* the function to compute the residual */ + struct TE *next; /* computes next entry in the sequence */ +} TblEntry; + +static bool mergeTwoLists(TblEntry **, TblEntry **, Comparison *); +static bool mergeFourLists(TblEntry **, TblEntry **, TblEntry **, TblEntry **); +static TblEntry *setList(Vector); + +static mpz_t tmp_z; + +/* + * Functions prefixed with 'r' compute the residual vectors. This is done + * "in-place". By this we mean that once we have decided that a sign 's' + * can be emitted from a vector 'v', then we overwrite 'v' with residual + * after 's' is emitted. + */ + +static void rNoOp(Vector); + +#define rgg rNoOp +#define regSPOS rNoOp +#define rgeSPOS rNoOp + +static void regSNEG(Vector); +static void regSZERO(Vector); +static void relSPOS(Vector); +static void relSNEG(Vector); +static void relSZERO(Vector); +static void rgeSNEG(Vector); +static void rgeSINF(Vector); +static void rleSPOS(Vector); +static void rleSNEG(Vector); +static void rleSINF(Vector); +static void rggeSZERO(Vector); +static void rggeSINF(Vector); +static void rgggSINF(Vector); +static void rgglSZERO(Vector); +static void rll(Vector); +static void rlleSZERO(Vector); +static void rlleSINF(Vector); +static void rllgSZERO(Vector); +static void rlllSINF(Vector); +static void rgl(Vector); +static void rgleSZERO(Vector); +static void rgleSINF(Vector); +static void rglgSZERO(Vector); +static void rgllSINF(Vector); +static void rlg(Vector); +static void rlgeSZERO(Vector); +static void rlgeSINF(Vector); +static void rlggSINF(Vector); +static void rlglSZERO(Vector); + +static TblEntry tee4 = {SNEG, EQ, rNoOp, NULL}; +static TblEntry tee3 = {SPOS, EQ, rNoOp, &tee4}; +static TblEntry tee2 = {SINF, EQ, rNoOp, &tee3}; +static TblEntry tee1 = {SZERO, EQ, rNoOp, &tee2}; + +static TblEntry teg3 = {SNEG, GT, regSNEG, NULL}; +static TblEntry teg2 = {SPOS, GT, regSPOS, &teg3}; +static TblEntry teg1 = {SZERO, GT, regSZERO, &teg2}; + +static TblEntry tel3 = {SNEG, LT, relSNEG, NULL}; +static TblEntry tel2 = {SPOS, LT, relSPOS, &tel3}; +static TblEntry tel1 = {SZERO, LT, relSZERO, &tel2}; + +static TblEntry tge3 = {SNEG, LT, rgeSNEG, NULL}; +static TblEntry tge2 = {SPOS, GT, rgeSPOS, &tge3}; +static TblEntry tge1 = {SINF, GT, rgeSINF, &tge2}; + +static TblEntry tle3 = {SNEG, GT, rleSNEG, NULL}; +static TblEntry tle2 = {SPOS, LT, rleSPOS, &tle3}; +static TblEntry tle1 = {SINF, LT, rleSINF, &tle2}; + +static TblEntry tgge3 = {SPOS, GT, rgg, NULL}; +static TblEntry tgge2 = {SINF, GT, rggeSINF, &tgge3}; +static TblEntry tgge1 = {SZERO, GT, rggeSZERO, &tgge2}; + +static TblEntry tggg2 = {SPOS, GT, rgg, NULL}; +static TblEntry tggg1 = {SINF, GT, rgggSINF, &tggg2}; + +static TblEntry tggl2 = {SPOS, GT, rgg, NULL}; +static TblEntry tggl1 = {SZERO, GT, rgglSZERO, &tggl2}; + +static TblEntry tlle3 = {SPOS, LT, rll, NULL}; +static TblEntry tlle2 = {SINF, LT, rlleSINF, &tlle3}; +static TblEntry tlle1 = {SZERO, LT, rlleSZERO, &tlle2}; + +static TblEntry tllg2 = {SPOS, LT, rll, NULL}; +static TblEntry tllg1 = {SZERO, LT, rllgSZERO, &tllg2}; + +static TblEntry tlll2 = {SPOS, LT, rll, NULL}; +static TblEntry tlll1 = {SINF, LT, rlllSINF, &tlll2}; + +static TblEntry tgle3 = {SNEG, LT, rgl, NULL}; +static TblEntry tgle2 = {SINF, GT, rgleSINF, &tgle3}; +static TblEntry tgle1 = {SZERO, LT, rgleSZERO, &tgle2}; + +static TblEntry tglg2 = {SNEG, LT, rgl, NULL}; +static TblEntry tglg1 = {SZERO, LT, rglgSZERO, &tglg2}; + +static TblEntry tgll2 = {SNEG, LT, rgl, NULL}; +static TblEntry tgll1 = {SINF, GT, rgllSINF, &tgll2}; + +static TblEntry tlge3 = {SNEG, GT, rlg, NULL}; +static TblEntry tlge2 = {SINF, LT, rlgeSINF, &tlge3}; +static TblEntry tlge1 = {SZERO, GT, rlgeSZERO, &tlge2}; + +static TblEntry tlgg2 = {SNEG, GT, rlg, NULL}; +static TblEntry tlgg1 = {SINF, LT, rlggSINF, &tlgg2}; + +static TblEntry tlgl2 = {SNEG, GT, rlg, NULL}; +static TblEntry tlgl1 = {SZERO, GT, rlglSZERO, &tlgl2}; + +void +initEmitSign(void) +{ + mpz_init(tmp_z); +} + +bool +emitSignFromVector(Vector v, Sign *s) +{ + TblEntry *e; + + if ((e = setList(v)) != NULL) { + (e->residual)(v); + *s = e->sign; +#ifdef TRACE + if (TRACE && !validVector(v)) { + printf("emitSignFromVector: bad vector, sign=%d\n", *s); + dumpVector(v); + Error(FATAL, E_INT, "emitSignFromVector", "bad vector"); + } +#endif + return TRUE; + } else + return FALSE; +} + +bool +emitSignFromMatrix(Matrix m, Sign *s) +{ + TblEntry *e0, *e1; + Comparison comp; /* not used */ + + e0 = setList(m[0]); + e1 = setList(m[1]); + + if (mergeTwoLists(&e0, &e1, &comp)) { + (e0->residual)(m[0]); + (e1->residual)(m[1]); + *s = e0->sign; +#ifdef TRACE + if (TRACE && !validMatrix(m)) { + printf("emitSignFromMatrix: bad matrix, sign=%d\n", *s); + dumpMatrix(m); + Error(FATAL, E_INT, "emitSignFromMatrix", "bad matrix"); + } +#endif + return TRUE; + } else + return FALSE; +} + +bool +emitSignFromTensor(Tensor t, Sign *s) +{ + TblEntry *e0, *e1, *e2, *e3; + + e0 = setList(t[0]); + e1 = setList(t[1]); + e2 = setList(t[2]); + e3 = setList(t[3]); + + if (mergeFourLists(&e0, &e1, &e2, &e3)) { + (e0->residual)(t[0]); + (e1->residual)(t[1]); + (e2->residual)(t[2]); + (e3->residual)(t[3]); + *s = e0->sign; +#ifdef TRACE + if (TRACE && !validTensor(t)) { + printf("emitSignFromTensor: bad tensor, sign=%d\n", *s); + dumpTensor(t); + Error(FATAL, E_INT, "emitSignFromTensor", "bad tensor"); + } +#endif + return TRUE; + } else + return FALSE; +} + +static bool +mergeFourLists(TblEntry **e0, TblEntry **e1, TblEntry **e2, TblEntry **e3) +{ + bool flag; + Comparison compLeft, compRight; + + flag = mergeTwoLists(e0, e1, &compLeft) + && mergeTwoLists(e2, e3, &compRight); + + while (flag) { + if ((*e0)->sign == (*e2)->sign) { + if (compLeft == compRight || compLeft == EQ || compRight == EQ) + return TRUE; + else { + *e0 = (*e0)->next; + *e1 = (*e1)->next; + *e2 = (*e2)->next; + *e3 = (*e3)->next; + flag = mergeTwoLists(e0, e1, &compLeft) + && mergeTwoLists(e2, e3, &compRight); + } + } + else { + if ((*e0)->sign < (*e2)->sign) { + *e0 = (*e0)->next; + *e1 = (*e1)->next; + flag = mergeTwoLists(e0, e1, &compLeft); + } + else { + *e2 = (*e2)->next; + *e3 = (*e3)->next; + flag = mergeTwoLists(e2, e3, &compRight); + } + } + } + return FALSE; +} + +static bool +mergeTwoLists(TblEntry **e0, TblEntry **e1, Comparison *comp) +{ + while (*e0 != NULL && *e1 != NULL) { + if ((*e0)->sign == (*e1)->sign) { + if ((*e0)->compare == (*e1)->compare) { + *comp = (*e0)->compare; + return TRUE; + } + else { + if ((*e0)->compare == EQ) { + *comp = (*e1)->compare; /* must be GT or LT */ + return TRUE; + } + else { + if ((*e1)->compare == EQ) { + *comp = (*e0)->compare; /* must be GT or LT */ + return TRUE; + } + else { /* one is LT and the other GT */ + *e0 = (*e0)->next; + *e1 = (*e1)->next; + } + } + } + } + else { + if ((*e0)->sign < (*e1)->sign) + *e0 = (*e0)->next; + else + *e1 = (*e1)->next; + } + } + return FALSE; +} + +/* + * Given a vector, this function returns a pointer to the table entry (a list) + * which classifies the vector in terms of what signs can be emitted for + * that vector. + */ +static TblEntry * +setList(Vector v) +{ + int sign_a, sign_b; + int tmp; + + sign_a = mpz_sgn(v[0]); + sign_b = mpz_sgn(v[1]); + switch (sign_a) { + case 0 : + switch (sign_b) { + case 0 : /* a == 0, b == 0 */ + return &tee1; + case 1 : /* a == 0, b > 0 */ + return &teg1; + case -1 : /* a == 0, b < 0 */ + return &tel1; + default : + Error(FATAL,E_INT,"signFromVector", + "bad value returned by mpz_sgn"); + return NULL; + } + break; + case 1 : + switch (sign_b) { + case 0 : /* a > 0, b == 0 */ + return &tge1; + + case 1 : /* a > 0, b > 0 */ + tmp = mpz_cmp(v[0], v[1]); + switch (MPZ_SIGN(tmp)) { + case 0 : /* a > 0, b > 0, a == b */ + return &tgge1; + case 1 : /* a > b, b > 0, a > b */ + return &tggg1; + case -1 : /* a > b, b > 0, a < b */ + return &tggl1; + } + break; + + case -1 : /* a > 0, b < 0 */ + /* it is faster to negate a twice */ + mpz_neg(v[0], v[0]); + tmp = mpz_cmp(v[0], v[1]); + mpz_neg(v[0], v[0]); + switch (MPZ_SIGN(tmp)) { + case 0 : /* a > 0, b < 0, -a == b */ + return &tgle1; + case 1 : /* a > b, b < 0, -a > b */ + return &tglg1; + case -1 : /* a > b, b < 0, -a < b */ + return &tgll1; + } + break; + + default : + Error(FATAL,E_INT,"signFromVector", + "bad value returned by mpz_sgn"); + return NULL; + } + break; + + case -1 : + switch (sign_b) { + case 0 : /* a < 0, b == 0 */ + return &tle1; + case 1 : /* a < 0, b > 0 */ + /* it is faster to negate a twice */ + mpz_neg(v[0], v[0]); + tmp = mpz_cmp(v[0], v[1]); + mpz_neg(v[0], v[0]); + switch (MPZ_SIGN(tmp)) { + case 0 : /* a < 0, b > 0, -a == b */ + return &tlge1; + case 1 : /* a < b, b > 0, -a > b */ + return &tlgg1; + case -1 : /* a < b, b > 0, -a < b */ + return &tlgl1; + } + break; + + case -1 : /* a < 0, b < 0 */ + tmp = mpz_cmp(v[0], v[1]); + switch (MPZ_SIGN(tmp)) { + case 0 : /* a < 0, b < 0, a == b */ + return &tlle1; + case 1 : /* a < b, b < 0, a > b */ + return &tllg1; + case -1 : /* a < b, b < 0, a < b */ + return &tlll1; + } + break; + + default : + Error(FATAL,E_INT,"signFromVector", + "bad value returned by mpz_sgn"); + return NULL; + } + break; + + default : + Error(FATAL,E_INT,"signFromVector","bad value returned by mpz_sgn"); + return NULL; + } + return NULL; +} + +static void +rNoOp(Vector v) +{ +} + +static void +regSNEG(Vector v) +{ + MPZ_SWAP(v[0], v[1]); +} + +static void +regSZERO(Vector v) +{ + mpz_set(v[0], v[1]); +} + +static void +relSPOS(Vector v) +{ + mpz_neg(v[1], v[1]); +} + +static void +relSNEG(Vector v) +{ + mpz_neg(v[1], v[1]); + MPZ_SWAP(v[0], v[1]); +} + +static void +relSZERO(Vector v) +{ + mpz_neg(v[1], v[1]); + mpz_set(v[0], v[1]); +} + +static void +rgeSNEG(Vector v) +{ + MPZ_SWAP(v[0], v[1]); +} + +static void +rgeSINF(Vector v) +{ + mpz_set(v[1], v[0]); +} + +static void +rleSPOS(Vector v) +{ + mpz_neg(v[0], v[0]); +} + +static void +rleSNEG(Vector v) +{ + mpz_neg(v[0], v[0]); + MPZ_SWAP(v[0], v[1]); +} + +static void +rleSINF(Vector v) +{ + mpz_neg(v[0], v[0]); + mpz_set(v[1], v[0]); +} + +static void +rggeSZERO(Vector v) +{ + mpz_mul_2exp(v[0], v[0], (unsigned long) 1); + mpz_set_ui(v[1], (unsigned long) 0); +} + +static void +rggeSINF(Vector v) +{ + mpz_mul_2exp(v[1], v[0], (unsigned long) 1); + mpz_set_ui(v[0], (unsigned long) 0); +} + +static void +rgggSINF(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_sub(v[0], v[0], v[1]); + mpz_add(v[1], tmp_z, v[1]); +} + +static void +rgglSZERO(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_add(v[0], v[0], v[1]); + mpz_sub(v[1], v[1], tmp_z); +} + +static void +rll(Vector v) +{ + mpz_neg(v[0], v[0]); + mpz_neg(v[1], v[1]); +} + +static void +rlleSZERO(Vector v) +{ + mpz_neg(v[0], v[0]); + mpz_mul_2exp(v[0], v[0], (unsigned long) 1); + mpz_set_ui(v[1], (unsigned long) 0); +} + +static void +rlleSINF(Vector v) +{ + mpz_neg(v[0], v[0]); + mpz_mul_2exp(v[1], v[0], (unsigned long) 1); + mpz_set_ui(v[0], (unsigned long) 0); +} + +static void +rllgSZERO(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_add(v[0], v[0], v[1]); + mpz_neg(v[0], v[0]); + mpz_sub(v[1], tmp_z, v[1]); +} + +static void +rlllSINF(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_sub(v[0], v[1], v[0]); + mpz_add(v[1], tmp_z, v[1]); + mpz_neg(v[1], v[1]); +} + +static void +rgl(Vector v) +{ + mpz_neg(v[1], v[1]); + MPZ_SWAP(v[0], v[1]); +} + +static void +rgleSZERO(Vector v) +{ + mpz_mul_2exp(v[1], v[0], (unsigned long) 1); + mpz_set_ui(v[0], (unsigned long) 0); +} + +static void +rgleSINF(Vector v) +{ + mpz_mul_2exp(v[0], v[0], (unsigned long) 1); + mpz_set_ui(v[1], (unsigned long) 0); +} + +static void +rglgSZERO(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_add(v[0], v[0], v[1]); + mpz_neg(v[0], v[0]); + mpz_sub(v[1], tmp_z, v[1]); /* - (b - a) */ +} + +static void +rgllSINF(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_sub(v[0], v[0], v[1]); + mpz_add(v[1], tmp_z, v[1]); +} + +static void +rlg(Vector v) +{ + mpz_neg(v[0], v[0]); + MPZ_SWAP(v[0], v[1]); +} + +static void +rlgeSZERO(Vector v) +{ + mpz_mul_2exp(v[1], v[1], (unsigned long) 1); + mpz_set_ui(v[0], (unsigned long) 0); +} + +static void +rlgeSINF(Vector v) +{ + mpz_mul_2exp(v[0], v[1], (unsigned long) 1); + mpz_set_ui(v[1], (unsigned long) 0); +} + +static void +rlggSINF(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_sub(v[0], v[1], v[0]); + mpz_add(v[1], tmp_z, v[1]); + mpz_neg(v[1], v[1]); +} + +static void +rlglSZERO(Vector v) +{ + mpz_set(tmp_z, v[0]); + mpz_add(v[0], v[0], v[1]); + mpz_sub(v[1], v[1], tmp_z); +} -- cgit v1.2.3