# Python code coverage for Modules/_decimal/libmpdec/literature/fnt.py

# | count | content |
---|---|---|

1 | n/a | # |

2 | n/a | # Copyright (c) 2008-2016 Stefan Krah. All rights reserved. |

3 | n/a | # |

4 | n/a | # Redistribution and use in source and binary forms, with or without |

5 | n/a | # modification, are permitted provided that the following conditions |

6 | n/a | # are met: |

7 | n/a | # |

8 | n/a | # 1. Redistributions of source code must retain the above copyright |

9 | n/a | # notice, this list of conditions and the following disclaimer. |

10 | n/a | # |

11 | n/a | # 2. Redistributions in binary form must reproduce the above copyright |

12 | n/a | # notice, this list of conditions and the following disclaimer in the |

13 | n/a | # documentation and/or other materials provided with the distribution. |

14 | n/a | # |

15 | n/a | # THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |

16 | n/a | # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |

17 | n/a | # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |

18 | n/a | # ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |

19 | n/a | # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |

20 | n/a | # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |

21 | n/a | # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |

22 | n/a | # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |

23 | n/a | # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |

24 | n/a | # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |

25 | n/a | # SUCH DAMAGE. |

26 | n/a | # |

27 | n/a | |

28 | n/a | |

29 | n/a | ###################################################################### |

30 | n/a | # This file lists and checks some of the constants and limits used # |

31 | n/a | # in libmpdec's Number Theoretic Transform. At the end of the file # |

32 | n/a | # there is an example function for the plain DFT transform. # |

33 | n/a | ###################################################################### |

34 | n/a | |

35 | n/a | |

36 | n/a | # |

37 | n/a | # Number theoretic transforms are done in subfields of F(p). P[i] |

38 | n/a | # are the primes, D[i] = P[i] - 1 are highly composite and w[i] |

39 | n/a | # are the respective primitive roots of F(p). |

40 | n/a | # |

41 | n/a | # The strategy is to convolute two coefficients modulo all three |

42 | n/a | # primes, then use the Chinese Remainder Theorem on the three |

43 | n/a | # result arrays to recover the result in the usual base RADIX |

44 | n/a | # form. |

45 | n/a | # |

46 | n/a | |

47 | n/a | # ====================================================================== |

48 | n/a | # Primitive roots |

49 | n/a | # ====================================================================== |

50 | n/a | |

51 | n/a | # |

52 | n/a | # Verify primitive roots: |

53 | n/a | # |

54 | n/a | # For a prime field, r is a primitive root if and only if for all prime |

55 | n/a | # factors f of p-1, r**((p-1)/f) =/= 1 (mod p). |

56 | n/a | # |

57 | n/a | def prod(F, E): |

58 | n/a | """Check that the factorization of P-1 is correct. F is the list of |

59 | n/a | factors of P-1, E lists the number of occurrences of each factor.""" |

60 | n/a | x = 1 |

61 | n/a | for y, z in zip(F, E): |

62 | n/a | x *= y**z |

63 | n/a | return x |

64 | n/a | |

65 | n/a | def is_primitive_root(r, p, factors, exponents): |

66 | n/a | """Check if r is a primitive root of F(p).""" |

67 | n/a | if p != prod(factors, exponents) + 1: |

68 | n/a | return False |

69 | n/a | for f in factors: |

70 | n/a | q, control = divmod(p-1, f) |

71 | n/a | if control != 0: |

72 | n/a | return False |

73 | n/a | if pow(r, q, p) == 1: |

74 | n/a | return False |

75 | n/a | return True |

76 | n/a | |

77 | n/a | |

78 | n/a | # ================================================================= |

79 | n/a | # Constants and limits for the 64-bit version |

80 | n/a | # ================================================================= |

81 | n/a | |

82 | n/a | RADIX = 10**19 |

83 | n/a | |

84 | n/a | # Primes P1, P2 and P3: |

85 | n/a | P = [2**64-2**32+1, 2**64-2**34+1, 2**64-2**40+1] |

86 | n/a | |

87 | n/a | # P-1, highly composite. The transform length d is variable and |

88 | n/a | # must divide D = P-1. Since all D are divisible by 3 * 2**32, |

89 | n/a | # transform lengths can be 2**n or 3 * 2**n (where n <= 32). |

90 | n/a | D = [2**32 * 3 * (5 * 17 * 257 * 65537), |

91 | n/a | 2**34 * 3**2 * (7 * 11 * 31 * 151 * 331), |

92 | n/a | 2**40 * 3**2 * (5 * 7 * 13 * 17 * 241)] |

93 | n/a | |

94 | n/a | # Prime factors of P-1 and their exponents: |

95 | n/a | F = [(2,3,5,17,257,65537), (2,3,7,11,31,151,331), (2,3,5,7,13,17,241)] |

96 | n/a | E = [(32,1,1,1,1,1), (34,2,1,1,1,1,1), (40,2,1,1,1,1,1)] |

97 | n/a | |

98 | n/a | # Maximum transform length for 2**n. Above that only 3 * 2**31 |

99 | n/a | # or 3 * 2**32 are possible. |

100 | n/a | MPD_MAXTRANSFORM_2N = 2**32 |

101 | n/a | |

102 | n/a | |

103 | n/a | # Limits in the terminology of Pollard's paper: |

104 | n/a | m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. |

