1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
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 <stdio.h>
#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;
}
|
