summaryrefslogtreecommitdiff
path: root/mpn/generic/sub_err3_n.c
blob: 0d665e139fde5b9411b1341a95755a865a65ee22 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
/* mpn_sub_err3_n -- sub_n with three error terms

   Contributed by David Harvey.

   THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   GUARANTEED THAT IT'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.

Copyright 2011 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see https://www.gnu.org/licenses/.  */

#include "gmp.h"
#include "gmp-impl.h"

/*
  Computes:

  (1) {rp,n} := {up,n} - {vp,n} (just like mpn_sub_n) with incoming borrow cy,
  return value is borrow out.

  (2) Let c[i+1] = borrow from i-th limb subtraction (c[0] = cy).
  Computes c[1]*yp1[n-1] + ... + c[n]*yp1[0],
           c[1]*yp2[n-1] + ... + c[n]*yp2[0],
           c[1]*yp3[n-1] + ... + c[n]*yp3[0],
  stores two-limb results at {ep,2}, {ep+2,2} and {ep+4,2} respectively.

  Requires n >= 1.

  None of the outputs may overlap each other or any of the inputs, except
  that {rp,n} may be equal to {up,n} or {vp,n}.
*/
mp_limb_t
mpn_sub_err3_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
                mp_ptr ep, mp_srcptr yp1, mp_srcptr yp2, mp_srcptr yp3,
                mp_size_t n, mp_limb_t cy)
{
  mp_limb_t el1, eh1, el2, eh2, el3, eh3, ul, vl, yl1, yl2, yl3, zl1, zl2, zl3, rl, sl, cy1, cy2;

  ASSERT (n >= 1);
  ASSERT (MPN_SAME_OR_SEPARATE_P (rp, up, n));
  ASSERT (MPN_SAME_OR_SEPARATE_P (rp, vp, n));
  ASSERT (! MPN_OVERLAP_P (rp, n, yp1, n));
  ASSERT (! MPN_OVERLAP_P (rp, n, yp2, n));
  ASSERT (! MPN_OVERLAP_P (rp, n, yp3, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, up, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, vp, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, yp1, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, yp2, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, yp3, n));
  ASSERT (! MPN_OVERLAP_P (ep, 6, rp, n));

  yp1 += n - 1;
  yp2 += n - 1;
  yp3 += n - 1;
  el1 = eh1 = 0;
  el2 = eh2 = 0;
  el3 = eh3 = 0;

  do
    {
      yl1 = *yp1--;
      yl2 = *yp2--;
      yl3 = *yp3--;
      ul = *up++;
      vl = *vp++;

      /* ordinary sub_n */
      SUBC_LIMB (cy1, sl, ul, vl);
      SUBC_LIMB (cy2, rl, sl, cy);
      cy = cy1 | cy2;
      *rp++ = rl;

      /* update (eh1:el1) */
      zl1 = (-cy) & yl1;
      el1 += zl1;
      eh1 += el1 < zl1;

      /* update (eh2:el2) */
      zl2 = (-cy) & yl2;
      el2 += zl2;
      eh2 += el2 < zl2;

      /* update (eh3:el3) */
      zl3 = (-cy) & yl3;
      el3 += zl3;
      eh3 += el3 < zl3;
    }
  while (--n);

#if GMP_NAIL_BITS != 0
  eh1 = (eh1 << GMP_NAIL_BITS) + (el1 >> GMP_NUMB_BITS);
  el1 &= GMP_NUMB_MASK;
  eh2 = (eh2 << GMP_NAIL_BITS) + (el2 >> GMP_NUMB_BITS);
  el2 &= GMP_NUMB_MASK;
  eh3 = (eh3 << GMP_NAIL_BITS) + (el3 >> GMP_NUMB_BITS);
  el3 &= GMP_NUMB_MASK;
#endif

  ep[0] = el1;
  ep[1] = eh1;
  ep[2] = el2;
  ep[3] = eh2;
  ep[4] = el3;
  ep[5] = eh3;

  return cy;
}