?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
8322 nl: misleading-indentation
[unleashed/tickless.git] / usr / src / lib / libbc / libc / gen / common / _times_power.c
blobe0fb6ae9f1d7fb1af8c80de6bb93d4e52c6e9705
1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22 /*
23 * Copyright 1995 Sun Microsystems, Inc. All rights reserved.
24 * Use is subject to license terms.
25 */
26
27 #include "base_conversion.h"
28 #include <malloc.h>
29
30 void
31 _copy_big_float_digits(_BIG_FLOAT_DIGIT *p1, _BIG_FLOAT_DIGIT *p2,
32 short unsigned n)
33 { /* Copies p1[n] = p2[n] */
34 short unsigned i;
35
36 for (i = 0; i < n; i++)
37 *p1++ = *p2++;
38 }
39
40 void
41 _free_big_float(_big_float *p)
42 {
43 /* Central routine to call free for base conversion. */
44
45 char *freearg = (char *) p;
46
47 (void) free(freearg);
48 #ifdef DEBUG
49 printf(" free called with %X \n", freearg);
50 #endif
51 }
52
53 void
54 _base_conversion_abort(int ern, char *bcastring)
55 {
56 char pstring[160];
57
58 errno = ern;
59 (void) sprintf(pstring, " libc base conversion file %s line %d: %s", __FILE__, __LINE__, bcastring);
60 perror(pstring);
61 abort();
62 }
63
64 /*
65 * function to multiply a big_float times a positive power of two or ten.
66 *
67 * Arguments
68 * pbf: Operand x, to be replaced by the product x * mult ** n.
69 * mult: if mult is two, x is base 10**4;
70 * if mult is ten, x is base 2**16
71 * n:
72 * precision: Number of bits of precision ultimately required
73 * (mult=10) or number of digits of precision ultimately
74 * required (mult=2).
75 * Extra bits are allowed internally to permit correct
76 * rounding.
77 * pnewbf: Return result *pnewbf is set to: pbf if uneventful
78 * BIG_FLOAT_TIMES_TOOBIG if n is bigger than the tables
79 * permit ;
80 * BIG_FLOAT_TIMES_NOMEM if pbf->blength was
81 * insufficient to hold the product, and malloc failed to
82 * produce a new block ;
83 * &newbf if pbf->blength was
84 * insufficient to hold the product, and a new _big_float
85 * was allocated by malloc. newbf holds the product.
86 * It's the caller's responsibility to free this space
87 * when no longer needed.
88 *
89 * if precision is < 0, then -pfb->bexponent/{4 or 16} digits are discarded
90 * on the last product.
91 */
92 void
93 _big_float_times_power(_big_float *pbf, int mult, int n, int precision,
94 _big_float **pnewbf)
95 {
96 short unsigned base, sumlz = 0;
97 unsigned short productsize, trailing_zeros_to_delete, needed_precision, *pp, *table[3], max[3], *start[3], *lz[3], tablepower[3];
98 int i, j, itlast;
99 _big_float *pbfold = pbf;
100 int discard;
101
102 if (precision >= 0)
103 discard = -1;
104 else {
105 precision = -precision;
106 discard = 0;
107 }
108 switch (mult) {
109 case 2: /* *pbf is in base 10**4 so multiply by a
110 * power of two */
111 base = 10;
112 max[0] = _max_tiny_powers_two;
113 max[1] = _max_small_powers_two;
114 max[2] = _max_big_powers_two;
115 table[0] = _tiny_powers_two;
116 table[1] = _small_powers_two;
117 table[2] = _big_powers_two;
118 lz[0] = 0;
119 lz[1] = 0;
120 lz[2] = 0;
121 start[0] = _start_tiny_powers_two;
122 start[1] = _start_small_powers_two;
123 start[2] = _start_big_powers_two;
124 needed_precision = 2 + (precision + 1) / 4; /* Precision is in base
125 * ten; counts round and
126 * sticky. */
127 break;
128 case 10: /* *pbf is in base 2**16 so multiply by a
129 * power of ten */
130 base = 2;
131 max[0] = _max_tiny_powers_ten;
132 max[1] = _max_small_powers_ten;
133 max[2] = _max_big_powers_ten;
134 table[0] = _tiny_powers_ten;
135 table[1] = _small_powers_ten;
136 table[2] = _big_powers_ten;
137 start[0] = _start_tiny_powers_ten;
138 start[1] = _start_small_powers_ten;
139 start[2] = _start_big_powers_ten;
140 lz[0] = _leading_zeros_tiny_powers_ten;
141 lz[1] = _leading_zeros_small_powers_ten;
142 lz[2] = _leading_zeros_big_powers_ten;
143 needed_precision = 2 + (precision + 1) / 16; /* Precision is in base
144 * two; counts round and
145 * sticky. */
146 break;
147 }
148 for (i = 0; i < 3; i++) {
149 tablepower[i] = n % max[i];
150 n = n / max[i];
151 }
152 for (itlast = 2; (itlast >= 0) && (tablepower[itlast] == 0); itlast--);
153 /* Determine last i; could be 0, 1, or 2. */
154 if (n > 0) { /* The tables aren't big enough to accomodate
155 * mult**n, but it doesn't matter since the
156 * result would undoubtedly overflow even
157 * binary quadruple precision format. Return
158 * an error code. */
159 (void) printf("\n _times_power failed due to exponent %d %d %d leftover: %d \n", tablepower[0], tablepower[1], tablepower[2], n);
160 *pnewbf = BIG_FLOAT_TIMES_TOOBIG;
161 goto ret;
162 }
163 productsize = pbf->blength;
164 for (i = 0; i < 3; i++)
165 productsize += (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]];
166
167 if (productsize < needed_precision)
168 needed_precision = productsize;
169
170 if (productsize <= pbf->bsize) {
171 *pnewbf = pbf; /* Work with *pnewbf from now on. */
172 } else { /* Need more significance than *pbf can hold. */
173 char *mallocresult;
174 int mallocarg;
175
176 mallocarg = sizeof(_big_float) + sizeof(_BIG_FLOAT_DIGIT) * (productsize - _BIG_FLOAT_SIZE);
177 mallocresult = malloc(mallocarg);
178 #ifdef DEBUG
179 printf(" malloc arg %X result %X \n", mallocarg, (int) mallocresult);
180 #endif
181 if (mallocresult == (char *) 0) { /* Not enough memory
182 * left, bail out. */
183 *pnewbf = BIG_FLOAT_TIMES_NOMEM;
184 goto ret;
185 }
186 *pnewbf = (_big_float *) mallocresult;
187 _copy_big_float_digits((*pnewbf)->bsignificand, pbf->bsignificand, pbf->blength);
188 (*pnewbf)->blength = pbf->blength;
189 (*pnewbf)->bexponent = pbf->bexponent;
190 pbf = *pnewbf;
191 pbf->bsize = productsize;
192 }
193
194 /* pbf now points to the input and the output big_floats. */
195
196 for (i = 0; i <= itlast; i++)
197 if (tablepower[i] != 0) { /* Step through each of the
198 * tables. */
199 unsigned lengthx, lengthp;
200
201 /* Powers of 10**4 have leading zeros in base 2**16. */
202 lengthp = (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]];
203 lengthx = pbf->blength;
204
205 if (discard >= 0)
206 switch (base) {
207 case 2:
208 discard = (-pbf->bexponent) / 16;
209 break;
210 case 10:
211 discard = (-pbf->bexponent) / 4;
212 break;
213 }
214
215 #ifdef DEBUG
216 {
217 long basexp;
218 int id;
219
220 printf(" step %d x operand length %d \n", i, lengthx);
221 _display_big_float(pbf, base);
222 printf(" step %d p operand length %d power %d \n", i, lengthp, tablepower[i]);
223 basexp = (base == 2) ? (lz[i])[tablepower[i]] : 0;
224 for (id = 0; id < lengthp; id++) {
225 printf("+ %d * ", (table[i])[id + (start[i])[tablepower[i]]]);
226 if (base == 2)
227 printf("2**%d", 16 * (basexp + id));
228 if (base == 10)
229 printf("10**%d", 4 * (basexp + id));
230 if ((id % 4) == 3)
231 printf("\n");
232 }
233 printf("\n");
234 }
235 if ((i == itlast) && (discard >= 0))
236 printf(" alternative discard %d digits \n", discard);
237 #endif
238
239 if (base == 2) {
240 sumlz += (lz[i])[tablepower[i]];
241 pbf->bexponent += 16 * (lz[i])[tablepower[i]];
242 }
243 if (lengthp == 1) { /* Special case - multiply by
244 * <= 10**4 or 2**13 */
245 switch (base) {
246 case 10:
247 _multiply_base_ten_by_two(pbf, tablepower[i]);
248 break;
249 case 2:
250 _multiply_base_two(pbf, (_BIG_FLOAT_DIGIT) ((table[i])[tablepower[i]]), (unsigned long) 0);
251 break;
252 }
253 #ifdef DEBUG
254 assert(pbf->blength <= pbf->bsize);
255 #endif
256 } else if (lengthx == 1) { /* Special case of short
257 * multiplicand. */
258 _BIG_FLOAT_DIGIT multiplier = pbf->bsignificand[0];
259
260 _copy_big_float_digits(pbf->bsignificand, (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]), lengthp);
261 pbf->blength = lengthp;
262 switch (base) {
263 case 10:
264 _multiply_base_ten(pbf, multiplier);
265 break;
266 case 2:
267 _multiply_base_two(pbf, multiplier, (unsigned long) 0);
268 break;
269 }
270 #ifdef DEBUG
271 assert(pbf->blength <= pbf->bsize);
272 #endif
273 } else {/* General case. */
274 short unsigned canquit;
275 short unsigned excess;
276
277 /*
278 * The result will be accumulated in *pbf
279 * from most significant to least
280 * significant.
281 */
282
283 /* Generate criterion for early termination. */
284 switch (base) {
285 case 2:
286 canquit = (short unsigned)65536;
287 break;
288 case 10:
289 canquit = 10000;
290 break;
291 }
292 canquit -= 3 + ((lengthx < lengthp) ? lengthx : lengthp);
293
294 pbf->bsignificand[lengthx + lengthp - 1] = 0; /* Only gets filled by
295 * carries. */
296 for (j = lengthx + lengthp - 2; j >= 0; j--) {
297 int istart = j - lengthp + 1, istop = lengthx - 1;
298 short unsigned lengthprod;
299 _BIG_FLOAT_DIGIT product[3];
300
301 pp = (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]);
302 if (j < istop)
303 istop = j;
304 if (0 > istart)
305 istart = 0;
306
307 switch (base) {
308 case 2:
309 _multiply_base_two_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product);
310 if (product[2] != 0)
311 _carry_propagate_two((unsigned long) product[2], &(pbf->bsignificand[j + 2]));
312 if (product[1] != 0)
313 _carry_propagate_two((unsigned long) product[1], &(pbf->bsignificand[j + 1]));
314 break;
315 case 10:
316 _multiply_base_ten_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product);
317 if (product[2] != 0)
318 _carry_propagate_ten((unsigned long) product[2], &(pbf->bsignificand[j + 2]));
319 if (product[1] != 0)
320 _carry_propagate_ten((unsigned long) product[1], &(pbf->bsignificand[j + 1]));
321 break;
322 }
323 pbf->bsignificand[j] = product[0];
324 lengthprod = lengthx + lengthp;
325 if (pbf->bsignificand[lengthprod - 1] == 0)
326 lengthprod--;
327 if (lengthprod > needed_precision)
328 excess = lengthprod - needed_precision;
329 else
330 excess = 0;
331 if ((i == itlast) && ((j + 2) <= excess) && (pbf->bsignificand[j + 1] <= canquit)
332 && ((pbf->bsignificand[j + 1] | pbf->bsignificand[j]) != 0)) {
333 /*
334 * On the last
335 * multiplication, it's not
336 * necessary to develop the
337 * entire product, if further
338 * digits can't possibly
339 * affect significant digits,
340 * unless there's a chance
341 * the product might be
342 * exact!
343 */
344 /*
345 * Note that the product
346 * might be exact if the j
347 * and j+1 terms are zero; if
348 * they are non-zero, then it
349 * won't be after they're
350 * discarded.
351 */
352
353 excess = j + 2; /* Can discard j+1, j,
354 * ... 0 */
355 #ifdef DEBUG
356 printf(" decided to quit early at j %d since s[j+1] is %d <= %d \n", j, pbf->bsignificand[j + 1], canquit);
357 printf(" s[j+2..j] are %d %d %d \n", pbf->bsignificand[j + 2], pbf->bsignificand[j + 1], pbf->bsignificand[j]);
358 printf(" requested precision %d needed_precision %d big digits out of %d \n", precision, needed_precision, lengthprod);
359 #endif
360 if ((discard >= 0) && ((j + 2) > (discard - (int) sumlz))) {
361 #ifdef DEBUG
362 printf(" early quit rejected because j+2 = %d > %d = discard \n", j + 2, discard);
363 #endif
364 goto pastdiscard;
365 }
366 pbf->bsignificand[excess] |= 1; /* Sticky bit on. */
367 #ifdef DEBUG
368 printf(" discard %d digits - last gets %d \n", excess, pbf->bsignificand[excess]);
369 #endif
370 trailing_zeros_to_delete = excess;
371 goto donegeneral;
372 }
373 pastdiscard: ;
374 #ifdef DEBUG
375 /*
376 * else { printf(" early termination
377 * rejected at j %d since s[j+1] =
378 * %d, canquit = %d \n", j,
379 * pbf->bsignificand[j + 1],
380 * canquit); printf(" s[j+2..j] are
381 * %d %d %d \n", pbf->bsignificand[j
382 * + 2], pbf->bsignificand[j + 1],
383 * pbf->bsignificand[j]); printf("
384 * requested precision %d
385 * needed_precision %d big digits out
386 * of %d \n", precision,
387 * needed_precision, lengthprod); }
388 */
389 #endif
390 }
391 trailing_zeros_to_delete = 0;
392 donegeneral:
393 pbf->blength = lengthx + lengthp;
394 if (pbf->bsignificand[pbf->blength - 1] == 0)
395 pbf->blength--;
396 for (; pbf->bsignificand[trailing_zeros_to_delete] == 0; trailing_zeros_to_delete++);
397 /*
398 * Look for additional trailing zeros to
399 * delete.
400 */
401
402 /*
403 * fix for bug 1070565; if too many trailing
404 * zeroes are deleted, we'll violate the
405 * assertion that bexponent is in [-3,+4]
406 */
407 if (base == 10) {
408 int deletelimit=(1-((pbf->bexponent+3)/4));
409
410 if ((int)trailing_zeros_to_delete > deletelimit) {
411 #ifdef DEBUG
412 printf("\n __x_power trailing zeros delete count lowered from %d to "
413 "%d \n", trailing_zeros_to_delete,deletelimit);
414 #endif
415
416 trailing_zeros_to_delete = deletelimit;
417 }
418 }
419
420
421 if (trailing_zeros_to_delete != 0) {
422 #ifdef DEBUG
423 printf(" %d trailing zeros deleted \n", trailing_zeros_to_delete);
424 #endif
425 _copy_big_float_digits(pbf->bsignificand, &(pbf->bsignificand[trailing_zeros_to_delete]), pbf->blength - trailing_zeros_to_delete);
426 pbf->blength -= trailing_zeros_to_delete;
427 switch (base) {
428 case 2:
429 pbf->bexponent += 16 * trailing_zeros_to_delete;
430 break;
431 case 10:
432 pbf->bexponent += 4 * trailing_zeros_to_delete;
433 break;
434 }
435 }
436 }
437 }
438 if ((pbfold != pbf) && (pbf->blength <= pbfold->bsize)) { /* Don't need that huge
439 * buffer after all! */
440 #ifdef DEBUG
441 printf(" free called from times_power because final length %d <= %d original size \n", pbf->blength, pbfold->bsize);
442 #endif
443
444 /* Copy product to original buffer. */
445 pbfold->blength = pbf->blength;
446 pbfold->bexponent = pbf->bexponent;
447 _copy_big_float_digits(pbfold->bsignificand, pbf->bsignificand, pbf->blength);
448 _free_big_float(*pnewbf); /* Free new buffer. */
449 *pnewbf = pbfold; /* New buffer pointer now agrees with
450 * original. */
451 }
452 ret:
453 return;
454 }
unleashed repo fork