1
2c3 = <float64ne> (_c3);
3c4 = <float64ne> (_c4);
4c5 = <float64ne> (_c5);
5c6 = <float64ne> (_c6);
6c7 = <float64ne> (_c7);
7
8
9E = 0; #MAPLE
10
11zh = <float64ne> (Z);
12zl = Z - zh; #MAPLE
13
14polyHorner <float64ne>= c3 + zh * (c4 + zh * (c5 + zh * (c6 + zh * c7)));
15
16ZhSquarehl = zh * zh; #MAPLE
17zhSquareh = <float64ne> (ZhSquarehl);
18zhSquarel = <float64ne> (ZhSquarehl - zhSquareh);
19
20zhSquareHalfh = zhSquareh * (-0.5); #MAPLE
21zhSquareHalfl = zhSquarel * (-0.5); #MAPLE
22ZhSquareHalfhl = ZhSquarehl * (-0.5); #MAPLE
23
24ZhCube <float64ne>= (zh * zhSquareh);
25polyUpper <float64ne>= polyHorner * ZhCube;
26
27temp = <float64ne> (zh * zl);
28T1hl = polyUpper - temp; #MAPLE
29t1h = <float64ne> (T1hl);
30t1l = <float64ne> (T1hl - t1h);
31
32T2 = Z + ZhSquareHalfhl; #MAPLE
33t2h = <float64ne> (T2hl);
34t2l = <float64ne> (T2hl - t2h);
35
36PE = T2hl + T1hl; #MAPLE
37ph = <float64ne> (Phl);
38pl = <float64ne> (Phl - ph);
39
40
41#We can simplify the computations in the function in this case as we know that
42#all operations (add, mult) on 0 (as a double double) are exact.
43
44Loghm = Phl; #MAPLE
45
46logh = <float64ne> (Loghm);
47logm = <float64ne> (Loghm - logh);
48
49#Mathematical definition of the logarithm and the polynomial
50
51Phigher = (c3 + Z * (c4 + Z * (c5 + Z * (c6 + Z * c7)))); #MAPLE
52ZZZ = Z*Z*Z; #MAPLE
53ZZZPhigher = ZZZ * Phigher; #MAPLE
54HZZ = (-0.5*Z*Z); #MAPLE
55ZpHZZ = Z + HZZ; #MAPLE
56P = ZpHZZ + ZZZPhigher; #MAPLE
57
58#We apply the same simplification on the mathematical definition of the log
59Log = Log1pZ; #MAPLE
60
61
62#Additional useful definitions
63
64ZZ = Z*Z; #MAPLE
65ZZZPhigherPzhzl = ZZZPhigher - zh * zl; #MAPLE
66
67HZ = -0.5*Z; #MAPLE
68
69Flzhzl = temp; #MAPLE
70
71{
72(T2hl - T2) / T2 in [-1b-103,1b-103]
73/\ (Phl - PE) / PE in [-1b-103,1b-103]
74/\ Z in [_zmin,_zmax]
75/\ (P - Log1pZ) / Log1pZ in [-_epsilonApproxQuick,_epsilonApproxQuick]
76/\ ((logh + logm) - Loghm) / Loghm in [-1b-106,1b-106]
77->
78((logh + logm) - Log) / Log in [-5b-65,5b-65]
79}
80
81T2hl - T2 -> ((T2hl - T2) / T2) * T2;
82T2hl -> (T2hl - T2) + T2;
83
84Phl - PE -> ((Phl - PE) / PE) * PE;
85Phl -> (Phl - PE) + PE;
86
87
88(ZhSquarehl - ZZ) / ZZ -> 2 * ((zh - Z) / Z) + ((zh - Z) / Z) * ((zh - Z) / Z);
89
90(zhSquareh - ZZ) / ZZ -> ((ZhSquarehl - ZZ) / ZZ) + ((zhSquareh - ZhSquarehl) / ZZ);
91
92(zhSquareh - ZhSquarehl) / ZZ -> ((zhSquareh - ZhSquarehl) / ZhSquarehl) * (ZhSquarehl / ZZ);
93
94ZhSquarehl / ZZ -> ((ZhSquarehl - ZZ) / ZZ) + 1;
95
96(ZhCube - ZZZ) / ZZZ -> (((zh * zhSquareh) - ZZZ) / ZZZ) + ((ZhCube - (zh * zhSquareh)) / ZZZ);
97
98((zh * zhSquareh) - ZZZ) / ZZZ -> (1 + ((zh - Z) / Z)) * (1 + ((zhSquareh - ZZ) / ZZ)) - 1;
99
100((ZhCube - (zh * zhSquareh)) / ZZZ) -> ((ZhCube - (zh * zhSquareh)) / (zh * zhSquareh)) * (((zh - Z) / Z) + 1) * (((zhSquareh - ZZ) / ZZ) + 1);
101
102polyHorner / Phigher -> ((polyHorner - Phigher) / Phigher) + 1;
103
104(polyUpper - ZZZPhigher) / ZZZPhigher -> ((polyHorner - Phigher) / Phigher) + ((ZhCube - ZZZ) / ZZZ) * (polyHorner / Phigher) +
105					  + ((polyUpper - (polyHorner * ZhCube)) / (polyHorner * ZhCube)) * (polyHorner / Phigher) +
106					  + ((ZhCube - ZZZ) / ZZZ) * ((polyUpper - (polyHorner * ZhCube)) / (polyHorner * ZhCube)) *
107					    (polyHorner / Phigher);
108
109
110((ZhSquareHalfhl - (zh * zl)) - HZZ) / HZZ -> - ((zh - Z) / Z) * ((zh - Z) / Z);
111
112(ZhSquareHalfhl - HZZ) / HZZ -> (ZhSquarehl - ZZ) / ZZ;
113
114((T2hl - (zh * zl)) - ZpHZZ) / ZpHZZ -> ((HZ * (((ZhSquareHalfhl - (zh * zl)) - HZZ) / HZZ)) + ((T2hl - T2) / T2)
115                                        + (HZ * ((T2hl - T2) / T2))
116					+ (HZ * ((ZhSquareHalfhl - HZZ) / HZZ) * ((T2hl - T2) / T2))) / (1 + HZ);
117
118(PE - P) / P -> (((1 + HZ) * (((T2hl - (zh * zl)) - ZpHZZ) / ZpHZZ)) +
119		((1 + ((zh - Z) / Z)) * (Z * ((zh - Z) / Z)) * ((Flzhzl - (zh * zl)) / (zh * zl)))
120	        + (ZZ * Phigher * ((polyUpper - ZZZPhigher) / ZZZPhigher))) / (1 + HZ + ZZ * Phigher);
121
122(Phl - P) / P -> ((PE - P) / P) + ((((PE - P) / P) + 1) * ((Phl - PE) / PE));
123
124(Loghm - Log) / Log -> ((Loghm - P) / P) + ((P - Log) / Log) + ((Loghm - P) / P) * ((P - Log) / Log);
125
126(((logh + logm) - Log) / Log) -> (((logh + logm) - Loghm) / Loghm) + ((Loghm - Log) / Log) + (((logh + logm) - Loghm) / Loghm) * ((Loghm - Log) / Log);
127