CoCalc -- mulmod-64.txt

GitHub Repository: allendowney/cpython
Path: blob/main/Modules/_decimal/libmpdec/literature/mulmod-64.txt
¹² views
1

2

3
(* Copyright (c) 2011-2020 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 unnecessary 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

128

129
Product

Resources

Company
