# Python code coverage for Modules/_decimal/libmpdec/basearith.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 <string.h> |

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

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

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

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

37 | n/a | |

38 | n/a | |

39 | n/a | /*********************************************************************/ |

40 | n/a | /* Calculations in base MPD_RADIX */ |

41 | n/a | /*********************************************************************/ |

42 | n/a | |

43 | n/a | |

44 | n/a | /* |

45 | n/a | * Knuth, TAOCP, Volume 2, 4.3.1: |

46 | n/a | * w := sum of u (len m) and v (len n) |

47 | n/a | * n > 0 and m >= n |

48 | n/a | * The calling function has to handle a possible final carry. |

49 | n/a | */ |

50 | n/a | mpd_uint_t |

51 | n/a | _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |

52 | n/a | mpd_size_t m, mpd_size_t n) |

53 | n/a | { |

54 | n/a | mpd_uint_t s; |

55 | n/a | mpd_uint_t carry = 0; |

56 | n/a | mpd_size_t i; |

57 | n/a | |

58 | n/a | assert(n > 0 && m >= n); |

59 | n/a | |

60 | n/a | /* add n members of u and v */ |

61 | n/a | for (i = 0; i < n; i++) { |

62 | n/a | s = u[i] + (v[i] + carry); |

63 | n/a | carry = (s < u[i]) | (s >= MPD_RADIX); |

64 | n/a | w[i] = carry ? s-MPD_RADIX : s; |

65 | n/a | } |

66 | n/a | /* if there is a carry, propagate it */ |

67 | n/a | for (; carry && i < m; i++) { |

68 | n/a | s = u[i] + carry; |

69 | n/a | carry = (s == MPD_RADIX); |

70 | n/a | w[i] = carry ? 0 : s; |

71 | n/a | } |

72 | n/a | /* copy the rest of u */ |

73 | n/a | for (; i < m; i++) { |

74 | n/a | w[i] = u[i]; |

75 | n/a | } |

76 | n/a | |

77 | n/a | return carry; |

78 | n/a | } |

79 | n/a | |

80 | n/a | /* |

81 | n/a | * Add the contents of u to w. Carries are propagated further. The caller |

82 | n/a | * has to make sure that w is big enough. |

83 | n/a | */ |

84 | n/a | void |

85 | n/a | _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |

86 | n/a | { |

87 | n/a | mpd_uint_t s; |

88 | n/a | mpd_uint_t carry = 0; |

89 | n/a | mpd_size_t i; |

90 | n/a | |

91 | n/a | if (n == 0) return; |

92 | n/a | |

93 | n/a | /* add n members of u to w */ |

94 | n/a | for (i = 0; i < n; i++) { |

95 | n/a | s = w[i] + (u[i] + carry); |

96 | n/a | carry = (s < w[i]) | (s >= MPD_RADIX); |

97 | n/a | w[i] = carry ? s-MPD_RADIX : s; |

98 | n/a | } |

99 | n/a | /* if there is a carry, propagate it */ |

100 | n/a | for (; carry; i++) { |

101 | n/a | s = w[i] + carry; |

102 | n/a | carry = (s == MPD_RADIX); |

103 | n/a | w[i] = carry ? 0 : s; |

104 | n/a | } |

105 | n/a | } |

106 | n/a | |

107 | n/a | /* |

108 | n/a | * Add v to w (len m). The calling function has to handle a possible |

109 | n/a | * final carry. Assumption: m > 0. |

110 | n/a | */ |

111 | n/a | mpd_uint_t |

112 | n/a | _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v) |

113 | n/a | { |

114 | n/a | mpd_uint_t s; |

115 | n/a | mpd_uint_t carry; |

116 | n/a | mpd_size_t i; |

117 | n/a | |

118 | n/a | assert(m > 0); |

119 | n/a | |

120 | n/a | /* add v to w */ |

121 | n/a | s = w[0] + v; |

122 | n/a | carry = (s < v) | (s >= MPD_RADIX); |

123 | n/a | w[0] = carry ? s-MPD_RADIX : s; |

124 | n/a | |

125 | n/a | /* if there is a carry, propagate it */ |

126 | n/a | for (i = 1; carry && i < m; i++) { |

127 | n/a | s = w[i] + carry; |

128 | n/a | carry = (s == MPD_RADIX); |

129 | n/a | w[i] = carry ? 0 : s; |

130 | n/a | } |

131 | n/a | |

132 | n/a | return carry; |

133 | n/a | } |

