# Python code coverage for Modules/_decimal/libmpdec/sixstep.c

# | 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 | #include "mpdecimal.h" |

30 | n/a | #include <stdio.h> |

31 | n/a | #include <stdlib.h> |

32 | n/a | #include <assert.h> |

33 | n/a | #include "bits.h" |

34 | n/a | #include "difradix2.h" |

35 | n/a | #include "numbertheory.h" |

36 | n/a | #include "transpose.h" |

37 | n/a | #include "umodarith.h" |

38 | n/a | #include "sixstep.h" |

39 | n/a | |

40 | n/a | |

41 | n/a | /* Bignum: Cache efficient Matrix Fourier Transform for arrays of the |

42 | n/a | form 2**n (See literature/six-step.txt). */ |

43 | n/a | |

44 | n/a | |

45 | n/a | /* forward transform with sign = -1 */ |

46 | n/a | int |

47 | n/a | six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum) |

48 | n/a | { |

49 | n/a | struct fnt_params *tparams; |

50 | n/a | mpd_size_t log2n, C, R; |

51 | n/a | mpd_uint_t kernel; |

52 | n/a | mpd_uint_t umod; |

53 | n/a | #ifdef PPRO |

54 | n/a | double dmod; |

55 | n/a | uint32_t dinvmod[3]; |

56 | n/a | #endif |

57 | n/a | mpd_uint_t *x, w0, w1, wstep; |

58 | n/a | mpd_size_t i, k; |

59 | n/a | |

60 | n/a | |

61 | n/a | assert(ispower2(n)); |

62 | n/a | assert(n >= 16); |

63 | n/a | assert(n <= MPD_MAXTRANSFORM_2N); |

64 | n/a | |

65 | n/a | log2n = mpd_bsr(n); |

66 | n/a | C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */ |

67 | n/a | R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */ |

68 | n/a | |

69 | n/a | |

70 | n/a | /* Transpose the matrix. */ |

71 | n/a | if (!transpose_pow2(a, R, C)) { |

72 | n/a | return 0; |

73 | n/a | } |

74 | n/a | |

75 | n/a | /* Length R transform on the rows. */ |

76 | n/a | if ((tparams = _mpd_init_fnt_params(R, -1, modnum)) == NULL) { |

77 | n/a | return 0; |

78 | n/a | } |

79 | n/a | for (x = a; x < a+n; x += R) { |

80 | n/a | fnt_dif2(x, R, tparams); |

81 | n/a | } |

82 | n/a | |

83 | n/a | /* Transpose the matrix. */ |

84 | n/a | if (!transpose_pow2(a, C, R)) { |

85 | n/a | mpd_free(tparams); |

86 | n/a | return 0; |

87 | n/a | } |

88 | n/a | |

89 | n/a | /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */ |

90 | n/a | SETMODULUS(modnum); |

91 | n/a | kernel = _mpd_getkernel(n, -1, modnum); |

92 | n/a | for (i = 1; i < R; i++) { |

93 | n/a | w0 = 1; /* r**(i*0): initial value for k=0 */ |

94 | n/a | w1 = POWMOD(kernel, i); /* r**(i*1): initial value for k=1 */ |

95 | n/a | wstep = MULMOD(w1, w1); /* r**(2*i) */ |

96 | n/a | for (k = 0; k < C; k += 2) { |

97 | n/a | mpd_uint_t x0 = a[i*C+k]; |

98 | n/a | mpd_uint_t x1 = a[i*C+k+1]; |

99 | n/a | MULMOD2(&x0, w0, &x1, w1); |

100 | n/a | MULMOD2C(&w0, &w1, wstep); /* r**(i*(k+2)) = r**(i*k) * r**(2*i) */ |

101 | n/a | a[i*C+k] = x0; |

102 | n/a | a[i*C+k+1] = x1; |

103 | n/a | } |

104 | n/a | } |

105 | n/a | |

106 | n/a | /* Length C transform on the rows. */ |

107 | n/a | if (C != R) { |

108 | n/a | mpd_free(tparams); |

109 | n/a | if ((tparams = _mpd_init_fnt_params(C, -1, modnum)) == NULL) { |

110 | n/a | return 0; |

111 | n/a | } |

112 | n/a | } |

113 | n/a | for (x = a; x < a+n; x += C) { |

114 | n/a | fnt_dif2(x, C, tparams); |

115 | n/a | } |

116 | n/a | mpd_free(tparams); |

117 | n/a | |

118 | n/a | #if 0 |

119 | n/a | /* An unordered transform is sufficient for convolution. */ |

120 | n/a | /* Transpose the matrix. */ |

121 | n/a | if (!transpose_pow2(a, R, C)) { |

122 | n/a | return 0; |

123 | n/a | } |

124 | n/a | #endif |

125 | n/a | |

126 | n/a | return 1; |

127 | n/a | } |

