xref: /onnv-gate/usr/src/lib/libbc/libc/gen/common/_times_power.c (revision 722:636b850d4ee9)
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