1#@ -Eprecision=108
2
3# In this Gappa file, the mathematical values are not
4# marked with a leading M but with a leading E (for exact)
5# The reason is that we have already variables starting
6# with the letter M
7
8
9# Rounding operators and sequences definition
10
11@double = float<ieee_64,ne>;
12@Add22 = add_rel<102>;
13@Mul22 = mul_rel<102>;
14
15
16# Instantiating of constants
17SQRTPOLYC0 = double(_SQRTPOLYC0);
18SQRTPOLYC1 = double(_SQRTPOLYC1);
19SQRTPOLYC2 = double(_SQRTPOLYC2);
20SQRTPOLYC3 = double(_SQRTPOLYC3);
21SQRTPOLYC4 = double(_SQRTPOLYC4);
22
23
24# Definition of what is a sqrt
25# We start with the reciprocal square root. The square root is its reciprocal.
26# The reduced argument m is the square of its root.
27
28ESqrtm = 1 / ERecSqrtm;
29m = ESqrtm * ESqrtm;
30ERecM = 1 / m;
31
32# Translation of the code
33r0 double= SQRTPOLYC0 + m * (SQRTPOLYC1 + m * (SQRTPOLYC2 + m * (SQRTPOLYC3 + m * SQRTPOLYC4)));
34
35r1 double= (1/2) * r0 * (3 - m * (r0 * r0));
36r2 double= (1/2) * r1 * (3 - m * (r1 * r1));
37
38R2Sqhl = r2 * r2;
39R2PHr2hl = (3/2) * r2;
40
41MMr2hl = m * r2;
42
43MMr2Chl = Mul22(MMr2hl, R2Sqhl);
44
45MHmMr2Chl = -(1/2) * MMr2Chl;
46
47R3hl = Add22(R2PHr2hl, MHmMr2Chl);
48
49R3Sqhl = Mul22(R3hl,R3hl);
50
51MMr3Sqhl = Mul22(m,R3Sqhl);
52
53# The next translation of the code is correct only if we can prove that
54# mMr3Sqh is exactly 1 in each case
55# We show this below
56
57T1 = 1 + -(1/2) * (MMr3Sqhl - 1);
58
59R4hl = Mul22(R3hl,T1);
60
61Srtmhl = Mul22(m,R4hl);
62
63# Show that mMr3Sqh is exactly 1 in each case
64# We round the double-double result MMr3Sqhl to double
65# precision and show that the result is in the singleton
66# interval [1,1]
67
68mMr3Sqh = double(MMr3Sqhl);
69
70# Mathematical equivalents (marked with E for exact)
71
72Er0 = SQRTPOLYC0 + m * (SQRTPOLYC1 + m * (SQRTPOLYC2 + m * (SQRTPOLYC3 + m * SQRTPOLYC4)));
73
74Er1 = (1/2) * Er0 * (3 - m * (Er0 * Er0));
75Er2 = (1/2) * Er1 * (3 - m * (Er1 * Er1));
76
77#Definition of the errors
78
79epsilon = (Srtmhl - ESqrtm) / ESqrtm;
80
81epsilonApproxMath0 = (Er0 - ERecSqrtm) / ERecSqrtm;
82epsilonApproxMath1 = (Er1 - ERecSqrtm) / ERecSqrtm;
83epsilonApproxMath2 = (Er2 - ERecSqrtm) / ERecSqrtm;
84
85epsilonApproxArith0 = (r0 - ERecSqrtm) / ERecSqrtm;
86epsilonApproxArith1 = (r1 - ERecSqrtm) / ERecSqrtm;
87epsilonApproxArith2 = (r2 - ERecSqrtm) / ERecSqrtm;
88epsilonApproxArith3 = (R3hl - ERecSqrtm) / ERecSqrtm;
89epsilonApproxArith4 = (R4hl - ERecSqrtm) / ERecSqrtm;
90
91epsilon0 = (r0 - Er0) / Er0;
92epsilon1 = (r1 - Er1) / Er1;
93epsilon2 = (r2 - Er2) / Er2;
94
95epsilon3 = ((R2PHr2hl + MHmMr2Chl) - ERecSqrtm) / ERecSqrtm;
96epsilon4 = ((R3hl * T1) - ERecSqrtm) / ERecSqrtm;
97
98epsilon5 = ((T1 - 3/2) - (-1/2 * m * R3hl * R3hl)) / (-1/2 * m * R3hl * R3hl);
99epsilon6 = (MMr3Sqhl - (m * (R3hl * R3hl))) / (m * (R3hl * R3hl));
100
101epsilon7 = (R3Sqhl - ERecM) / ERecM;
102epsilon8 = (MMr3Sqhl - 1) / 1;
103
104epsilon9 = ((R3hl * R3hl) - ERecM) / ERecM;
105epsilon10 = ((m * R3Sqhl) - 1) / 1;
106
107epsilon11 = ((m * R4hl) - ESqrtm) / ESqrtm;
108
109epsilonOp0 = (MMr2Chl - (MMr2hl * R2Sqhl)) / (MMr2hl * R2Sqhl);
110epsilonOp1 = (R3hl - (R2PHr2hl + MHmMr2Chl)) / (R2PHr2hl + MHmMr2Chl);
111epsilonOp2 = (R3Sqhl - (R3hl * R3hl)) / (R3hl * R3hl);
112epsilonOp3 = (MMr3Sqhl - (m * R3Sqhl)) / (m * R3Sqhl);
113epsilonOp4 = (R4hl - (R3hl * T1)) / (R3hl * T1);
114epsilonOp5 = (Srtmhl - (m * R4hl)) / (m * R4hl);
115
116
117#Logical implication to prove
118
119{
120(
121   ERecSqrtm in [_ERecSqrtmMin,_ERecSqrtmMax]
122/\ epsilonApproxMath0 in [-_epsilonApprox,_epsilonApprox]
123)
124->
125(
126   epsilon in [-1b-100,1b-100]               # Accuracy bound
127/\ mMr3Sqh in [1,1]                          # For showing that mMr3Sqh is always equal to 1 (see above)
128)
129}
130
131# Hints for bounding the variables in the code (not the errors)
132
133Er0 -> ERecSqrtm * epsilonApproxMath0 + ERecSqrtm;
134Er1 -> ERecSqrtm * epsilonApproxMath1 + ERecSqrtm;
135Er2 -> ERecSqrtm * epsilonApproxMath2 + ERecSqrtm;
136
137r0 -> Er0 * epsilon0 + Er0;
138r1 -> Er1 * epsilon1 + Er1;
139r2 -> Er2 * epsilon2 + Er2;
140
141(R2PHr2hl + MHmMr2Chl) -> ERecSqrtm * epsilon3 + ERecSqrtm;
142
143MMr3Sqhl -> 1 + epsilon8;
144
145
146# Hints for the Newton iteration
147
148epsilonApproxMath1 -> -3/2 * epsilonApproxMath0 * epsilonApproxMath0 -
149                      1/2 * epsilonApproxMath0 * epsilonApproxMath0 * epsilonApproxMath0;
150
151epsilonApproxMath2 -> -3/2 * epsilonApproxMath1 * epsilonApproxMath1 -
152                      1/2 * epsilonApproxMath1 * epsilonApproxMath1 * epsilonApproxMath1;
153
154
155# Hints for the Newton iteration with integrated arithmetical errors
156
157epsilon3 -> -3/2 * epsilonApproxArith2 * epsilonApproxArith2
158            -1/2 * epsilonApproxArith2 * epsilonApproxArith2 * epsilonApproxArith2
159            -1/2 * epsilonOp0
160            -3/2 * epsilonApproxArith2 * epsilonOp0
161            -3/2 * epsilonApproxArith2 * epsilonApproxArith2 * epsilonOp0
162            -1/2 * epsilonApproxArith2 * epsilonApproxArith2 * epsilonApproxArith2 * epsilonOp0;
163
164epsilon5 -> epsilon6;
165
166epsilon4 -> -3/2 * epsilonApproxArith3 * epsilonApproxArith3
167            -1/2 * epsilonApproxArith3 * epsilonApproxArith3 * epsilonApproxArith3
168            -1/2 * epsilon5
169            -3/2 * epsilon5 * epsilonApproxArith3
170            -3/2 * epsilon5 * epsilonApproxArith3 * epsilonApproxArith3
171            -1/2 * epsilon5 * epsilonApproxArith3 * epsilonApproxArith3 * epsilonApproxArith3;
172
173epsilon9 -> 2 * epsilonApproxArith3 + epsilonApproxArith3 * epsilonApproxArith3;
174
175epsilon10 -> epsilon7;
176
177epsilon11 -> epsilonApproxArith4;
178
179# Meta-Hints
180
181r0 ~ Er0;
182r1 ~ Er1;
183r2 ~ Er2;
184