134 | n/a | |

135 | n/a | /* Increment u. The calling function has to handle a possible carry. */ |

136 | n/a | mpd_uint_t |

137 | n/a | _mpd_baseincr(mpd_uint_t *u, mpd_size_t n) |

138 | n/a | { |

139 | n/a | mpd_uint_t s; |

140 | n/a | mpd_uint_t carry = 1; |

141 | n/a | mpd_size_t i; |

142 | n/a | |

143 | n/a | assert(n > 0); |

144 | n/a | |

145 | n/a | /* if there is a carry, propagate it */ |

146 | n/a | for (i = 0; carry && i < n; i++) { |

147 | n/a | s = u[i] + carry; |

148 | n/a | carry = (s == MPD_RADIX); |

149 | n/a | u[i] = carry ? 0 : s; |

150 | n/a | } |

151 | n/a | |

152 | n/a | return carry; |

153 | n/a | } |

154 | n/a | |

155 | n/a | /* |

156 | n/a | * Knuth, TAOCP, Volume 2, 4.3.1: |

157 | n/a | * w := difference of u (len m) and v (len n). |

158 | n/a | * number in u >= number in v; |

159 | n/a | */ |

160 | n/a | void |

161 | n/a | _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |

162 | n/a | mpd_size_t m, mpd_size_t n) |

163 | n/a | { |

164 | n/a | mpd_uint_t d; |

165 | n/a | mpd_uint_t borrow = 0; |

166 | n/a | mpd_size_t i; |

167 | n/a | |

168 | n/a | assert(m > 0 && n > 0); |

169 | n/a | |

170 | n/a | /* subtract n members of v from u */ |

171 | n/a | for (i = 0; i < n; i++) { |

172 | n/a | d = u[i] - (v[i] + borrow); |

173 | n/a | borrow = (u[i] < d); |

174 | n/a | w[i] = borrow ? d + MPD_RADIX : d; |

175 | n/a | } |

176 | n/a | /* if there is a borrow, propagate it */ |

177 | n/a | for (; borrow && i < m; i++) { |

178 | n/a | d = u[i] - borrow; |

179 | n/a | borrow = (u[i] == 0); |

180 | n/a | w[i] = borrow ? MPD_RADIX-1 : d; |

181 | n/a | } |

182 | n/a | /* copy the rest of u */ |

183 | n/a | for (; i < m; i++) { |

184 | n/a | w[i] = u[i]; |

185 | n/a | } |

186 | n/a | } |

187 | n/a | |

188 | n/a | /* |

189 | n/a | * Subtract the contents of u from w. w is larger than u. Borrows are |

190 | n/a | * propagated further, but eventually w can absorb the final borrow. |

191 | n/a | */ |

192 | n/a | void |

193 | n/a | _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n) |

194 | n/a | { |

195 | n/a | mpd_uint_t d; |

196 | n/a | mpd_uint_t borrow = 0; |

197 | n/a | mpd_size_t i; |

198 | n/a | |

199 | n/a | if (n == 0) return; |

200 | n/a | |

201 | n/a | /* subtract n members of u from w */ |

202 | n/a | for (i = 0; i < n; i++) { |

203 | n/a | d = w[i] - (u[i] + borrow); |

204 | n/a | borrow = (w[i] < d); |

205 | n/a | w[i] = borrow ? d + MPD_RADIX : d; |

206 | n/a | } |

207 | n/a | /* if there is a borrow, propagate it */ |

208 | n/a | for (; borrow; i++) { |

209 | n/a | d = w[i] - borrow; |

210 | n/a | borrow = (w[i] == 0); |

211 | n/a | w[i] = borrow ? MPD_RADIX-1 : d; |

212 | n/a | } |

213 | n/a | } |

