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/epsDel.c | 133 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 ic-reals-6.3/base/epsDel.c (limited to 'ic-reals-6.3/base/epsDel.c') diff --git a/ic-reals-6.3/base/epsDel.c b/ic-reals-6.3/base/epsDel.c new file mode 100644 index 0000000..3f4fe08 --- /dev/null +++ b/ic-reals-6.3/base/epsDel.c @@ -0,0 +1,133 @@ +/* + * 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" + +/* + * Functions dealing with the epsilon-delta analysis for matrices and + * tensors. + */ + +/* +static mpz_t vec0_sum, vec1_sum, vec2_sum, vec3_sum; +*/ + +#define vec0_sum tmpc_z +#define vec1_sum tmpd_z +#define vec2_sum tmpe_z +#define vec3_sum tmpf_z + +static int EpsDelVectorPair(Vector, Vector, mpz_t, mpz_t, int); + +void +initEpsDel() +{ +/* + * These aren't needed anymore as temporary storage is shared amongst + * the entire library + * + mpz_init(vec0_sum); + mpz_init(vec1_sum); + mpz_init(vec2_sum); + mpz_init(vec3_sum); +*/ +} + +/* + * Given a matrix and an epsilon (number of digits needed), this returns + * the number of digits the matrix must consume from its argument + */ +int +epsDelMatrix(Matrix mat, int epsilon) +{ + mpz_add(vec0_sum, mat[0][0], mat[0][1]); + mpz_add(vec1_sum, mat[1][0], mat[1][1]); + return EpsDelVectorPair(mat[0], mat[1], vec0_sum, vec1_sum, epsilon); +} + +/* + * As for matrices above, but this gives the number of digits needed from + * the two arguments of a tensor. + */ +void +epsDelTensor(Tensor ten, int epsilon, int *nLeft, int *nRight) +{ + int r, s, p, q; + + mpz_add(vec0_sum, ten[0][0], ten[0][1]); + mpz_abs(vec0_sum, vec0_sum); + mpz_add(vec1_sum, ten[1][0], ten[1][1]); + mpz_abs(vec1_sum, vec1_sum); + mpz_add(vec2_sum, ten[2][0], ten[2][1]); + mpz_abs(vec2_sum, vec2_sum); + mpz_add(vec3_sum, ten[3][0], ten[3][1]); + mpz_abs(vec3_sum, vec3_sum); + + r = EpsDelVectorPair(ten[0], ten[2], vec0_sum, vec2_sum, epsilon); + s = EpsDelVectorPair(ten[1], ten[3], vec1_sum, vec3_sum, epsilon); + p = EpsDelVectorPair(ten[0], ten[1], vec0_sum, vec1_sum, epsilon); + q = EpsDelVectorPair(ten[2], ten[3], vec2_sum, vec3_sum, epsilon); + + *nLeft = MIN(r,s); + *nRight = MIN(p,q); +} + +/* + * As above for tensors, but we are interested only in the number + * of arguments needed from the x (left) argument. + */ +int +epsDelTensorX(Tensor ten, int epsilon) +{ + int r, s; + + mpz_add(vec0_sum, ten[0][0], ten[0][1]); + mpz_abs(vec0_sum, vec0_sum); + mpz_add(vec1_sum, ten[1][0], ten[1][1]); + mpz_abs(vec1_sum, vec1_sum); + mpz_add(vec2_sum, ten[2][0], ten[2][1]); + mpz_abs(vec2_sum, vec2_sum); + mpz_add(vec3_sum, ten[3][0], ten[3][1]); + mpz_abs(vec3_sum, vec3_sum); + + r = EpsDelVectorPair(ten[0], ten[2], vec0_sum, vec2_sum, epsilon); + s = EpsDelVectorPair(ten[1], ten[3], vec1_sum, vec3_sum, epsilon); + + return MIN(r,s); +} + +static int +EpsDelVectorPair(Vector vec0, Vector vec1, + mpz_t vec0_sum, + mpz_t vec1_sum, + int epsilon) +{ + int logdet, logalsq, bitcount; + + mpz_mul(tmpa_z, vec0[0], vec1[1]); + mpz_mul(tmpb_z, vec0[1], vec1[0]); + mpz_sub(tmpa_z, tmpa_z, tmpb_z); + mpz_abs(tmpa_z, tmpa_z); + logdet = mpz_sizeinbase(tmpa_z, 2); + + if (mpz_cmp(vec0_sum, vec1_sum) >= 0) { + bitcount = mpz_popcount(vec0_sum); + mpz_mul(tmpa_z, vec0_sum, vec0_sum); + } + else { + bitcount = mpz_popcount(vec1_sum); + mpz_mul(tmpa_z, vec1_sum, vec1_sum); + } + logalsq = mpz_sizeinbase(tmpa_z, 2); + if (bitcount == 1) + return epsilon + logdet - logalsq; + else + return epsilon + logdet - logalsq - 1; +} -- cgit v1.2.3