Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/aorsmul_i.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz_addmul_ui, mpz_submul_ui - add or subtract small multiple.
2
3
   THE mpz_aorsmul_1 FUNCTION IN THIS FILE IS FOR INTERNAL USE ONLY AND IS
4
   ALMOST CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR
5
   COMPLETELY IN FUTURE GNU MP RELEASES.
6
7
Copyright 2001, 2002, 2004, 2005, 2012, 2021, 2022 Free Software
8
Foundation, Inc.
9
10
This file is part of the GNU MP Library.
11
12
The GNU MP Library is free software; you can redistribute it and/or modify
13
it under the terms of either:
14
15
  * the GNU Lesser General Public License as published by the Free
16
    Software Foundation; either version 3 of the License, or (at your
17
    option) any later version.
18
19
or
20
21
  * the GNU General Public License as published by the Free Software
22
    Foundation; either version 2 of the License, or (at your option) any
23
    later version.
24
25
or both in parallel, as here.
26
27
The GNU MP Library is distributed in the hope that it will be useful, but
28
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
30
for more details.
31
32
You should have received copies of the GNU General Public License and the
33
GNU Lesser General Public License along with the GNU MP Library.  If not,
34
see https://www.gnu.org/licenses/.  */
35
36
#include "gmp-impl.h"
37
38
39
#if HAVE_NATIVE_mpn_mul_1c
40
#define MPN_MUL_1C(cout, dst, src, size, n, cin)        \
41
15.4k
  do {                                                  \
42
15.4k
    (cout) = mpn_mul_1c (dst, src, size, n, cin);       \
43
15.4k
  } while (0)
44
#else
45
#define MPN_MUL_1C(cout, dst, src, size, n, cin)        \
46
  do {                                                  \
47
    mp_limb_t __cy;                                     \
48
    __cy = mpn_mul_1 (dst, src, size, n);               \
49
    (cout) = __cy + mpn_add_1 (dst, dst, size, cin);    \
50
  } while (0)
51
#endif
52
53
54
/* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
55
56
   All that's needed to account for negative w or x is to flip "sub".
57
58
   The final w will retain its sign, unless an underflow occurs in a submul
59
   of absolute values, in which case it's flipped.
60
61
   If x has more limbs than w, then mpn_submul_1 followed by mpn_com is
62
   used.  The alternative would be mpn_mul_1 into temporary space followed
63
   by mpn_sub_n.  Avoiding temporary space seem good, and submul+com stands
64
   a chance of being faster since it involves only one set of carry
65
   propagations, not two.  Note that doing an addmul_1 with a
66
   twos-complement negative y doesn't work, because it effectively adds an
67
   extra x * 2^GMP_LIMB_BITS.  */