214 | n/a | |

215 | n/a | /* w := product of u (len n) and v (single word) */ |

216 | n/a | void |

217 | n/a | _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |

218 | n/a | { |

219 | n/a | mpd_uint_t hi, lo; |

220 | n/a | mpd_uint_t carry = 0; |

221 | n/a | mpd_size_t i; |

222 | n/a | |

223 | n/a | assert(n > 0); |

224 | n/a | |

225 | n/a | for (i=0; i < n; i++) { |

226 | n/a | |

227 | n/a | _mpd_mul_words(&hi, &lo, u[i], v); |

228 | n/a | lo = carry + lo; |

229 | n/a | if (lo < carry) hi++; |

230 | n/a | |

231 | n/a | _mpd_div_words_r(&carry, &w[i], hi, lo); |

232 | n/a | } |

233 | n/a | w[i] = carry; |

234 | n/a | } |

235 | n/a | |

236 | n/a | /* |

237 | n/a | * Knuth, TAOCP, Volume 2, 4.3.1: |

238 | n/a | * w := product of u (len m) and v (len n) |

239 | n/a | * w must be initialized to zero |

240 | n/a | */ |

241 | n/a | void |

242 | n/a | _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |

243 | n/a | mpd_size_t m, mpd_size_t n) |

244 | n/a | { |

245 | n/a | mpd_uint_t hi, lo; |

246 | n/a | mpd_uint_t carry; |

247 | n/a | mpd_size_t i, j; |

248 | n/a | |

249 | n/a | assert(m > 0 && n > 0); |

250 | n/a | |

251 | n/a | for (j=0; j < n; j++) { |

252 | n/a | carry = 0; |

253 | n/a | for (i=0; i < m; i++) { |

254 | n/a | |

255 | n/a | _mpd_mul_words(&hi, &lo, u[i], v[j]); |

256 | n/a | lo = w[i+j] + lo; |

257 | n/a | if (lo < w[i+j]) hi++; |

258 | n/a | lo = carry + lo; |

259 | n/a | if (lo < carry) hi++; |

260 | n/a | |

261 | n/a | _mpd_div_words_r(&carry, &w[i+j], hi, lo); |

262 | n/a | } |

263 | n/a | w[j+m] = carry; |

264 | n/a | } |

265 | n/a | } |

266 | n/a | |

267 | n/a | /* |

268 | n/a | * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |

269 | n/a | * w := quotient of u (len n) divided by a single word v |

270 | n/a | */ |

271 | n/a | mpd_uint_t |

272 | n/a | _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |

273 | n/a | { |

274 | n/a | mpd_uint_t hi, lo; |

275 | n/a | mpd_uint_t rem = 0; |

276 | n/a | mpd_size_t i; |

277 | n/a | |

278 | n/a | assert(n > 0); |

279 | n/a | |

280 | n/a | for (i=n-1; i != MPD_SIZE_MAX; i--) { |

281 | n/a | |

282 | n/a | _mpd_mul_words(&hi, &lo, rem, MPD_RADIX); |

283 | n/a | lo = u[i] + lo; |

284 | n/a | if (lo < u[i]) hi++; |

285 | n/a | |

286 | n/a | _mpd_div_words(&w[i], &rem, hi, lo, v); |

287 | n/a | } |

288 | n/a | |

289 | n/a | return rem; |

290 | n/a | } |

291 | n/a | |

292 | n/a | /* |

293 | n/a | * Knuth, TAOCP Volume 2, 4.3.1: |

294 | n/a | * q, r := quotient and remainder of uconst (len nplusm) |

295 | n/a | * divided by vconst (len n) |

296 | n/a | * nplusm >= n |

297 | n/a | * |

298 | n/a | * If r is not NULL, r will contain the remainder. If r is NULL, the |

299 | n/a | * return value indicates if there is a remainder: 1 for true, 0 for |

300 | n/a | * false. A return value of -1 indicates an error. |

301 | n/a | */ |

302 | n/a | int |

303 | n/a | _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r, |

304 | n/a | const mpd_uint_t *uconst, const mpd_uint_t *vconst, |

