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/force_R.c | 231 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100644 ic-reals-6.3/base/force_R.c (limited to 'ic-reals-6.3/base/force_R.c') diff --git a/ic-reals-6.3/base/force_R.c b/ic-reals-6.3/base/force_R.c new file mode 100644 index 0000000..c67c0f3 --- /dev/null +++ b/ic-reals-6.3/base/force_R.c @@ -0,0 +1,231 @@ +/* + * 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 + +/* + * This is the toplevel function called by the user to force information + * from a real stream. + */ +void +force_R_Digs(Real x, int digitsNeeded) +{ + DigsX *digsX; + void runStack(); + void force_To_Alt_Entry(); + + if (x->gen.tag.type == ALT) { + if (x->alt.redirect == NULL) { + PUSH_2(force_To_Alt_Entry, x); + runStack(); + } + force_R_Digs(x->alt.redirect, digitsNeeded); + return; + } + + if (x->gen.tag.type == CLOSURE) { + if (x->cls.redirect == NULL) { + PUSH_2(x->cls.force, x); + runStack(); + } + force_R_Digs(x->cls.redirect, digitsNeeded); + return; + } + + x = makeStream(x); + + /* + * So now we know that x is a real stream. If x is signed and + * the sign has not been determined, we force it. + */ + if (x->gen.tag.type == SIGNX) { + if (x->signX.tag.value == SIGN_UNKN) { + PUSH_2(x->signX.force, x); + runStack(); + } + digsX = (DigsX *) x->signX.x; + } + else + digsX = (DigsX *) x; + + if (digsX->tag.type != DIGSX) + Error(FATAL, E_INT, "force_R_Digs", "badly formed real"); + + /* + * Now, if there is not enough information available, we force + * the number of digits needed. + */ + if (digsX->count < (unsigned int)digitsNeeded) { + PUSH_3(digsX->force, digsX, digitsNeeded - digsX->count); + runStack(); + } +} + +/* force_R_Dec + * July 2000, Marko Krznaric + * + * force_R_Digs(x,digitsNeeded) is a function which forces an + * emission of digitsNeeded digit matrices from x (of the + * type Real). + * + * On the other hand, force_R_Dec(x,decimalPrecision) guarantees + * that enough digit matrices, say digitsNeeded, will be emitted + * from a real number x in order to have required absolute + * decimal precision, i.e. it guarantees that the result will be + * accurate within 10^(-decimalPrecision). + * + * NOTES: + * - initGuess = an initial guess (minimum number) of digit + * matrices which has to be emitted from a real number x. + * - digitsNeeded = number of digit matrices required. + * - e = number of 'bad' digits. + * - using function retrieveInfo, we can extract the sign (=sign), + * the number of digits emitted so far (=count) and the + * compressed digits (=digits) for a Real x. + * - 3.322 = an upper bound for log2(10). + * + * + * PROBLEMS: + * - should decimalPrecision & digitsNeeded be of type long int? + * + */ +void +force_R_Dec(Real x, int decimalPrecision) +{ + int initGuess = ceil(3.322 * decimalPrecision) + 2; + int digitsNeeded; + int e = 0; + + mpz_t digits; + Sign sign; + int count; + + mpz_init(digits); + + digitsNeeded = initGuess; + force_R_Digs(x, digitsNeeded); + retrieveInfo(x, &sign, &count, digits); + + switch (sign) { + case SZERO: + /* + * SZERO: every digit matrix will half the interval + * (starting with [-1,1]). Easy to determine digitsNeeded. + */ + digitsNeeded = initGuess - 1; + break; + + case SPOS: + case SNEG: + /* + * SPOS, SNEG: e, the number of 'bad' digits is actually + * the number of leading 1s in the binary representation + * of digits (the value of compressed digits). + */ + while (digitsNeeded < forceDecUpperBound) { + force_R_Digs(x, digitsNeeded); + retrieveInfo(x, &sign, &count, digits); + if (sign == SNEG) + mpz_neg(digits, digits); + e = leadingOnes (digits); + digitsNeeded = 2 * e + 2 + initGuess; + if (count > e) + /* not all of the extracted digit matrices are 'bad', + * i.e. the interval is bounded. Therefore, we can + * leave the loop. + */ + break; + } + if (digitsNeeded >= forceDecUpperBound) + Error(FATAL, E_INT, "force_R_Dec", + "forceDecUpperBound reached (1)"); + break; + + case SINF: + /* SINF: we are still dealing with the unbounded interval, + * i.e. the interval contains infinity, as far as digits + * (the value of compressed digits matrices) is either + * -1, 0 or 1. As soon as we get something greater (in + * absolute value) than 1, the interval doesn't contain + * the infinity, and we can calculate digitsNeeded. + */ + while (digitsNeeded < forceDecUpperBound) { + force_R_Digs(x, digitsNeeded); + retrieveInfo(x, &sign, &count, digits); + if (mpz_cmpabs_ui(digits, 1) > 0) { + e = count - mpz_sizeinbase(digits, 2); + digitsNeeded = 2 * e + 2 + initGuess; + break; + } + else { + e = count; + digitsNeeded = 2 * e + 2+ initGuess; + } + } + if (digitsNeeded >= forceDecUpperBound) + Error(FATAL, E_INT, "force_R_Dec", + "forceDecUpperBound reached (2)"); + break; + + case SIGN_UNKN: + Error(FATAL, E_INT, "force_R_Dec", "argument is signed"); + break; + + default: + Error(FATAL, E_INT, "force_R_Dec", "bad sign"); + break; + } + + force_R_Digs(x, digitsNeeded); + mpz_clear(digits); +} + +/* leadingOnes + * July 2000, Marko Krznaric + * + * For input variable c (of type mpz_t) check how many 1s + * are leading the binary representation of c. + * + * We use this function as an auxiliary function to + * force_R_Dec. We have to check how many D+ (D-) digit + * matrices follow S+ (S-). This seems to be the fastest + * way to do it. + * + * If c is zero or negative - result = 0. + * If there are no zeros at all - result = binary size of c. + * Otherwise, we scan from left-most digit and count until + * we reach 0. + * + * NOTES: + * - The index of the right-most digit is 0, while the index + * of the left-most digit is (size-1). + */ +int +leadingOnes(mpz_t c) +{ + int size = mpz_sizeinbase(c,2); + int i; + int count = 0; + + if (mpz_sgn(c) <= 0) + return 0; + + if (mpz_scan0(c, 0) == (long unsigned int)size) + return size; + + for (i = size - 1; i >= 0; i--) { + if (mpz_scan1(c, i) == (long unsigned int)i) + count++; + else + return count; + } + return 0; +} -- cgit v1.2.3