/src/gmp-6.2.1/mpn/sbpi1_bdiv_qr.c
Line | Count | Source |
1 | | /* mpn_sbpi1_bdiv_qr -- schoolbook Hensel division with precomputed inverse, |
2 | | returning quotient and remainder. |
3 | | |
4 | | Contributed to the GNU project by Niels Möller and Torbjörn Granlund. |
5 | | |
6 | | THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE INTERFACES. |
7 | | IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS |
8 | | ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. |
9 | | |
10 | | Copyright 2006, 2009, 2011, 2012, 2017 Free Software Foundation, Inc. |
11 | | |
12 | | This file is part of the GNU MP Library. |
13 | | |
14 | | The GNU MP Library is free software; you can redistribute it and/or modify |
15 | | it under the terms of either: |
16 | | |
17 | | * the GNU Lesser General Public License as published by the Free |
18 | | Software Foundation; either version 3 of the License, or (at your |
19 | | option) any later version. |
20 | | |
21 | | or |
22 | | |
23 | | * the GNU General Public License as published by the Free Software |
24 | | Foundation; either version 2 of the License, or (at your option) any |
25 | | later version. |
26 | | |
27 | | or both in parallel, as here. |
28 | | |
29 | | The GNU MP Library is distributed in the hope that it will be useful, but |
30 | | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
31 | | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
32 | | for more details. |
33 | | |
34 | | You should have received copies of the GNU General Public License and the |
35 | | GNU Lesser General Public License along with the GNU MP Library. If not, |
36 | | see https://www.gnu.org/licenses/. */ |
37 | | |
38 | | #include "gmp-impl.h" |
39 | | |
40 | | |
41 | | /* Computes a binary quotient of size qn = un - dn. |
42 | | Output: |
43 | | |
44 | | Q = -U * D^{-1} mod B^qn, |
45 | | |
46 | | R = (U + Q * D) * B^(-qn) |
47 | | |
48 | | Stores the dn least significant limbs of R at {up + un - dn, dn}, |
49 | | and returns the carry from the addition N + Q*D. |
50 | | |
51 | | D must be odd. dinv is (-D)^-1 mod B. */ |
52 | | |
53 | | mp_limb_t |
54 | | mpn_sbpi1_bdiv_qr (mp_ptr qp, |
55 | | mp_ptr up, mp_size_t un, |
56 | | mp_srcptr dp, mp_size_t dn, mp_limb_t dinv) |
57 | 831 | { |
58 | 831 | mp_size_t i; |
59 | 831 | mp_limb_t cy; |
60 | | |
61 | 831 | ASSERT (dn > 0); |
62 | 831 | ASSERT (un > dn); |
63 | 831 | ASSERT ((dp[0] & 1) != 0); |
64 | 831 | ASSERT (-(dp[0] * dinv) == 1); |
65 | 831 | ASSERT (up == qp || !MPN_OVERLAP_P (up, un, qp, un - dn)); |
66 | | |
67 | 7.57k | for (i = un - dn, cy = 0; i != 0; i--) |
68 | 6.74k | { |
69 | 6.74k | mp_limb_t q = dinv * up[0]; |
70 | 6.74k | mp_limb_t hi = mpn_addmul_1 (up, dp, dn, q); |
71 | 6.74k | *qp++ = q; |
72 | | |
73 | 6.74k | hi += cy; |
74 | 6.74k | cy = hi < cy; |
75 | 6.74k | hi += up[dn]; |
76 | 6.74k | cy += hi < up[dn]; |
77 | 6.74k | up[dn] = hi; |
78 | 6.74k | up++; |
79 | 6.74k | } |
80 | | |
81 | 831 | return cy; |
82 | 831 | } |