![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to addmul_ui.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz|
Return to addmul_ui.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpz |
1.1 maekawa 1: /* mpz_addmul_ui(prodsum, multiplier, small_multiplicand) --
2: Add MULTIPLICATOR times SMALL_MULTIPLICAND to PRODSUM.
3:
4: Copyright (C) 1997, 2000 Free Software Foundation, Inc.
5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
9: it under the terms of the GNU Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 of the License, or (at your
11: option) any later version.
12:
13: The GNU MP Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Lesser General Public License
19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25:
26: static mp_limb_t mpn_neg1 _PROTO ((mp_ptr, mp_size_t));
27:
28: #if 0
29: #undef MPN_NORMALIZE
30: #define MPN_NORMALIZE(DST, NLIMBS) \
31: do { \
32: while (--(NLIMBS) >= 0 && (DST)[NLIMBS] == 0) \
33: ; \
34: (NLIMBS)++; \
35: } while (0)
36: #undef MPN_NORMALIZE_NOT_ZERO
37: #define MPN_NORMALIZE_NOT_ZERO(DST, NLIMBS) \
38: do { \
39: while ((DST)[--(NLIMBS)] == 0) \
40: ; \
41: (NLIMBS)++; \
42: } while (0)
43: #endif
44:
45: void
46: #if __STDC__
47: mpz_addmul_ui (mpz_ptr rz, mpz_srcptr az, unsigned long int bu)
48: #else
49: mpz_addmul_ui (rz, az, bu)
50: mpz_ptr rz;
51: mpz_srcptr az;
52: unsigned long int bu;
53: #endif
54: {
55: mp_size_t rn, an;
56: mp_ptr rp, ap;
57:
58: an = SIZ (az);
59:
60: /* If either multiplier is zero, result is unaffected. */
61: if (bu == 0 || an == 0)
62: return;
63:
64: rn = SIZ (rz);
65:
66: if (rn == 0)
67: {
68: mp_limb_t cy;
69:
70: an = ABS (an);
71: if (ALLOC (rz) <= an)
72: _mpz_realloc (rz, an + 1);
73: rp = PTR (rz);
74: ap = PTR (az);
75: cy = mpn_mul_1 (rp, ap, an, (mp_limb_t) bu);
76: rp[an] = cy;
77: an += cy != 0;
78: SIZ (rz) = SIZ (az) >= 0 ? an : -an;
79: return;
80: }
81:
82: if ((an ^ rn) >= 0)
83: {
84: /* Sign of operands are the same--really add. */
85: an = ABS (an);
86: rn = ABS (rn);
87: if (rn > an)
88: {
89: mp_limb_t cy;
90: if (ALLOC (rz) <= rn)
91: _mpz_realloc (rz, rn + 1);
92: rp = PTR (rz);
93: ap = PTR (az);
94: cy = mpn_addmul_1 (rp, ap, an, (mp_limb_t) bu);
95: cy = mpn_add_1 (rp + an, rp + an, rn - an, cy);
96: rp[rn] = cy;
97: rn += cy != 0;
98: SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn;
99: return;
100: }
101: else
102: {
103: mp_limb_t cy;
104: if (ALLOC (rz) <= an)
105: _mpz_realloc (rz, an + 1);
106: rp = PTR (rz);
107: ap = PTR (az);
108: cy = mpn_addmul_1 (rp, ap, rn, (mp_limb_t) bu);
109: if (an != rn)
110: {
111: mp_limb_t cy2;
112: cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu);
113: cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy);
114: }
115: rn = an;
116: rp[rn] = cy;
117: rn += cy != 0;
118: SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn;
119: return;
120: }
121: }
122: else
123: {
124: /* Sign of operands are different--actually subtract. */
125: an = ABS (an);
126: rn = ABS (rn);
127: if (rn > an)
128: {
129: mp_limb_t cy;
130: rp = PTR (rz);
131: ap = PTR (az);
132: cy = mpn_submul_1 (rp, ap, an, (mp_limb_t) bu);
133: cy = mpn_sub_1 (rp + an, rp + an, rn - an, cy);
134: if (cy != 0)
135: {
136: mpn_neg1 (rp, rn);
137: MPN_NORMALIZE_NOT_ZERO (rp, rn);
138: }
139: else
140: {
141: MPN_NORMALIZE (rp, rn);
142: rn = -rn;
143: }
144:
145: SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn;
146: return;
147: }
148: else
149: {
150: /* Tricky case. We need to subtract an operand that might be larger
151: than the minuend. To avoid allocating temporary space, we compute
152: a*b-r instead of r-a*b and then negate. */
153: mp_limb_t cy;
154: if (ALLOC (rz) <= an)
155: _mpz_realloc (rz, an + 1);
156: rp = PTR (rz);
157: ap = PTR (az);
158: cy = mpn_submul_1 (rp, ap, rn, (mp_limb_t) bu);
159: if (an != rn)
160: {
161: mp_limb_t cy2;
162: cy -= mpn_neg1 (rp, rn);
163: cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu);
164: if (cy == ~(mp_limb_t) 0)
165: cy = cy2 - mpn_sub_1 (rp + rn, rp + rn, an - rn, (mp_limb_t) 1);
166: else
167: cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy);
168: rp[an] = cy;
169: rn = an + (cy != 0);
170: rn -= rp[rn - 1] == 0;
171: }
172: else if (cy != 0)
173: {
174: cy -= mpn_neg1 (rp, rn);
175: rp[an] = cy;
176: rn = an + 1;
177: MPN_NORMALIZE_NOT_ZERO (rp, rn);
178: }
179: else
180: {
181: rn = an;
182: MPN_NORMALIZE (rp, rn);
183: rn = -rn;
184: }
185:
186: SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn;
187: return;
188: }
189: }
190: }
191:
192: static mp_limb_t
193: #if __STDC__
194: mpn_neg1 (mp_ptr rp, mp_size_t rn)
195: #else
196: mpn_neg1 (rp, rn)
197: mp_ptr rp;
198: mp_size_t rn;
199: #endif
200: {
201: mp_size_t i;
202:
203: while (rn != 0 && rp[0] == 0)
204: rp++, rn--;
205:
206: if (rn != 0)
207: {
208: rp[0] = -rp[0];
209: for (i = 1; i < rn; i++)
210: rp[i] = ~rp[i];
211: return 1;
212: }
213: return 0;
214: }