1 /* doaddsub.c: bcmath library file. */
2 /*
3 Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
4 Copyright (C) 2000 Philip A. Nelson
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2 of the License, or (at your option) any later version.
10
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details. (LICENSE)
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to:
18
19 The Free Software Foundation, Inc.
20 59 Temple Place, Suite 330
21 Boston, MA 02111-1307 USA.
22
23 You may contact the author by:
24 e-mail: philnelson@acm.org
25 us-mail: Philip A. Nelson
26 Computer Science Department, 9062
27 Western Washington University
28 Bellingham, WA 98226-9062
29
30 *************************************************************************/
31
32 #include "bcmath.h"
33 #include "private.h"
34 #include <stddef.h>
35
36
37 /* Perform addition: N1 is added to N2 and the value is
38 returned. The signs of N1 and N2 are ignored.
39 SCALE_MIN is to set the minimum scale of the result. */
40
_bc_do_add(bc_num n1,bc_num n2)41 bc_num _bc_do_add(bc_num n1, bc_num n2)
42 {
43 bc_num sum;
44 size_t sum_len = MAX(n1->n_len, n2->n_len) + 1;
45 size_t sum_scale = MAX(n1->n_scale, n2->n_scale);
46 size_t min_len = MIN (n1->n_len, n2->n_len);
47 size_t min_scale = MIN(n1->n_scale, n2->n_scale);
48 size_t min_bytes = min_len + min_scale;
49 char *n1ptr, *n2ptr, *sumptr;
50 bool carry = 0;
51 size_t count;
52
53 /* Prepare sum. */
54 sum = bc_new_num_nonzeroed(sum_len, sum_scale);
55
56 /* Start with the fraction part. Initialize the pointers. */
57 n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale - 1);
58 n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale - 1);
59 sumptr = (char *) (sum->n_value + sum_scale + sum_len - 1);
60
61 /* Add the fraction part. First copy the longer fraction.*/
62 if (n1->n_scale != min_scale) {
63 /* n1 has the longer scale */
64 for (count = n1->n_scale - min_scale; count > 0; count--) {
65 *sumptr-- = *n1ptr--;
66 }
67 } else {
68 /* n2 has the longer scale */
69 for (count = n2->n_scale - min_scale; count > 0; count--) {
70 *sumptr-- = *n2ptr--;
71 }
72 }
73
74 /* Now add the remaining fraction part and equal size integer parts. */
75 count = 0;
76 /* Uses SIMD to perform calculations at high speed. */
77 if (min_bytes >= sizeof(BC_VECTOR)) {
78 sumptr++;
79 n1ptr++;
80 n2ptr++;
81 while (count + sizeof(BC_VECTOR) <= min_bytes) {
82 sumptr -= sizeof(BC_VECTOR);
83 n1ptr -= sizeof(BC_VECTOR);
84 n2ptr -= sizeof(BC_VECTOR);
85
86 BC_VECTOR n1bytes;
87 BC_VECTOR n2bytes;
88 memcpy(&n1bytes, n1ptr, sizeof(n1bytes));
89 memcpy(&n2bytes, n2ptr, sizeof(n2bytes));
90
91 #if BC_LITTLE_ENDIAN
92 /* Little endian requires changing the order of bytes. */
93 n1bytes = BC_BSWAP(n1bytes);
94 n2bytes = BC_BSWAP(n2bytes);
95 #endif
96
97 /*
98 * In order to add 1 to the "next digit" when a carry occurs, adjust it so that it
99 * overflows when add 10.
100 * e.g.
101 * 00001001(9) + 00000001(1) = 00001010(10) to
102 * 11111111 + 00000001 = 00000000(0) and carry 1
103 */
104 n1bytes += SWAR_REPEAT(0xF6) + n2bytes + carry;
105 /* If the most significant bit is 0, a carry has occurred. */
106 carry = !(n1bytes & ((BC_VECTOR) 1 << (8 * sizeof(BC_VECTOR) - 1)));
107
108 /*
109 * The calculation result is a mixture of bytes that have been carried and bytes that have not.
110 * The most significant bit of each byte is 0 if it is carried forward, and 1 if it is not.
111 * Using this, subtract the 0xF6 added for adjustment from the byte that has not been carried
112 * over to return it to the correct value as a decimal number.
113 */
114 BC_VECTOR sum_mask = ((n1bytes & SWAR_REPEAT(0x80)) >> 7) * 0xF6;
115 n1bytes -= sum_mask;
116
117 #if BC_LITTLE_ENDIAN
118 /* Little endian requires changing the order of bytes back. */
119 n1bytes = BC_BSWAP(n1bytes);
120 #endif
121
122 memcpy(sumptr, &n1bytes, sizeof(n1bytes));
123
124 count += sizeof(BC_VECTOR);
125 }
126 sumptr--;
127 n1ptr--;
128 n2ptr--;
129 }
130
131 for (; count < min_bytes; count++) {
132 *sumptr = *n1ptr-- + *n2ptr-- + carry;
133 if (*sumptr >= BASE) {
134 *sumptr -= BASE;
135 carry = 1;
136 } else {
137 carry = 0;
138 }
139 sumptr--;
140 }
141
142 /* Now add carry the longer integer part. */
143 if (n1->n_len != n2->n_len) {
144 if (n2->n_len > n1->n_len) {
145 n1ptr = n2ptr;
146 }
147 for (count = sum_len - min_len; count > 1; count--) {
148 *sumptr = *n1ptr-- + carry;
149 if (*sumptr >= BASE) {
150 *sumptr -= BASE;
151 carry = 1;
152 } else {
153 carry = 0;
154 }
155 sumptr--;
156 }
157 }
158
159 /* Set final carry. */
160 *sumptr = carry;
161
162 /* Adjust sum and return. */
163 _bc_rm_leading_zeros(sum);
164 return sum;
165 }
166
167
168 /* Perform subtraction: N2 is subtracted from N1 and the value is
169 returned. The signs of N1 and N2 are ignored. Also, N1 is
170 assumed to be larger than N2. SCALE_MIN is the minimum scale
171 of the result. */
_bc_do_sub(bc_num n1,bc_num n2)172 bc_num _bc_do_sub(bc_num n1, bc_num n2)
173 {
174 bc_num diff;
175 /* The caller is guaranteed that n1 is always large. */
176 size_t diff_len = EXPECTED(n1->n_len >= n2->n_len) ? n1->n_len : n2->n_len;
177 size_t diff_scale = MAX(n1->n_scale, n2->n_scale);
178 /* Same condition as EXPECTED before, but using EXPECTED again will make it slower. */
179 size_t min_len = n1->n_len >= n2->n_len ? n2->n_len : n1->n_len;
180 size_t min_scale = MIN(n1->n_scale, n2->n_scale);
181 size_t min_bytes = min_len + min_scale;
182 size_t borrow = 0;
183 size_t count;
184 int val;
185 char *n1ptr, *n2ptr, *diffptr;
186
187 /* Allocate temporary storage. */
188 diff = bc_new_num_nonzeroed(diff_len, diff_scale);
189
190 /* Initialize the subtract. */
191 n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale - 1);
192 n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale - 1);
193 diffptr = (char *) (diff->n_value + diff_len + diff_scale - 1);
194
195 /* Take care of the longer scaled number. */
196 if (n1->n_scale != min_scale) {
197 /* n1 has the longer scale */
198 for (count = n1->n_scale - min_scale; count > 0; count--) {
199 *diffptr-- = *n1ptr--;
200 }
201 } else {
202 /* n2 has the longer scale */
203 for (count = n2->n_scale - min_scale; count > 0; count--) {
204 val = -*n2ptr-- - borrow;
205 if (val < 0) {
206 val += BASE;
207 borrow = 1;
208 } else {
209 borrow = 0;
210 }
211 *diffptr-- = val;
212 }
213 }
214
215 /* Now do the equal length scale and integer parts. */
216 count = 0;
217 /* Uses SIMD to perform calculations at high speed. */
218 if (min_bytes >= sizeof(BC_VECTOR)) {
219 diffptr++;
220 n1ptr++;
221 n2ptr++;
222 while (count + sizeof(BC_VECTOR) <= min_bytes) {
223 diffptr -= sizeof(BC_VECTOR);
224 n1ptr -= sizeof(BC_VECTOR);
225 n2ptr -= sizeof(BC_VECTOR);
226
227 BC_VECTOR n1bytes;
228 BC_VECTOR n2bytes;
229 memcpy(&n1bytes, n1ptr, sizeof(n1bytes));
230 memcpy(&n2bytes, n2ptr, sizeof(n2bytes));
231
232 #if BC_LITTLE_ENDIAN
233 /* Little endian requires changing the order of bytes. */
234 n1bytes = BC_BSWAP(n1bytes);
235 n2bytes = BC_BSWAP(n2bytes);
236 #endif
237
238 n1bytes -= n2bytes + borrow;
239 /* If the most significant bit is 1, a carry down has occurred. */
240 bool tmp_borrow = n1bytes & ((BC_VECTOR) 1 << (8 * sizeof(BC_VECTOR) - 1));
241
242 /*
243 * Check the most significant bit of each of the bytes, and if it is 1, a carry down has
244 * occurred. When carrying down occurs, due to the difference between decimal and hexadecimal
245 * numbers, an extra 6 is added to the lower 4 bits.
246 * Therefore, for a byte that has been carried down, set all the upper 4 bits to 0 and subtract
247 * 6 from the lower 4 bits to adjust it to the correct value as a decimal number.
248 */
249 BC_VECTOR borrow_mask = ((n1bytes & SWAR_REPEAT(0x80)) >> 7) * 0x06;
250 n1bytes = (n1bytes & SWAR_REPEAT(0x0F)) - borrow_mask;
251
252 #if BC_LITTLE_ENDIAN
253 /* Little endian requires changing the order of bytes back. */
254 n1bytes = BC_BSWAP(n1bytes);
255 #endif
256
257 memcpy(diffptr, &n1bytes, sizeof(n1bytes));
258
259 borrow = tmp_borrow;
260 count += sizeof(BC_VECTOR);
261 }
262 diffptr--;
263 n1ptr--;
264 n2ptr--;
265 }
266
267 /* Calculate the remaining bytes that are less than the size of BC_VECTOR using a normal loop. */
268 for (; count < min_bytes; count++) {
269 val = *n1ptr-- - *n2ptr-- - borrow;
270 if (val < 0) {
271 val += BASE;
272 borrow = 1;
273 } else {
274 borrow = 0;
275 }
276 *diffptr-- = val;
277 }
278
279 /* If n1 has more digits than n2, we now do that subtract. */
280 if (diff_len != min_len) {
281 for (count = diff_len - min_len; count > 0; count--) {
282 val = *n1ptr-- - borrow;
283 if (val < 0) {
284 val += BASE;
285 borrow = 1;
286 } else {
287 borrow = 0;
288 }
289 *diffptr-- = val;
290 }
291 }
292
293 /* Clean up and return. */
294 _bc_rm_leading_zeros(diff);
295 return diff;
296 }
297
