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
|
// Implement WorstCase functions to compute the worst case for x mod C, with
// the exponent of x ranges from emin to emax, and precision of x is p.
// Adapted to Sollya from the Maple function in
// J-M. Muller, "Elementary Functions", 3rd ed, Section 11.3.2.
//
// Some examples:
//
// 1) Worst case for trig range reduction fast passes:
//
// Single precision
// > WorstCase(24, -6, 32, pi/32, 128);
// numbermin : 10741887
// expmin : 7
// Worst case: 0x1.47d0fep30
// numberofdigits : -32.888
//
// Double precision
// > WorstCase(53, -8, 32, pi/128, 256);
// numbermin : 6411027962775774
// expmin : -53
// Worst case: 0x1.6c6cbc45dc8dep-1
// numberofdigits : -66.4867
//
// 2) Worst case for exponential function range reduction:
//
// Single precision
// > WorstCase(24, 1, 8, log(2), 128);
// numbermin : 2907270
// expmin : -22
// Worst case: 0x1.62e43p-1
// numberofdigits : -28.9678
//
// Double precision
// > WorstCase(53, 0, 11, log(2), 128);
// numbermin : 7804143460206699
// expmin : -51
// Worst case: 0x1.bb9d3beb8c86bp1
// numberofdigits : -57.4931
//
verbosity=0;
procedure WorstCase(p, emin, emax, C, ndigits) {
epsilonmin = 12345.0;
Digits = ndigits;
powerofBoverC = 2^(emin - p) / C;
for e from emin - p + 1 to emax - p + 1 do {
powerofBoverC = 2 * powerofBoverC;
a = floor(powerofBoverC);
Plast = a;
r = round(1/(powerofBoverC - a), ndigits, RN);
a = floor(r);
Qlast = 1;
Q = a;
P = Plast * a + 1;
while (Q < 2^p - 1) do {
r = round(1/(r - a), ndigits, RN);
a = floor(r);
NewQ = Q * a + Qlast;
NewP = P * a + Plast;
Qlast = Q;
Plast = P;
Q = NewQ;
P = NewP;
};
epsilon = C * abs(Plast - Qlast * powerofBoverC);
if (epsilon < epsilonmin) then {
epsilonmin = epsilon;
numbermin = Qlast;
expmin = e;
};
};
display=decimal!;
print("numbermin : ", numbermin);
print("expmin : ", expmin);
display=hexadecimal!;
print("Worst case: ", numbermin * 2^expmin);
display=decimal!;
ell = round(log2(epsilonmin), ndigits, RN);
print("numberofdigits : ", ell);
};
|