305 | n/a | mpd_size_t nplusm, mpd_size_t n) |

306 | n/a | { |

307 | n/a | mpd_uint_t ustatic[MPD_MINALLOC_MAX]; |

308 | n/a | mpd_uint_t vstatic[MPD_MINALLOC_MAX]; |

309 | n/a | mpd_uint_t *u = ustatic; |

310 | n/a | mpd_uint_t *v = vstatic; |

311 | n/a | mpd_uint_t d, qhat, rhat, w2[2]; |

312 | n/a | mpd_uint_t hi, lo, x; |

313 | n/a | mpd_uint_t carry; |

314 | n/a | mpd_size_t i, j, m; |

315 | n/a | int retval = 0; |

316 | n/a | |

317 | n/a | assert(n > 1 && nplusm >= n); |

318 | n/a | m = sub_size_t(nplusm, n); |

319 | n/a | |

320 | n/a | /* D1: normalize */ |

321 | n/a | d = MPD_RADIX / (vconst[n-1] + 1); |

322 | n/a | |

323 | n/a | if (nplusm >= MPD_MINALLOC_MAX) { |

324 | n/a | if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) { |

325 | n/a | return -1; |

326 | n/a | } |

327 | n/a | } |

328 | n/a | if (n >= MPD_MINALLOC_MAX) { |

329 | n/a | if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) { |

330 | n/a | mpd_free(u); |

331 | n/a | return -1; |

332 | n/a | } |

333 | n/a | } |

334 | n/a | |

335 | n/a | _mpd_shortmul(u, uconst, nplusm, d); |

336 | n/a | _mpd_shortmul(v, vconst, n, d); |

337 | n/a | |

338 | n/a | /* D2: loop */ |

339 | n/a | for (j=m; j != MPD_SIZE_MAX; j--) { |

340 | n/a | |

341 | n/a | /* D3: calculate qhat and rhat */ |

342 | n/a | rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]); |

343 | n/a | qhat = w2[1] * MPD_RADIX + w2[0]; |

344 | n/a | |

345 | n/a | while (1) { |

346 | n/a | if (qhat < MPD_RADIX) { |

347 | n/a | _mpd_singlemul(w2, qhat, v[n-2]); |

348 | n/a | if (w2[1] <= rhat) { |

349 | n/a | if (w2[1] != rhat || w2[0] <= u[j+n-2]) { |

350 | n/a | break; |

351 | n/a | } |

352 | n/a | } |

353 | n/a | } |

354 | n/a | qhat -= 1; |

355 | n/a | rhat += v[n-1]; |

356 | n/a | if (rhat < v[n-1] || rhat >= MPD_RADIX) { |

357 | n/a | break; |

358 | n/a | } |

359 | n/a | } |

360 | n/a | /* D4: multiply and subtract */ |

361 | n/a | carry = 0; |

362 | n/a | for (i=0; i <= n; i++) { |

363 | n/a | |

364 | n/a | _mpd_mul_words(&hi, &lo, qhat, v[i]); |

365 | n/a | |

366 | n/a | lo = carry + lo; |

367 | n/a | if (lo < carry) hi++; |

368 | n/a | |

369 | n/a | _mpd_div_words_r(&hi, &lo, hi, lo); |

370 | n/a | |

371 | n/a | x = u[i+j] - lo; |

372 | n/a | carry = (u[i+j] < x); |

373 | n/a | u[i+j] = carry ? x+MPD_RADIX : x; |

374 | n/a | carry += hi; |

375 | n/a | } |

376 | n/a | q[j] = qhat; |

377 | n/a | /* D5: test remainder */ |

378 | n/a | if (carry) { |

379 | n/a | q[j] -= 1; |

380 | n/a | /* D6: add back */ |

381 | n/a | (void)_mpd_baseadd(u+j, u+j, v, n+1, n); |

382 | n/a | } |

383 | n/a | } |

384 | n/a | |

385 | n/a | /* D8: unnormalize */ |

386 | n/a | if (r != NULL) { |

387 | n/a | _mpd_shortdiv(r, u, n, d); |

388 | n/a | /* we are not interested in the return value here */ |

389 | n/a | retval = 0; |

390 | n/a | } |

