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