Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(41461)

Side by Side Diff: Modules/_decimal/libmpdec/literature/mulmod-64.txt

Issue 7652: Merge C version of decimal into py3k.
Patch Set: Created 7 years, 6 months ago
Left:
Right:
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View unified diff | Download patch
OLDNEW
(Empty)
1
2
3 (* Copyright (c) 2011 Stefan Krah. All rights reserved. *)
4
5
6 ==========================================================================
7 Calculate (a * b) % p using special primes
8 ==========================================================================
9
10 A description of the algorithm can be found in the apfloat manual by
11 Tommila [1].
12
13
14 Definitions:
15 ------------
16
17 In the whole document, "==" stands for "is congruent with".
18
19 Result of a * b in terms of high/low words:
20
21 (1) hi * 2**64 + lo = a * b
22
23 Special primes:
24
25 (2) p = 2**64 - z + 1, where z = 2**n
26
27 Single step modular reduction:
28
29 (3) R(hi, lo) = hi * z - hi + lo
30
31
32 Strategy:
33 ---------
34
35 a) Set (hi, lo) to the result of a * b.
36
37 b) Set (hi', lo') to the result of R(hi, lo).
38
39 c) Repeat step b) until 0 <= hi' * 2**64 + lo' < 2*p.
40
41 d) If the result is less than p, return lo'. Otherwise return lo' - p.
42
43
44 The reduction step b) preserves congruence:
45 -------------------------------------------
46
47 hi * 2**64 + lo == hi * z - hi + lo (mod p)
48
49 Proof:
50 ~~~~~~
51
52 hi * 2**64 + lo = (2**64 - z + 1) * hi + z * hi - hi + lo
53
54 = p * hi + z * hi - hi + lo
55
56 == z * hi - hi + lo (mod p)
57
58
59 Maximum numbers of step b):
60 ---------------------------
61
62 # To avoid unneccessary formalism, define:
63
64 def R(hi, lo, z):
65 return divmod(hi * z - hi + lo, 2**64)
66
67 # For simplicity, assume hi=2**64-1, lo=2**64-1 after the
68 # initial multiplication a * b. This is of course impossible
69 # but certainly covers all cases.
70
71 # Then, for p1:
72 hi=2**64-1; lo=2**64-1; z=2**32
73 p1 = 2**64 - z + 1
74
75 hi, lo = R(hi, lo, z) # First reduction
76 hi, lo = R(hi, lo, z) # Second reduction
77 hi * 2**64 + lo < 2 * p1 # True
78
79 # For p2:
80 hi=2**64-1; lo=2**64-1; z=2**34
81 p2 = 2**64 - z + 1
82
83 hi, lo = R(hi, lo, z) # First reduction
84 hi, lo = R(hi, lo, z) # Second reduction
85 hi, lo = R(hi, lo, z) # Third reduction
86 hi * 2**64 + lo < 2 * p2 # True
87
88 # For p3:
89 hi=2**64-1; lo=2**64-1; z=2**40
90 p3 = 2**64 - z + 1
91
92 hi, lo = R(hi, lo, z) # First reduction
93 hi, lo = R(hi, lo, z) # Second reduction
94 hi, lo = R(hi, lo, z) # Third reduction
95 hi * 2**64 + lo < 2 * p3 # True
96
97
98 Step d) preserves congruence and yields a result < p:
99 -----------------------------------------------------
100
101 Case hi = 0:
102
103 Case lo < p: trivial.
104
105 Case lo >= p:
106
107 lo == lo - p (mod p) # result is congruent
108
109 p <= lo < 2*p -> 0 <= lo - p < p # result is in the correct range
110
111 Case hi = 1:
112
113 p < 2**64 /\ 2**64 + lo < 2*p -> lo < p # lo is always less than p
114
115 2**64 + lo == 2**64 + (lo - p) (mod p) # result is congruent
116
117 = lo - p # exactly the same value as the previous RHS
118 # in uint64_t arithmetic.
119
120 p < 2**64 + lo < 2*p -> 0 < 2**64 + (lo - p) < p # correct range
121
122
123
124 [1] http://www.apfloat.org/apfloat/2.40/apfloat.pdf
125
126
127
OLDNEW
« no previous file with comments | « Modules/_decimal/libmpdec/literature/matrix-transform.txt ('k') | Modules/_decimal/libmpdec/literature/mulmod-ppro.txt » ('j') | no next file with comments »

RSS Feeds Recent Issues | This issue
This is Rietveld 894c83f36cb7+