391 | n/a | else { |

392 | n/a | retval = !_mpd_isallzero(u, n); |

393 | n/a | } |

394 | n/a | |

395 | n/a | |

396 | n/a | if (u != ustatic) mpd_free(u); |

397 | n/a | if (v != vstatic) mpd_free(v); |

398 | n/a | return retval; |

399 | n/a | } |

400 | n/a | |

401 | n/a | /* |

402 | n/a | * Left shift of src by 'shift' digits; src may equal dest. |

403 | n/a | * |

404 | n/a | * dest := area of n mpd_uint_t with space for srcdigits+shift digits. |

405 | n/a | * src := coefficient with length m. |

406 | n/a | * |

407 | n/a | * The case splits in the function are non-obvious. The following |

408 | n/a | * equations might help: |

409 | n/a | * |

410 | n/a | * Let msdigits denote the number of digits in the most significant |

411 | n/a | * word of src. Then 1 <= msdigits <= rdigits. |

412 | n/a | * |

413 | n/a | * 1) shift = q * rdigits + r |

414 | n/a | * 2) srcdigits = qsrc * rdigits + msdigits |

415 | n/a | * 3) destdigits = shift + srcdigits |

416 | n/a | * = q * rdigits + r + qsrc * rdigits + msdigits |

417 | n/a | * = q * rdigits + (qsrc * rdigits + (r + msdigits)) |

418 | n/a | * |

419 | n/a | * The result has q zero words, followed by the coefficient that |

420 | n/a | * is left-shifted by r. The case r == 0 is trivial. For r > 0, it |

421 | n/a | * is important to keep in mind that we always read m source words, |

422 | n/a | * but write m+1 destination words if r + msdigits > rdigits, m words |

423 | n/a | * otherwise. |

424 | n/a | */ |

425 | n/a | void |

426 | n/a | _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m, |

427 | n/a | mpd_size_t shift) |

428 | n/a | { |

429 | n/a | #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |

430 | n/a | /* spurious uninitialized warnings */ |

431 | n/a | mpd_uint_t l=l, lprev=lprev, h=h; |

432 | n/a | #else |

433 | n/a | mpd_uint_t l, lprev, h; |

434 | n/a | #endif |

435 | n/a | mpd_uint_t q, r; |

436 | n/a | mpd_uint_t ph; |

437 | n/a | |

438 | n/a | assert(m > 0 && n >= m); |

439 | n/a | |

440 | n/a | _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |

441 | n/a | |

442 | n/a | if (r != 0) { |

443 | n/a | |

444 | n/a | ph = mpd_pow10[r]; |

445 | n/a | |

446 | n/a | --m; --n; |

447 | n/a | _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r); |

448 | n/a | if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */ |

449 | n/a | dest[n--] = h; |

450 | n/a | } |

451 | n/a | /* write m-1 shifted words */ |

452 | n/a | for (; m != MPD_SIZE_MAX; m--,n--) { |

453 | n/a | _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r); |

454 | n/a | dest[n] = ph * lprev + h; |

455 | n/a | lprev = l; |

456 | n/a | } |

457 | n/a | /* write least significant word */ |

458 | n/a | dest[q] = ph * lprev; |

459 | n/a | } |

460 | n/a | else { |

461 | n/a | while (--m != MPD_SIZE_MAX) { |

462 | n/a | dest[m+q] = src[m]; |

463 | n/a | } |

464 | n/a | } |

465 | n/a | |

466 | n/a | mpd_uint_zero(dest, q); |

467 | n/a | } |

468 | n/a | |

