| Austin Schuh | dace2a6 | 2020-08-18 10:56:48 -0700 | [diff] [blame] | 1 | /* invert.c -- Compute floor((B^{2n}-1)/U) - B^n. |
| 2 | |
| 3 | Contributed to the GNU project by Marco Bodrato. |
| 4 | |
| 5 | THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY |
| 6 | SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST |
| 7 | GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. |
| 8 | |
| 9 | Copyright (C) 2007, 2009, 2010, 2012, 2014-2016 Free Software Foundation, Inc. |
| 10 | |
| 11 | This file is part of the GNU MP Library. |
| 12 | |
| 13 | The GNU MP Library is free software; you can redistribute it and/or modify |
| 14 | it under the terms of either: |
| 15 | |
| 16 | * the GNU Lesser General Public License as published by the Free |
| 17 | Software Foundation; either version 3 of the License, or (at your |
| 18 | option) any later version. |
| 19 | |
| 20 | or |
| 21 | |
| 22 | * the GNU General Public License as published by the Free Software |
| 23 | Foundation; either version 2 of the License, or (at your option) any |
| 24 | later version. |
| 25 | |
| 26 | or both in parallel, as here. |
| 27 | |
| 28 | The GNU MP Library is distributed in the hope that it will be useful, but |
| 29 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| 30 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| 31 | for more details. |
| 32 | |
| 33 | You should have received copies of the GNU General Public License and the |
| 34 | GNU Lesser General Public License along with the GNU MP Library. If not, |
| 35 | see https://www.gnu.org/licenses/. */ |
| 36 | |
| 37 | #include "gmp-impl.h" |
| 38 | #include "longlong.h" |
| 39 | |
| 40 | void |
| 41 | mpn_invert (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch) |
| 42 | { |
| 43 | ASSERT (n > 0); |
| 44 | ASSERT (dp[n-1] & GMP_NUMB_HIGHBIT); |
| 45 | ASSERT (! MPN_OVERLAP_P (ip, n, dp, n)); |
| 46 | ASSERT (! MPN_OVERLAP_P (ip, n, scratch, mpn_invertappr_itch(n))); |
| 47 | ASSERT (! MPN_OVERLAP_P (dp, n, scratch, mpn_invertappr_itch(n))); |
| 48 | |
| 49 | if (n == 1) |
| 50 | invert_limb (*ip, *dp); |
| 51 | else if (BELOW_THRESHOLD (n, INV_APPR_THRESHOLD)) |
| 52 | { |
| 53 | /* Maximum scratch needed by this branch: 2*n */ |
| 54 | mp_ptr xp; |
| 55 | |
| 56 | xp = scratch; /* 2 * n limbs */ |
| 57 | /* n > 1 here */ |
| 58 | MPN_FILL (xp, n, GMP_NUMB_MAX); |
| 59 | mpn_com (xp + n, dp, n); |
| 60 | if (n == 2) { |
| 61 | mpn_divrem_2 (ip, 0, xp, 4, dp); |
| 62 | } else { |
| 63 | gmp_pi1_t inv; |
| 64 | invert_pi1 (inv, dp[n-1], dp[n-2]); |
| 65 | /* FIXME: should we use dcpi1_div_q, for big sizes? */ |
| 66 | mpn_sbpi1_div_q (ip, xp, 2 * n, dp, n, inv.inv32); |
| 67 | } |
| 68 | } |
| 69 | else { /* Use approximated inverse; correct the result if needed. */ |
| 70 | mp_limb_t e; /* The possible error in the approximate inverse */ |
| 71 | |
| 72 | ASSERT ( mpn_invert_itch (n) >= mpn_invertappr_itch (n) ); |
| 73 | e = mpn_ni_invertappr (ip, dp, n, scratch); |
| 74 | |
| 75 | if (UNLIKELY (e)) { /* Assume the error can only be "0" (no error) or "1". */ |
| 76 | /* Code to detect and correct the "off by one" approximation. */ |
| 77 | mpn_mul_n (scratch, ip, dp, n); |
| 78 | e = mpn_add_n (scratch, scratch, dp, n); /* FIXME: we only need e.*/ |
| 79 | if (LIKELY(e)) /* The high part can not give a carry by itself. */ |
| 80 | e = mpn_add_nc (scratch + n, scratch + n, dp, n, e); /* FIXME:e */ |
| 81 | /* If the value was wrong (no carry), correct it (increment). */ |
| 82 | e ^= CNST_LIMB (1); |
| 83 | MPN_INCR_U (ip, n, e); |
| 84 | } |
| 85 | } |
| 86 | } |