128 | n/a | |

129 | n/a | |

130 | n/a | /* reverse transform, sign = 1 */ |

131 | n/a | int |

132 | n/a | inv_six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum) |

133 | n/a | { |

134 | n/a | struct fnt_params *tparams; |

135 | n/a | mpd_size_t log2n, C, R; |

136 | n/a | mpd_uint_t kernel; |

137 | n/a | mpd_uint_t umod; |

138 | n/a | #ifdef PPRO |

139 | n/a | double dmod; |

140 | n/a | uint32_t dinvmod[3]; |

141 | n/a | #endif |

142 | n/a | mpd_uint_t *x, w0, w1, wstep; |

143 | n/a | mpd_size_t i, k; |

144 | n/a | |

145 | n/a | |

146 | n/a | assert(ispower2(n)); |

147 | n/a | assert(n >= 16); |

148 | n/a | assert(n <= MPD_MAXTRANSFORM_2N); |

149 | n/a | |

150 | n/a | log2n = mpd_bsr(n); |

151 | n/a | C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */ |

152 | n/a | R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */ |

153 | n/a | |

154 | n/a | |

155 | n/a | #if 0 |

156 | n/a | /* An unordered transform is sufficient for convolution. */ |

157 | n/a | /* Transpose the matrix, producing an R*C matrix. */ |

158 | n/a | if (!transpose_pow2(a, C, R)) { |

159 | n/a | return 0; |

160 | n/a | } |

161 | n/a | #endif |

162 | n/a | |

163 | n/a | /* Length C transform on the rows. */ |

164 | n/a | if ((tparams = _mpd_init_fnt_params(C, 1, modnum)) == NULL) { |

165 | n/a | return 0; |

166 | n/a | } |

167 | n/a | for (x = a; x < a+n; x += C) { |

168 | n/a | fnt_dif2(x, C, tparams); |

169 | n/a | } |

170 | n/a | |

171 | n/a | /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */ |

172 | n/a | SETMODULUS(modnum); |

173 | n/a | kernel = _mpd_getkernel(n, 1, modnum); |

174 | n/a | for (i = 1; i < R; i++) { |

175 | n/a | w0 = 1; |

176 | n/a | w1 = POWMOD(kernel, i); |

177 | n/a | wstep = MULMOD(w1, w1); |

178 | n/a | for (k = 0; k < C; k += 2) { |

179 | n/a | mpd_uint_t x0 = a[i*C+k]; |

180 | n/a | mpd_uint_t x1 = a[i*C+k+1]; |

181 | n/a | MULMOD2(&x0, w0, &x1, w1); |

182 | n/a | MULMOD2C(&w0, &w1, wstep); |

183 | n/a | a[i*C+k] = x0; |

184 | n/a | a[i*C+k+1] = x1; |

185 | n/a | } |

186 | n/a | } |

187 | n/a | |

188 | n/a | /* Transpose the matrix. */ |

189 | n/a | if (!transpose_pow2(a, R, C)) { |

190 | n/a | mpd_free(tparams); |

191 | n/a | return 0; |

192 | n/a | } |

193 | n/a | |

194 | n/a | /* Length R transform on the rows. */ |

195 | n/a | if (R != C) { |

196 | n/a | mpd_free(tparams); |

197 | n/a | if ((tparams = _mpd_init_fnt_params(R, 1, modnum)) == NULL) { |

198 | n/a | return 0; |

199 | n/a | } |

200 | n/a | } |

201 | n/a | for (x = a; x < a+n; x += R) { |

202 | n/a | fnt_dif2(x, R, tparams); |

203 | n/a | } |

204 | n/a | mpd_free(tparams); |

205 | n/a | |

206 | n/a | /* Transpose the matrix. */ |

207 | n/a | if (!transpose_pow2(a, C, R)) { |

208 | n/a | return 0; |

209 | n/a | } |

210 | n/a | |

211 | n/a | return 1; |

212 | n/a | } |

213 | n/a | |

214 | n/a |