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