105 | n/a | M1 = M2 = RADIX-1 # Maximum value per single word. |

106 | n/a | L = m2 * M1 * M2 |

107 | n/a | P[0] * P[1] * P[2] > 2 * L |

108 | n/a | |

109 | n/a | |

110 | n/a | # Primitive roots of F(P1), F(P2) and F(P3): |

111 | n/a | w = [7, 10, 19] |

112 | n/a | |

113 | n/a | # The primitive roots are correct: |

114 | n/a | for i in range(3): |

115 | n/a | if not is_primitive_root(w[i], P[i], F[i], E[i]): |

116 | n/a | print("FAIL") |

117 | n/a | |

118 | n/a | |

119 | n/a | # ================================================================= |

120 | n/a | # Constants and limits for the 32-bit version |

121 | n/a | # ================================================================= |

122 | n/a | |

123 | n/a | RADIX = 10**9 |

124 | n/a | |

125 | n/a | # Primes P1, P2 and P3: |

126 | n/a | P = [2113929217, 2013265921, 1811939329] |

127 | n/a | |

128 | n/a | # P-1, highly composite. All D = P-1 are divisible by 3 * 2**25, |

129 | n/a | # allowing for transform lengths up to 3 * 2**25 words. |

130 | n/a | D = [2**25 * 3**2 * 7, |

131 | n/a | 2**27 * 3 * 5, |

132 | n/a | 2**26 * 3**3] |

133 | n/a | |

134 | n/a | # Prime factors of P-1 and their exponents: |

135 | n/a | F = [(2,3,7), (2,3,5), (2,3)] |

136 | n/a | E = [(25,2,1), (27,1,1), (26,3)] |

137 | n/a | |

138 | n/a | # Maximum transform length for 2**n. Above that only 3 * 2**24 or |

139 | n/a | # 3 * 2**25 are possible. |

140 | n/a | MPD_MAXTRANSFORM_2N = 2**25 |

141 | n/a | |

142 | n/a | |

143 | n/a | # Limits in the terminology of Pollard's paper: |

144 | n/a | m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. |

145 | n/a | M1 = M2 = RADIX-1 # Maximum value per single word. |

146 | n/a | L = m2 * M1 * M2 |

147 | n/a | P[0] * P[1] * P[2] > 2 * L |

148 | n/a | |

149 | n/a | |

150 | n/a | # Primitive roots of F(P1), F(P2) and F(P3): |

151 | n/a | w = [5, 31, 13] |

152 | n/a | |

153 | n/a | # The primitive roots are correct: |

154 | n/a | for i in range(3): |

155 | n/a | if not is_primitive_root(w[i], P[i], F[i], E[i]): |

156 | n/a | print("FAIL") |

157 | n/a | |

158 | n/a | |

159 | n/a | # ====================================================================== |

160 | n/a | # Example transform using a single prime |

161 | n/a | # ====================================================================== |

162 | n/a | |

163 | n/a | def ntt(lst, dir): |

164 | n/a | """Perform a transform on the elements of lst. len(lst) must |

165 | n/a | be 2**n or 3 * 2**n, where n <= 25. This is the slow DFT.""" |

166 | n/a | p = 2113929217 # prime |

167 | n/a | d = len(lst) # transform length |

168 | n/a | d_prime = pow(d, (p-2), p) # inverse of d |

169 | n/a | xi = (p-1)//d |

170 | n/a | w = 5 # primitive root of F(p) |

171 | n/a | r = pow(w, xi, p) # primitive root of the subfield |

172 | n/a | r_prime = pow(w, (p-1-xi), p) # inverse of r |

173 | n/a | if dir == 1: # forward transform |

174 | n/a | a = lst # input array |

175 | n/a | A = [0] * d # transformed values |

176 | n/a | for i in range(d): |

177 | n/a | s = 0 |

178 | n/a | for j in range(d): |

179 | n/a | s += a[j] * pow(r, i*j, p) |

180 | n/a | A[i] = s % p |

181 | n/a | return A |

182 | n/a | elif dir == -1: # backward transform |

183 | n/a | A = lst # input array |

184 | n/a | a = [0] * d # transformed values |

185 | n/a | for j in range(d): |

186 | n/a | s = 0 |

187 | n/a | for i in range(d): |

188 | n/a | s += A[i] * pow(r_prime, i*j, p) |

189 | n/a | a[j] = (d_prime * s) % p |

190 | n/a | return a |

191 | n/a | |

192 | n/a | def ntt_convolute(a, b): |

193 | n/a | """convolute arrays a and b.""" |

194 | n/a | assert(len(a) == len(b)) |

195 | n/a | x = ntt(a, 1) |

196 | n/a | y = ntt(b, 1) |

197 | n/a | for i in range(len(a)): |

198 | n/a | y[i] = y[i] * x[i] |

199 | n/a | r = ntt(y, -1) |

200 | n/a | return r |

201 | n/a | |

202 | n/a | |

203 | n/a | # Example: Two arrays representing 21 and 81 in little-endian: |

204 | n/a | a = [1, 2, 0, 0] |

205 | n/a | b = [1, 8, 0, 0] |

206 | n/a | |

207 | n/a | assert(ntt_convolute(a, b) == [1, 10, 16, 0]) |

208 | n/a | assert(21 * 81 == (1*10**0 + 10*10**1 + 16*10**2 + 0*10**3)) |