1 /* mpz_addmul_ui, mpz_submul_ui - add or subtract small multiple.
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.
7 Copyright 2001, 2002, 2004, 2005 Free Software Foundation, Inc.
9 This file is part of the GNU MP Library.
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
28 #if HAVE_NATIVE_mpn_mul_1c
29 #define MPN_MUL_1C(cout, dst, src, size, n, cin) \
31 (cout) = mpn_mul_1c (dst, src, size, n, cin); \
34 #define MPN_MUL_1C(cout, dst, src, size, n, cin) \
37 __cy = mpn_mul_1 (dst, src, size, n); \
38 (cout) = __cy + mpn_add_1 (dst, dst, size, cin); \
43 /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
45 All that's needed to account for negative w or x is to flip "sub".
47 The final w will retain its sign, unless an underflow occurs in a submul
48 of absolute values, in which case it's flipped.
50 If x has more limbs than w, then mpn_submul_1 followed by mpn_com is
51 used. The alternative would be mpn_mul_1 into temporary space followed
52 by mpn_sub_n. Avoiding temporary space seem good, and submul+com stands
53 a chance of being faster since it involves only one set of carry
54 propagations, not two. Note that doing an addmul_1 with a
55 twos-complement negative y doesn't work, because it effectively adds an
56 extra x * 2^GMP_LIMB_BITS. */
59 mpz_aorsmul_1 (mpz_ptr w
, mpz_srcptr x
, mp_limb_t y
, mp_size_t sub
)
61 mp_size_t xsize
, wsize
, wsize_signed
, new_wsize
, min_size
, dsize
;
66 /* w unaffected if x==0 or y==0 */
68 if (xsize
== 0 || y
== 0)
74 wsize_signed
= SIZ (w
);
75 if (wsize_signed
== 0)
77 /* nothing to add to, just set x*y, "sub" gives the sign */
78 MPZ_REALLOC (w
, xsize
+1);
80 cy
= mpn_mul_1 (wp
, PTR(x
), xsize
, y
);
83 SIZ (w
) = (sub
>= 0 ? xsize
: -xsize
);
88 wsize
= ABS (wsize_signed
);
90 new_wsize
= MAX (wsize
, xsize
);
91 MPZ_REALLOC (w
, new_wsize
+1);
94 min_size
= MIN (wsize
, xsize
);
98 /* addmul of absolute values */
100 cy
= mpn_addmul_1 (wp
, xp
, min_size
, y
);
104 dsize
= xsize
- wsize
;
105 #if HAVE_NATIVE_mpn_mul_1c
107 cy
= mpn_mul_1c (wp
, xp
, dsize
, y
, cy
);
111 cy
= mpn_add_1 (wp
, wp
, dsize
, cy
);
118 cy2
= mpn_mul_1 (wp
, xp
, dsize
, y
);
124 cy
= cy2
+ mpn_add_1 (wp
, wp
, dsize
, cy
);
129 new_wsize
+= (cy
!= 0);
133 /* submul of absolute values */
135 cy
= mpn_submul_1 (wp
, xp
, min_size
, y
);
138 /* if w bigger than x, then propagate borrow through it */
140 cy
= mpn_sub_1 (wp
+xsize
, wp
+xsize
, wsize
-xsize
, cy
);
144 /* Borrow out of w, take twos complement negative to get
145 absolute value, flip sign of w. */
146 wp
[new_wsize
] = ~-cy
; /* extra limb is 0-cy */
147 mpn_com (wp
, wp
, new_wsize
);
149 MPN_INCR_U (wp
, new_wsize
, CNST_LIMB(1));
150 wsize_signed
= -wsize_signed
;
153 else /* wsize < xsize */
155 /* x bigger than w, so want x*y-w. Submul has given w-x*y, so
156 take twos complement and use an mpn_mul_1 for the rest. */
160 /* -(-cy*b^n + w-x*y) = (cy-1)*b^n + ~(w-x*y) + 1 */
161 mpn_com (wp
, wp
, wsize
);
162 cy
+= mpn_add_1 (wp
, wp
, wsize
, CNST_LIMB(1));
165 /* If cy-1 == -1 then hold that -1 for latter. mpn_submul_1 never
166 returns cy==MP_LIMB_T_MAX so that value always indicates a -1. */
167 cy2
= (cy
== MP_LIMB_T_MAX
);
169 MPN_MUL_1C (cy
, wp
+wsize
, xp
+wsize
, xsize
-wsize
, y
, cy
);
171 new_wsize
+= (cy
!= 0);
173 /* Apply any -1 from above. The value at wp+wsize is non-zero
174 because y!=0 and the high limb of x will be non-zero. */
176 MPN_DECR_U (wp
+wsize
, new_wsize
-wsize
, CNST_LIMB(1));
178 wsize_signed
= -wsize_signed
;
181 /* submul can produce high zero limbs due to cancellation, both when w
182 has more limbs or x has more */
183 MPN_NORMALIZE (wp
, new_wsize
);
186 SIZ (w
) = (wsize_signed
>= 0 ? new_wsize
: -new_wsize
);
188 ASSERT (new_wsize
== 0 || PTR(w
)[new_wsize
-1] != 0);
193 mpz_addmul_ui (mpz_ptr w
, mpz_srcptr x
, unsigned long y
)
195 #if BITS_PER_ULONG > GMP_NUMB_BITS
196 if (UNLIKELY (y
> GMP_NUMB_MAX
&& SIZ(x
) != 0))
204 MPZ_TMP_INIT (t
, ABS (xn
) + 1);
207 MPN_COPY (tp
+ 1, PTR(x
), ABS (xn
));
208 SIZ(t
) = xn
>= 0 ? xn
+ 1 : xn
- 1;
209 mpz_aorsmul_1 (w
, t
, (mp_limb_t
) y
>> GMP_NUMB_BITS
, (mp_size_t
) 0);
212 mpz_aorsmul_1 (w
, t
, (mp_limb_t
) y
& GMP_NUMB_MASK
, (mp_size_t
) 0);
217 mpz_aorsmul_1 (w
, x
, (mp_limb_t
) y
, (mp_size_t
) 0);
221 mpz_submul_ui (mpz_ptr w
, mpz_srcptr x
, unsigned long y
)
223 #if BITS_PER_ULONG > GMP_NUMB_BITS
224 if (y
> GMP_NUMB_MAX
&& SIZ(x
) != 0)
232 MPZ_TMP_INIT (t
, ABS (xn
) + 1);
235 MPN_COPY (tp
+ 1, PTR(x
), ABS (xn
));
236 SIZ(t
) = xn
>= 0 ? xn
+ 1 : xn
- 1;
237 mpz_aorsmul_1 (w
, t
, (mp_limb_t
) y
>> GMP_NUMB_BITS
, (mp_size_t
) -1);
240 mpz_aorsmul_1 (w
, t
, (mp_limb_t
) y
& GMP_NUMB_MASK
, (mp_size_t
) -1);
245 mpz_aorsmul_1 (w
, x
, (mp_limb_t
) y
& GMP_NUMB_MASK
, (mp_size_t
) -1);