469 | n/a | /* |

470 | n/a | * Right shift of src by 'shift' digits; src may equal dest. |

471 | n/a | * Assumption: srcdigits-shift > 0. |

472 | n/a | * |

473 | n/a | * dest := area with space for srcdigits-shift digits. |

474 | n/a | * src := coefficient with length 'slen'. |

475 | n/a | * |

476 | n/a | * The case splits in the function rely on the following equations: |

477 | n/a | * |

478 | n/a | * Let msdigits denote the number of digits in the most significant |

479 | n/a | * word of src. Then 1 <= msdigits <= rdigits. |

480 | n/a | * |

481 | n/a | * 1) shift = q * rdigits + r |

482 | n/a | * 2) srcdigits = qsrc * rdigits + msdigits |

483 | n/a | * 3) destdigits = srcdigits - shift |

484 | n/a | * = qsrc * rdigits + msdigits - (q * rdigits + r) |

485 | n/a | * = (qsrc - q) * rdigits + msdigits - r |

486 | n/a | * |

487 | n/a | * Since destdigits > 0 and 1 <= msdigits <= rdigits: |

488 | n/a | * |

489 | n/a | * 4) qsrc >= q |

490 | n/a | * 5) qsrc == q ==> msdigits > r |

491 | n/a | * |

492 | n/a | * The result has slen-q words if msdigits > r, slen-q-1 words otherwise. |

493 | n/a | */ |

494 | n/a | mpd_uint_t |

495 | n/a | _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen, |

496 | n/a | mpd_size_t shift) |

497 | n/a | { |

498 | n/a | #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__) |

499 | n/a | /* spurious uninitialized warnings */ |

500 | n/a | mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */ |

501 | n/a | #else |

502 | n/a | mpd_uint_t l, h, hprev; /* low, high, previous high */ |

503 | n/a | #endif |

504 | n/a | mpd_uint_t rnd, rest; /* rounding digit, rest */ |

505 | n/a | mpd_uint_t q, r; |

506 | n/a | mpd_size_t i, j; |

507 | n/a | mpd_uint_t ph; |

508 | n/a | |

509 | n/a | assert(slen > 0); |

510 | n/a | |

511 | n/a | _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS); |

512 | n/a | |

513 | n/a | rnd = rest = 0; |

514 | n/a | if (r != 0) { |

515 | n/a | |

516 | n/a | ph = mpd_pow10[MPD_RDIGITS-r]; |

517 | n/a | |

518 | n/a | _mpd_divmod_pow10(&hprev, &rest, src[q], r); |

519 | n/a | _mpd_divmod_pow10(&rnd, &rest, rest, r-1); |

520 | n/a | |

521 | n/a | if (rest == 0 && q > 0) { |

522 | n/a | rest = !_mpd_isallzero(src, q); |

523 | n/a | } |

524 | n/a | /* write slen-q-1 words */ |

525 | n/a | for (j=0,i=q+1; i<slen; i++,j++) { |

526 | n/a | _mpd_divmod_pow10(&h, &l, src[i], r); |

527 | n/a | dest[j] = ph * l + hprev; |

528 | n/a | hprev = h; |

529 | n/a | } |

530 | n/a | /* write most significant word */ |

531 | n/a | if (hprev != 0) { /* always the case if slen==q-1 */ |

532 | n/a | dest[j] = hprev; |

533 | n/a | } |

534 | n/a | } |

535 | n/a | else { |

536 | n/a | if (q > 0) { |

537 | n/a | _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1); |

538 | n/a | /* is there any non-zero digit below rnd? */ |

539 | n/a | if (rest == 0) rest = !_mpd_isallzero(src, q-1); |

540 | n/a | } |

541 | n/a | for (j = 0; j < slen-q; j++) { |

542 | n/a | dest[j] = src[q+j]; |

543 | n/a | } |

544 | n/a | } |

545 | n/a | |

546 | n/a | /* 0-4 ==> rnd+rest < 0.5 */ |

547 | n/a | /* 5 ==> rnd+rest == 0.5 */ |

548 | n/a | /* 6-9 ==> rnd+rest > 0.5 */ |

549 | n/a | return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd; |

550 | n/a | } |

551 | n/a | |

552 | n/a | |

553 | n/a | /*********************************************************************/ |

554 | n/a | /* Calculations in base b */ |

555 | n/a | /*********************************************************************/ |

556 | n/a | |

557 | n/a | /* |

558 | n/a | * Add v to w (len m). The calling function has to handle a possible |

559 | n/a | * final carry. Assumption: m > 0. |

560 | n/a | */ |

