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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
/*
* 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"
static mpz_t k;
Real Pi;
static void piCont();
static void nextMatrix(Tensor, int);
void
initPi()
{
Real sq, cls;
DigsX *digsX;
TenXY *tenXY;
static int doneInit;
void force_To_DigsX_From_TenXY_Entry();
if (!doneInit) {
registerForceFunc(piCont, "piCont", 2);
doneInit++;
}
mpz_init_set_str(k, "10939058860032000", 10);
sq = sqrt_QInt(10005, 1);
cls = (Real) allocCls(piCont, (void *) 0);
cls->gen.tag.isSigned = FALSE;
tenXY = (TenXY *) div_R_R(sq, cls);
tenXY->forceY = piCont;
digsX = allocDigsX();
digsX->x = (Real) tenXY;
digsX->force = force_To_DigsX_From_TenXY_Entry;
if (DAVINCI) {
beginGraphUpdate();
newEdgeToOnlyChild(digsX, tenXY);
endGraphUpdate();
}
Pi = (Real) digsX;
}
/*
* It might be possible compute the next matrix in the sequence
* in terms of the previous matrix. This would mean having much more
* state but presumably much quicker.
*/
static void
piCont()
{
TenXY *tenXY;
Cls *cls;
int digitsNeeded;
int n;
tenXY = (TenXY *) POP;
digitsNeeded = (int) POP;
cls = (Cls *) tenXY->y;
n = (int) cls->userData;
/*
* The two constants below have been determined empirically. The second
* may be off by one (ie 47) in some circumstances.
*/
while (digitsNeeded > 0) {
nextMatrix(tenXY->ten, n);
if (n == 0)
digitsNeeded -= 27;
else
digitsNeeded -= 46;
n +=1;
}
cls->userData = (void *) n;
}
static void
nextMatrix(Tensor ten, int n)
{
SmallMatrix smallMat;
if (n == 0) {
smallMat[0][0] = 6795705;
smallMat[0][1] = 213440;
smallMat[1][0] = 6795704;
smallMat[1][1] = 213440;
multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat);
multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat);
}
else {
/* b = 2n - 1 */
mpz_set_ui(tmpb_z, n);
mpz_mul_2exp(tmpb_z, tmpb_z, 1);
mpz_sub_ui(tmpb_z, tmpb_z, 1);
/* tmp = 6n - 5, tmpf_z = 6n */
mpz_set_ui(tmpa_z, n);
mpz_mul_ui(tmpa_z, tmpa_z, 6);
mpz_set(tmpf_z, tmpa_z);
mpz_sub_ui(tmpa_z, tmpa_z, 5);
/* tmpf_z = 6n - 1 */
mpz_sub_ui(tmpf_z, tmpf_z, 1);
/* b = (2n - 1) * (6n - 5) * (6n - 1) */
mpz_mul(tmpb_z, tmpb_z, tmpa_z);
mpz_mul(tmpb_z, tmpb_z, tmpf_z);
/* c = b*(545140134n + 13591409) */
mpz_set_ui(tmpa_z, 545140134);
mpz_mul_ui(tmpa_z, tmpa_z, n);
mpz_add_ui(tmpa_z, tmpa_z, 13591409);
mpz_mul(tmpc_z, tmpb_z, tmpa_z);
/* d = b * (n + 1) */
mpz_set_ui(tmpd_z, n + 1);
mpz_mul(tmpd_z, tmpd_z, tmpb_z);
/* e = 10939058860032000*n^4 */
mpz_set_ui(tmpe_z, n);
mpz_mul(tmpe_z, tmpe_z, tmpe_z);
mpz_mul(tmpe_z, tmpe_z, tmpe_z);
mpz_mul(tmpe_z, tmpe_z, k);
mpz_sub(bigTmpMat[0][0], tmpe_z, tmpd_z);
mpz_set(bigTmpMat[1][1], bigTmpMat[0][0]);
mpz_add(bigTmpMat[0][1], tmpe_z, tmpd_z);
mpz_set(bigTmpMat[1][0], bigTmpMat[0][1]);
mpz_sub(bigTmpMat[0][0], bigTmpMat[0][0], tmpc_z);
mpz_add(bigTmpMat[0][1], bigTmpMat[0][1], tmpc_z);
mpz_sub(bigTmpMat[1][0], bigTmpMat[1][0], tmpc_z);
mpz_add(bigTmpMat[1][1], bigTmpMat[1][1], tmpc_z);
multVectorPairTimesMatrix(ten[0], ten[1], bigTmpMat);
multVectorPairTimesMatrix(ten[2], ten[3], bigTmpMat);
}
}
|
