1 1.1 mrg /* mpn_add_err3_n -- add_n with three error terms 2 1.1 mrg 3 1.1 mrg Contributed by David Harvey. 4 1.1 mrg 5 1.1 mrg THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY 6 1.1 mrg SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 7 1.1 mrg GUARANTEED THAT IT'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. 8 1.1 mrg 9 1.1 mrg Copyright 2011 Free Software Foundation, Inc. 10 1.1 mrg 11 1.1 mrg This file is part of the GNU MP Library. 12 1.1 mrg 13 1.1 mrg The GNU MP Library is free software; you can redistribute it and/or modify 14 1.1.1.2 mrg it under the terms of either: 15 1.1.1.2 mrg 16 1.1.1.2 mrg * the GNU Lesser General Public License as published by the Free 17 1.1.1.2 mrg Software Foundation; either version 3 of the License, or (at your 18 1.1.1.2 mrg option) any later version. 19 1.1.1.2 mrg 20 1.1.1.2 mrg or 21 1.1.1.2 mrg 22 1.1.1.2 mrg * the GNU General Public License as published by the Free Software 23 1.1.1.2 mrg Foundation; either version 2 of the License, or (at your option) any 24 1.1.1.2 mrg later version. 25 1.1.1.2 mrg 26 1.1.1.2 mrg or both in parallel, as here. 27 1.1 mrg 28 1.1 mrg The GNU MP Library is distributed in the hope that it will be useful, but 29 1.1 mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 30 1.1.1.2 mrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 31 1.1.1.2 mrg for more details. 32 1.1 mrg 33 1.1.1.2 mrg You should have received copies of the GNU General Public License and the 34 1.1.1.2 mrg GNU Lesser General Public License along with the GNU MP Library. If not, 35 1.1.1.2 mrg see https://www.gnu.org/licenses/. */ 36 1.1 mrg 37 1.1 mrg #include "gmp-impl.h" 38 1.1 mrg 39 1.1 mrg /* 40 1.1 mrg Computes: 41 1.1 mrg 42 1.1 mrg (1) {rp,n} := {up,n} + {vp,n} (just like mpn_add_n) with incoming carry cy, 43 1.1 mrg return value is carry out. 44 1.1 mrg 45 1.1 mrg (2) Let c[i+1] = carry from i-th limb addition (c[0] = cy). 46 1.1 mrg Computes c[1]*yp1[n-1] + ... + c[n]*yp1[0], 47 1.1 mrg c[1]*yp2[n-1] + ... + c[n]*yp2[0], 48 1.1 mrg c[1]*yp3[n-1] + ... + c[n]*yp3[0], 49 1.1 mrg stores two-limb results at {ep,2}, {ep+2,2} and {ep+4,2} respectively. 50 1.1 mrg 51 1.1 mrg Requires n >= 1. 52 1.1 mrg 53 1.1 mrg None of the outputs may overlap each other or any of the inputs, except 54 1.1 mrg that {rp,n} may be equal to {up,n} or {vp,n}. 55 1.1 mrg */ 56 1.1 mrg mp_limb_t 57 1.1 mrg mpn_add_err3_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, 58 1.1 mrg mp_ptr ep, mp_srcptr yp1, mp_srcptr yp2, mp_srcptr yp3, 59 1.1 mrg mp_size_t n, mp_limb_t cy) 60 1.1 mrg { 61 1.1 mrg mp_limb_t el1, eh1, el2, eh2, el3, eh3, ul, vl, yl1, yl2, yl3, zl1, zl2, zl3, rl, sl, cy1, cy2; 62 1.1 mrg 63 1.1 mrg ASSERT (n >= 1); 64 1.1 mrg ASSERT (MPN_SAME_OR_SEPARATE_P (rp, up, n)); 65 1.1 mrg ASSERT (MPN_SAME_OR_SEPARATE_P (rp, vp, n)); 66 1.1 mrg ASSERT (! MPN_OVERLAP_P (rp, n, yp1, n)); 67 1.1 mrg ASSERT (! MPN_OVERLAP_P (rp, n, yp2, n)); 68 1.1 mrg ASSERT (! MPN_OVERLAP_P (rp, n, yp3, n)); 69 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, up, n)); 70 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, vp, n)); 71 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, yp1, n)); 72 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, yp2, n)); 73 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, yp3, n)); 74 1.1 mrg ASSERT (! MPN_OVERLAP_P (ep, 6, rp, n)); 75 1.1 mrg 76 1.1 mrg yp1 += n - 1; 77 1.1 mrg yp2 += n - 1; 78 1.1 mrg yp3 += n - 1; 79 1.1 mrg el1 = eh1 = 0; 80 1.1 mrg el2 = eh2 = 0; 81 1.1 mrg el3 = eh3 = 0; 82 1.1 mrg 83 1.1 mrg do 84 1.1 mrg { 85 1.1 mrg yl1 = *yp1--; 86 1.1 mrg yl2 = *yp2--; 87 1.1 mrg yl3 = *yp3--; 88 1.1 mrg ul = *up++; 89 1.1 mrg vl = *vp++; 90 1.1 mrg 91 1.1 mrg /* ordinary add_n */ 92 1.1 mrg ADDC_LIMB (cy1, sl, ul, vl); 93 1.1 mrg ADDC_LIMB (cy2, rl, sl, cy); 94 1.1 mrg cy = cy1 | cy2; 95 1.1 mrg *rp++ = rl; 96 1.1 mrg 97 1.1 mrg /* update (eh1:el1) */ 98 1.1 mrg zl1 = (-cy) & yl1; 99 1.1 mrg el1 += zl1; 100 1.1 mrg eh1 += el1 < zl1; 101 1.1 mrg 102 1.1 mrg /* update (eh2:el2) */ 103 1.1 mrg zl2 = (-cy) & yl2; 104 1.1 mrg el2 += zl2; 105 1.1 mrg eh2 += el2 < zl2; 106 1.1 mrg 107 1.1 mrg /* update (eh3:el3) */ 108 1.1 mrg zl3 = (-cy) & yl3; 109 1.1 mrg el3 += zl3; 110 1.1 mrg eh3 += el3 < zl3; 111 1.1 mrg } 112 1.1 mrg while (--n); 113 1.1 mrg 114 1.1 mrg #if GMP_NAIL_BITS != 0 115 1.1 mrg eh1 = (eh1 << GMP_NAIL_BITS) + (el1 >> GMP_NUMB_BITS); 116 1.1 mrg el1 &= GMP_NUMB_MASK; 117 1.1 mrg eh2 = (eh2 << GMP_NAIL_BITS) + (el2 >> GMP_NUMB_BITS); 118 1.1 mrg el2 &= GMP_NUMB_MASK; 119 1.1 mrg eh3 = (eh3 << GMP_NAIL_BITS) + (el3 >> GMP_NUMB_BITS); 120 1.1 mrg el3 &= GMP_NUMB_MASK; 121 1.1 mrg #endif 122 1.1 mrg 123 1.1 mrg ep[0] = el1; 124 1.1 mrg ep[1] = eh1; 125 1.1 mrg ep[2] = el2; 126 1.1 mrg ep[3] = eh2; 127 1.1 mrg ep[4] = el3; 128 1.1 mrg ep[5] = eh3; 129 1.1 mrg 130 1.1 mrg return cy; 131 1.1 mrg } 132