68
69
REGPARM_ATTR(1) void
70
mpz_aorsmul_1 (mpz_ptr w, mpz_srcptr x, mp_limb_t y, mp_size_t sub)
71
398k
{
72
398k
  mp_size_t  xsize, wsize, wsize_signed, new_wsize, min_size, dsize;
73
398k
  mp_srcptr  xp;
74
398k
  mp_ptr     wp;
75
398k
  mp_limb_t  cy;
76
77
  /* w unaffected if x==0 or y==0 */
78
398k
  xsize = SIZ (x);
79
398k
  if (xsize == 0 || y == 0)
80
3
    return;
81
82
398k
  sub ^= xsize;
83
398k
  xsize = ABS (xsize);
84
85
398k
  wsize_signed = SIZ (w);
86
398k
  if (wsize_signed == 0)
87
1
    {
88
      /* nothing to add to, just set x*y, "sub" gives the sign */
89
1
      wp = MPZ_NEWALLOC (w, xsize+1);
90
1
      cy = mpn_mul_1 (wp, PTR(x), xsize, y);
91
1
      wp[xsize] = cy;
92
1
      xsize += (cy != 0);
93
1
      SIZ (w) = (sub >= 0 ? xsize : -xsize);
94
1
      return;
95
1
    }
96
97
398k
  sub ^= wsize_signed;
98
398k
  wsize = ABS (wsize_signed);
99
100
398k
  new_wsize = MAX (wsize, xsize);
101
398k
  wp = MPZ_REALLOC (w, new_wsize+1);
102
398k
  xp = PTR (x);
103
398k
  min_size = MIN (wsize, xsize);
104
105
398k
  if (sub >= 0)
106
152k
    {
107
      /* addmul of absolute values */
108
109
152k
      cy = mpn_addmul_1 (wp, xp, min_size, y);
110
152k
      wp += min_size;
111
152k
      xp += min_size;
112
113
152k
      dsize = xsize - wsize;
114
152k
#if HAVE_NATIVE_mpn_mul_1c
115
152k
      if (dsize > 0)
116
2.47k
  cy = mpn_mul_1c (wp, xp, dsize, y, cy);
117
150k
      else if (dsize < 0)
118
2.76k
  {
119
2.76k
    dsize = -dsize;
120
2.76k
    cy = mpn_add_1 (wp, wp, dsize, cy);
121
2.76k
  }
122
#else
123
      if (dsize != 0)
124
  {
125
    mp_limb_t  cy2;
126
    if (dsize > 0)
127
      cy2 = mpn_mul_1 (wp, xp, dsize, y);
128
    else
129
      {
130
        dsize = -dsize;
131
        cy2 = 0;
132
      }
133
    cy = cy2 + mpn_add_1 (wp, wp, dsize, cy);
134
  }
135
#endif
136
137
152k
      wp[dsize] = cy;
138
152k
      new_wsize += (cy != 0);
139
152k
    }
140
246k
  else
141
246k
    {
142
      /* submul of absolute values */
143
144
246k
      cy = mpn_submul_1 (wp, xp, min_size, y);
145
246k
      if (wsize >= xsize)
146
230k
  {
147
    /* if w bigger than x, then propagate borrow through it */
148
230k
    if (wsize != xsize)
149
15.9k
      cy = mpn_sub_1 (wp+xsize, wp+xsize, wsize-xsize, cy);
150
151
230k
    if (cy != 0)
152
149k
      {
153
        /* Borrow out of w, take twos complement negative to get
154
     absolute value, flip sign of w.  */
155
149k
        cy -= mpn_neg (wp, wp, new_wsize);
156
149k
        wp[new_wsize] = cy;
157
149k
        new_wsize += (cy != 0);
158
149k
        wsize_signed = -wsize_signed;
159
149k
      }
160
230k
  }
161
15.4k
      else /* wsize < xsize */
162
15.4k
  {
163
    /* x bigger than w, so want x*y-w.  Submul has given w-x*y, so
164
       take twos complement and use an mpn_mul_1 for the rest.  */
165
166
15.4k
    mp_limb_t  cy2;
167
168
    /* -(-cy*b^n + w-x*y) = (cy-1)*b^n + ~(w-x*y) + 1 */
169
15.4k
    cy -= mpn_neg (wp, wp, wsize);
170
171
    /* If cy-1 == -1 then hold that -1 for latter.  mpn_submul_1 never
172
       returns cy==MP_LIMB_T_MAX so that value always indicates a -1. */
173
15.4k
    cy2 = (cy == MP_LIMB_T_MAX);
174
15.4k
    cy += cy2;
175
15.4k
    MPN_MUL_1C (cy, wp+wsize, xp+wsize, xsize-wsize, y, cy);
176
15.4k
    wp[new_wsize] = cy;
177
15.4k
    new_wsize += (cy != 0);
178
179
    /* Apply any -1 from above.  The value at wp+wsize is non-zero
180
       because y!=0 and the high limb of x will be non-zero.  */
181
15.4k
    if (cy2)
182
2.68k
      MPN_DECR_U (wp+wsize, new_wsize-wsize, CNST_LIMB(1));
183
184
15.4k
    wsize_signed = -wsize_signed;
185
15.4k
  }
186
187
      /* submul can produce high zero limbs due to cancellation, both when w
188
   has more limbs or x has more  */
189
246k
      MPN_NORMALIZE (wp, new_wsize);
190
246k
    }
191
192
398k
  SIZ (w) = (wsize_signed >= 0 ? new_wsize : -new_wsize);
193
194
398k
  ASSERT (new_wsize == 0 || PTR(w)[new_wsize-1] != 0);
195
398k
}
196
197
198
void
199
mpz_addmul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
200
152k
{
201
#if BITS_PER_ULONG > GMP_NUMB_BITS
202
  if (UNLIKELY (y > GMP_NUMB_MAX))
203
    {
204
      mpz_t t;
205
      mp_ptr tp;
206
      mp_size_t xn;
207
      TMP_DECL;
208
      TMP_MARK;
209
      xn = SIZ (x);
210
      if (xn == 0) return;
211
      MPZ_TMP_INIT (t, ABS (xn) + 1);
212
      tp = PTR (t);
213
      tp[0] = 0;
214
      MPN_COPY (tp + 1, PTR(x), ABS (xn));
215
      SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
216
      mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) 0);
217
      PTR(t) = tp + 1;
218
      SIZ(t) = xn;
219
      mpz_aorsmul_1 (w, t, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) 0);
220
      TMP_FREE;
221
      return;
222
    }
223
#endif
224
152k
  mpz_aorsmul_1 (w, x, (mp_limb_t) y, (mp_size_t) 0);
225
152k
}
226
227
void
228
mpz_submul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
229
246k
{
230
#if BITS_PER_ULONG > GMP_NUMB_BITS
231
  if (y > GMP_NUMB_MAX)
232
    {
233
      mpz_t t;
234
      mp_ptr tp;
235
      mp_size_t xn;
236
      TMP_DECL;
237
      TMP_MARK;
238
      xn = SIZ (x);
239
      if (xn == 0) return;
240
      MPZ_TMP_INIT (t, ABS (xn) + 1);
241
      tp = PTR (t);
242
      tp[0] = 0;
243
      MPN_COPY (tp + 1, PTR(x), ABS (xn));
244
      SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
245
      mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) -1);
246
      PTR(t) = tp + 1;
247
      SIZ(t) = xn;
248
      mpz_aorsmul_1 (w, t, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) -1);
249
      TMP_FREE;
250
      return;
251
    }
252
#endif
253
246k
  mpz_aorsmul_1 (w, x, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) -1);
254
246k
}