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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
|
(* Copyright (c) 2011-2020 Stefan Krah. All rights reserved. *)
========================================================================
Calculate (a * b) % p using the 80-bit x87 FPU
========================================================================
A description of the algorithm can be found in the apfloat manual by
Tommila [1].
The proof follows an argument made by Granlund/Montgomery in [2].
Definitions and assumptions:
----------------------------
The 80-bit extended precision format uses 64 bits for the significand:
(1) F = 64
The modulus is prime and less than 2**31:
(2) 2 <= p < 2**31
The factors are less than p:
(3) 0 <= a < p
(4) 0 <= b < p
The product a * b is less than 2**62 and is thus exact in 64 bits:
(5) n = a * b
The product can be represented in terms of quotient and remainder:
(6) n = q * p + r
Using (3), (4) and the fact that p is prime, the remainder is always
greater than zero:
(7) 0 <= q < p /\ 1 <= r < p
Strategy:
---------
Precalculate the 80-bit long double inverse of p, with a maximum
relative error of 2**(1-F):
(8) pinv = (long double)1.0 / p
Calculate an estimate for q = floor(n/p). The multiplication has another
maximum relative error of 2**(1-F):
(9) qest = n * pinv
If we can show that q < qest < q+1, then trunc(qest) = q. It is then
easy to recover the remainder r. The complete algorithm is:
a) Set the control word to 64-bit precision and truncation mode.
b) n = a * b # Calculate exact product.
c) qest = n * pinv # Calculate estimate for the quotient.
d) q = (qest+2**63)-2**63 # Truncate qest to the exact quotient.
f) r = n - q * p # Calculate remainder.
Proof for q < qest < q+1:
-------------------------
Using the cumulative error, the error bounds for qest are:
n n * (1 + 2**(1-F))**2
(9) --------------------- <= qest <= ---------------------
p * (1 + 2**(1-F))**2 p
Lemma 1:
--------
n q * p + r
(10) q < --------------------- = ---------------------
p * (1 + 2**(1-F))**2 p * (1 + 2**(1-F))**2
Proof:
~~~~~~
(I) q * p * (1 + 2**(1-F))**2 < q * p + r
(II) q * p * 2**(2-F) + q * p * 2**(2-2*F) < r
Using (1) and (7), it is sufficient to show that:
(III) q * p * 2**(-62) + q * p * 2**(-126) < 1 <= r
(III) can easily be verified by substituting the largest possible
values p = 2**31-1 and q = 2**31-2.
The critical cases occur when r = 1, n = m * p + 1. These cases
can be exhaustively verified with a test program.
Lemma 2:
--------
n * (1 + 2**(1-F))**2 (q * p + r) * (1 + 2**(1-F))**2
(11) --------------------- = ------------------------------- < q + 1
p p
Proof:
~~~~~~
(I) (q * p + r) + (q * p + r) * 2**(2-F) + (q * p + r) * 2**(2-2*F) < q * p + p
(II) (q * p + r) * 2**(2-F) + (q * p + r) * 2**(2-2*F) < p - r
Using (1) and (7), it is sufficient to show that:
(III) (q * p + r) * 2**(-62) + (q * p + r) * 2**(-126) < 1 <= p - r
(III) can easily be verified by substituting the largest possible
values p = 2**31-1, q = 2**31-2 and r = 2**31-2.
The critical cases occur when r = (p - 1), n = m * p - 1. These cases
can be exhaustively verified with a test program.
[1] http://www.apfloat.org/apfloat/2.40/apfloat.pdf
[2] http://gmplib.org/~tege/divcnst-pldi94.pdf
[Section 7: "Use of floating point"]
(* Coq proof for (10) and (11) *)
Require Import ZArith.
Require Import QArith.
Require Import Qpower.
Require Import Qabs.
Require Import Psatz.
Open Scope Q_scope.
Ltac qreduce T :=
rewrite <- (Qred_correct (T)); simpl (Qred (T)).
Theorem Qlt_move_right :
forall x y z:Q, x + z < y <-> x < y - z.
Proof.
intros.
split.
intros.
psatzl Q.
intros.
psatzl Q.
Qed.
Theorem Qlt_mult_by_z :
forall x y z:Q, 0 < z -> (x < y <-> x * z < y * z).
Proof.
intros.
split.
intros.
apply Qmult_lt_compat_r. trivial. trivial.
intros.
rewrite <- (Qdiv_mult_l x z). rewrite <- (Qdiv_mult_l y z).
apply Qmult_lt_compat_r.
apply Qlt_shift_inv_l.
trivial. psatzl Q. trivial. psatzl Q. psatzl Q.
Qed.
Theorem Qle_mult_quad :
forall (a b c d:Q),
0 <= a -> a <= c ->
0 <= b -> b <= d ->
a * b <= c * d.
intros.
psatz Q.
Qed.
Theorem q_lt_qest:
forall (p q r:Q),
(0 < p) -> (p <= (2#1)^31 - 1) ->
(0 <= q) -> (q <= p - 1) ->
(1 <= r) -> (r <= p - 1) ->
q < (q * p + r) / (p * (1 + (2#1)^(-63))^2).
Proof.
intros.
rewrite Qlt_mult_by_z with (z := (p * (1 + (2#1)^(-63))^2)).
unfold Qdiv.
rewrite <- Qmult_assoc.
rewrite (Qmult_comm (/ (p * (1 + (2 # 1) ^ (-63)) ^ 2)) (p * (1 + (2 # 1) ^ (-63)) ^ 2)).
rewrite Qmult_inv_r.
rewrite Qmult_1_r.
assert (q * (p * (1 + (2 # 1) ^ (-63)) ^ 2) == q * p + (q * p) * ((2 # 1) ^ (-62) + (2 # 1) ^ (-126))).
qreduce ((1 + (2 # 1) ^ (-63)) ^ 2).
qreduce ((2 # 1) ^ (-62) + (2 # 1) ^ (-126)).
ring_simplify.
reflexivity.
rewrite H5.
rewrite Qplus_comm.
rewrite Qlt_move_right.
ring_simplify (q * p + r - q * p).
qreduce ((2 # 1) ^ (-62) + (2 # 1) ^ (-126)).
apply Qlt_le_trans with (y := 1).
rewrite Qlt_mult_by_z with (z := 85070591730234615865843651857942052864 # 18446744073709551617).
ring_simplify.
apply Qle_lt_trans with (y := ((2 # 1) ^ 31 - (2#1)) * ((2 # 1) ^ 31 - 1)).
apply Qle_mult_quad.
assumption. psatzl Q. psatzl Q. psatzl Q. psatzl Q. psatzl Q. assumption. psatzl Q. psatzl Q.
Qed.
Theorem qest_lt_qplus1:
forall (p q r:Q),
(0 < p) -> (p <= (2#1)^31 - 1) ->
(0 <= q) -> (q <= p - 1) ->
(1 <= r) -> (r <= p - 1) ->
((q * p + r) * (1 + (2#1)^(-63))^2) / p < q + 1.
Proof.
intros.
rewrite Qlt_mult_by_z with (z := p).
unfold Qdiv.
rewrite <- Qmult_assoc.
rewrite (Qmult_comm (/ p) p).
rewrite Qmult_inv_r.
rewrite Qmult_1_r.
assert ((q * p + r) * (1 + (2 # 1) ^ (-63)) ^ 2 == q * p + r + (q * p + r) * ((2 # 1) ^ (-62) + (2 # 1) ^ (-126))).
qreduce ((1 + (2 # 1) ^ (-63)) ^ 2).
qreduce ((2 # 1) ^ (-62) + (2 # 1) ^ (-126)).
ring_simplify. reflexivity.
rewrite H5.
rewrite <- Qplus_assoc. rewrite <- Qplus_comm. rewrite Qlt_move_right.
ring_simplify ((q + 1) * p - q * p).
rewrite <- Qplus_comm. rewrite Qlt_move_right.
apply Qlt_le_trans with (y := 1).
qreduce ((2 # 1) ^ (-62) + (2 # 1) ^ (-126)).
rewrite Qlt_mult_by_z with (z := 85070591730234615865843651857942052864 # 18446744073709551617).
ring_simplify.
ring_simplify in H0.
apply Qle_lt_trans with (y := (2147483646 # 1) * (2147483647 # 1) + (2147483646 # 1)).
apply Qplus_le_compat.
apply Qle_mult_quad.
assumption. psatzl Q. auto with qarith. assumption. psatzl Q.
auto with qarith. auto with qarith.
psatzl Q. psatzl Q. assumption.
Qed.
|
