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/print.c | 321 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 321 insertions(+) create mode 100644 ic-reals-6.3/base/print.c (limited to 'ic-reals-6.3/base/print.c') diff --git a/ic-reals-6.3/base/print.c b/ic-reals-6.3/base/print.c new file mode 100644 index 0000000..8f001cd --- /dev/null +++ b/ic-reals-6.3/base/print.c @@ -0,0 +1,321 @@ +/* + * 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 + +void +print_R_Dec(Real x, int precision) +{ + force_R_Dec(x, precision); + print_R(x); +} + +void +print_R_Digs(Real x, int digitsNeeded) +{ + force_R_Digs(x, digitsNeeded); + print_R(x); +} + +/* + * The following two functions are pretty much the same and the common + * parts should be reconciled into a single function. The argument must + * have been forced so that it dereferences to a stream. + */ +void +print_R(Real x) +{ + SmallMatrix smallAccumMat; + int tmp; + mpf_t a, b, c, d, diff, mid; + int count; + Sign sign; + DigsX *digsX; + Real derefToStrm(Real); + + x = derefToStrm(x); + + if (x == NULL) { + printf("[-inf,+inf]"); + return; + } + + if (x->gen.tag.type == SIGNX) { + sign = x->signX.tag.value; + digsX = (DigsX *) x->signX.x; + } + else { + sign = SPOS; + digsX = (DigsX *) x; + } + + if (digsX->tag.type != DIGSX) + Error(FATAL, E_INT, "print_R", "real argument not a stream"); + + mpf_set_default_prec(digsX->count); + mpf_init(a); + mpf_init(b); + mpf_init(c); + mpf_init(d); + mpf_init(diff); + mpf_init(mid); + +#ifdef PACK_DIGITS + if (digsX->count <= DIGITS_PER_WORD) { + makeSmallMatrixFromDigits(smallAccumMat, digsX); + switch (sign) { + case SIGN_UNKN : + Error(FATAL, E_INT, "print_R", "sign not valid (1)\n"); + case SPOS : + break; + case SNEG : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = -smallAccumMat[0][1]; + smallAccumMat[0][1] = tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = -smallAccumMat[1][1]; + smallAccumMat[1][1] = tmp; + break; + case SINF : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = smallAccumMat[0][0] + smallAccumMat[0][1]; + smallAccumMat[0][1] = smallAccumMat[0][1] - tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = smallAccumMat[1][0] + smallAccumMat[1][1]; + smallAccumMat[1][1] = smallAccumMat[1][1] - tmp; + break; + case SZERO : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = smallAccumMat[0][0] - smallAccumMat[0][1]; + smallAccumMat[0][1] = smallAccumMat[0][1] + tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = smallAccumMat[1][0] - smallAccumMat[1][1]; + smallAccumMat[1][1] = smallAccumMat[1][1] + tmp; + break; + } +/* + printf("a=%d\n", smallAccumMat[0][0]); + printf("b=%d\n", smallAccumMat[0][1]); + printf("c=%d\n", smallAccumMat[1][0]); + printf("d=%d\n", smallAccumMat[1][1]); +*/ + mpf_set_si(a, smallAccumMat[0][0]); + mpf_set_si(b, smallAccumMat[0][1]); + mpf_set_si(c, smallAccumMat[1][0]); + mpf_set_si(d, smallAccumMat[1][1]); + } + else { +#endif + makeMatrixFromDigits(bigTmpMat, digsX); + switch (sign) { + case SIGN_UNKN : + Error(FATAL, E_INT, "print_R", "sign not valid (2)\n"); + case SPOS : + break; + case SNEG : + MPZ_SWAP(bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_neg(bigTmpMat[0][0], bigTmpMat[0][0]); + MPZ_SWAP(bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_neg(bigTmpMat[1][0], bigTmpMat[1][0]); + break; + case SINF : + mpz_set(tmpa_z, bigTmpMat[0][0]); + mpz_add(bigTmpMat[0][0], bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_sub(bigTmpMat[0][1], bigTmpMat[0][1], tmpa_z); + mpz_set(tmpa_z, bigTmpMat[1][0]); + mpz_add(bigTmpMat[1][0], bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_sub(bigTmpMat[1][1], bigTmpMat[1][1], tmpa_z); + break; + case SZERO : + mpz_set(tmpa_z, bigTmpMat[0][0]); + mpz_sub(bigTmpMat[0][0], bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_add(bigTmpMat[0][1], bigTmpMat[0][1], tmpa_z); + mpz_set(tmpa_z, bigTmpMat[1][0]); + mpz_sub(bigTmpMat[1][0], bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_add(bigTmpMat[1][1], bigTmpMat[1][1], tmpa_z); + break; + } +/* + dumpMatrix(bigTmpMat); +*/ + mpf_set_z(a, bigTmpMat[0][0]); + mpf_set_z(b, bigTmpMat[0][1]); + mpf_set_z(c, bigTmpMat[1][0]); + mpf_set_z(d, bigTmpMat[1][1]); +#ifdef PACK_DIGITS + } +#endif + + count = floor(digsX->count * log10(2.0)); + + if (mpf_cmp_ui(d,0) != 0 && mpf_cmp_ui(b,0) != 0) { + mpf_div(c, c, d); + mpf_div(a, a, b); + mpf_add(mid, c, a); + mpf_div_ui(mid, mid, 2); + mpf_sub(diff, mid, a); + mpf_abs(diff, diff); + mpf_out_str(stdout, 10, count + 2, mid); + printf(" +-"); + mpf_out_str(stdout, 10, 2, diff); + } + else { + printf("["); + if (mpf_cmp_ui(d,0) == 0) + printf("NaN,"); + else { + mpf_div(c, c, d); + mpf_out_str(stdout, 10, count + 2, c); + printf(","); + } + + if (mpf_cmp_ui(b,0) == 0) + printf("NaN]"); + else { + mpf_div(a, a, b); + mpf_out_str(stdout, 10, count + 2, a); + } + } + + mpf_clear(a); + mpf_clear(b); + mpf_clear(c); + mpf_clear(d); + mpf_clear(diff); + mpf_clear(mid); +} + +double +realToDouble(Real x) +{ + SmallMatrix smallAccumMat; + int tmp; + mpf_t a, b, c, d, diff; + double f; + Sign sign; + DigsX *digsX; + Real derefToStrm(Real); + + x = derefToStrm(x); + + if (x == NULL) + Error(FATAL, E_INT, "realToDouble", "real argument not a stream"); + + if (x->gen.tag.type == SIGNX) { + sign = x->signX.tag.value; + digsX = (DigsX *) x->signX.x; + } + else { + sign = SPOS; + digsX = (DigsX *) x; + } + + if (digsX->tag.type != DIGSX) + Error(FATAL, E_INT, "realToDouble", "real argument not a stream"); + + mpf_set_default_prec(digsX->count); + mpf_init(a); + mpf_init(b); + mpf_init(c); + mpf_init(d); + mpf_init(diff); + +#ifdef PACK_DIGITS + if (digsX->count <= DIGITS_PER_WORD) { + makeSmallMatrixFromDigits(smallAccumMat, digsX); + switch (sign) { + case SIGN_UNKN : + Error(FATAL, E_INT, "realToDouble", "sign not valid (1)\n"); + case SPOS : + break; + case SNEG : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = -smallAccumMat[0][1]; + smallAccumMat[0][1] = tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = -smallAccumMat[1][1]; + smallAccumMat[1][1] = tmp; + break; + case SINF : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = smallAccumMat[0][0] + smallAccumMat[0][1]; + smallAccumMat[0][1] = smallAccumMat[0][1] - tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = smallAccumMat[1][0] + smallAccumMat[1][1]; + smallAccumMat[1][1] = smallAccumMat[1][1] - tmp; + break; + case SZERO : + tmp = smallAccumMat[0][0]; + smallAccumMat[0][0] = smallAccumMat[0][0] - smallAccumMat[0][1]; + smallAccumMat[0][1] = smallAccumMat[0][1] + tmp; + tmp = smallAccumMat[1][0]; + smallAccumMat[1][0] = smallAccumMat[1][0] - smallAccumMat[1][1]; + smallAccumMat[1][1] = smallAccumMat[1][1] + tmp; + break; + } + mpf_set_si(a, smallAccumMat[0][0]); + mpf_set_si(b, smallAccumMat[0][1]); + mpf_set_si(c, smallAccumMat[1][0]); + mpf_set_si(d, smallAccumMat[1][1]); + } + else { +#endif + makeMatrixFromDigits(bigTmpMat, digsX); + switch (sign) { + case SIGN_UNKN : + Error(FATAL, E_INT, "realToDouble", "sign not valid (2)\n"); + case SPOS : + break; + case SNEG : + MPZ_SWAP(bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_neg(bigTmpMat[0][0], bigTmpMat[0][0]); + MPZ_SWAP(bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_neg(bigTmpMat[1][0], bigTmpMat[1][0]); + break; + case SINF : + mpz_set(tmpa_z, bigTmpMat[0][0]); + mpz_add(bigTmpMat[0][0], bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_sub(bigTmpMat[0][1], bigTmpMat[0][1], tmpa_z); + mpz_set(tmpa_z, bigTmpMat[1][0]); + mpz_add(bigTmpMat[1][0], bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_sub(bigTmpMat[1][1], bigTmpMat[1][1], tmpa_z); + break; + case SZERO : + mpz_set(tmpa_z, bigTmpMat[0][0]); + mpz_sub(bigTmpMat[0][0], bigTmpMat[0][0], bigTmpMat[0][1]); + mpz_add(bigTmpMat[0][1], bigTmpMat[0][1], tmpa_z); + mpz_set(tmpa_z, bigTmpMat[1][0]); + mpz_sub(bigTmpMat[1][0], bigTmpMat[1][0], bigTmpMat[1][1]); + mpz_add(bigTmpMat[1][1], bigTmpMat[1][1], tmpa_z); + break; + } + mpf_set_z(a, bigTmpMat[0][0]); + mpf_set_z(b, bigTmpMat[0][1]); + mpf_set_z(c, bigTmpMat[1][0]); + mpf_set_z(d, bigTmpMat[1][1]); +#ifdef PACK_DIGITS + } +#endif + if (mpf_cmp_ui(b,0) == 0 || mpf_cmp_ui(d,0) == 0) { + Error(FATAL, E_INT, "realToDouble", "NaN"); + f = 0; + } else { + mpf_div(a, a, b); + f = mpf_get_d(a); + } + mpf_clear(a); + mpf_clear(b); + mpf_clear(c); + mpf_clear(d); + mpf_clear(diff); + return f; +} + -- cgit v1.2.3