561 | n/a | mpd_uint_t |

562 | n/a | _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b) |

563 | n/a | { |

564 | n/a | mpd_uint_t s; |

565 | n/a | mpd_uint_t carry; |

566 | n/a | mpd_size_t i; |

567 | n/a | |

568 | n/a | assert(m > 0); |

569 | n/a | |

570 | n/a | /* add v to w */ |

571 | n/a | s = w[0] + v; |

572 | n/a | carry = (s < v) | (s >= b); |

573 | n/a | w[0] = carry ? s-b : s; |

574 | n/a | |

575 | n/a | /* if there is a carry, propagate it */ |

576 | n/a | for (i = 1; carry && i < m; i++) { |

577 | n/a | s = w[i] + carry; |

578 | n/a | carry = (s == b); |

579 | n/a | w[i] = carry ? 0 : s; |

580 | n/a | } |

581 | n/a | |

582 | n/a | return carry; |

583 | n/a | } |

584 | n/a | |

585 | n/a | /* w := product of u (len n) and v (single word). Return carry. */ |

586 | n/a | mpd_uint_t |

587 | n/a | _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v) |

588 | n/a | { |

589 | n/a | mpd_uint_t hi, lo; |

590 | n/a | mpd_uint_t carry = 0; |

591 | n/a | mpd_size_t i; |

592 | n/a | |

593 | n/a | assert(n > 0); |

594 | n/a | |

595 | n/a | for (i=0; i < n; i++) { |

596 | n/a | |

597 | n/a | _mpd_mul_words(&hi, &lo, u[i], v); |

598 | n/a | lo = carry + lo; |

599 | n/a | if (lo < carry) hi++; |

600 | n/a | |

601 | n/a | _mpd_div_words_r(&carry, &w[i], hi, lo); |

602 | n/a | } |

603 | n/a | |

604 | n/a | return carry; |

605 | n/a | } |

606 | n/a | |

607 | n/a | /* w := product of u (len n) and v (single word) */ |

608 | n/a | mpd_uint_t |

609 | n/a | _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |

610 | n/a | mpd_uint_t v, mpd_uint_t b) |

611 | n/a | { |

612 | n/a | mpd_uint_t hi, lo; |

613 | n/a | mpd_uint_t carry = 0; |

614 | n/a | mpd_size_t i; |

615 | n/a | |

616 | n/a | assert(n > 0); |

617 | n/a | |

618 | n/a | for (i=0; i < n; i++) { |

619 | n/a | |

620 | n/a | _mpd_mul_words(&hi, &lo, u[i], v); |

621 | n/a | lo = carry + lo; |

622 | n/a | if (lo < carry) hi++; |

623 | n/a | |

624 | n/a | _mpd_div_words(&carry, &w[i], hi, lo, b); |

625 | n/a | } |

626 | n/a | |

627 | n/a | return carry; |

628 | n/a | } |

629 | n/a | |

630 | n/a | /* |

631 | n/a | * Knuth, TAOCP Volume 2, 4.3.1, exercise 16: |

632 | n/a | * w := quotient of u (len n) divided by a single word v |

633 | n/a | */ |

634 | n/a | mpd_uint_t |

635 | n/a | _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |

636 | n/a | mpd_uint_t v, mpd_uint_t b) |

637 | n/a | { |

638 | n/a | mpd_uint_t hi, lo; |

639 | n/a | mpd_uint_t rem = 0; |

640 | n/a | mpd_size_t i; |

641 | n/a | |

642 | n/a | assert(n > 0); |

643 | n/a | |

644 | n/a | for (i=n-1; i != MPD_SIZE_MAX; i--) { |

645 | n/a | |

646 | n/a | _mpd_mul_words(&hi, &lo, rem, b); |

647 | n/a | lo = u[i] + lo; |

648 | n/a | if (lo < u[i]) hi++; |

649 | n/a | |

650 | n/a | _mpd_div_words(&w[i], &rem, hi, lo, v); |

651 | n/a | } |

652 | n/a | |

653 | n/a | return rem; |

654 | n/a | } |

655 | n/a | |

656 | n/a | |

657 | n/a |