File: common\bigflt.c
1 /* bigflt.c - C routines for big floating point numbers */
2
3 /*
4 Wesley Loewer's Big Numbers. (C) 1994-95, Wesley B. Loewer
5 */
6
7 #include <string.h>
8 /* see Fractint.c for a description of the "include" hierarchy */
9 #include "port.h"
10 #include "big.h"
11 #ifndef BIG_ANSI_C
12 #include <malloc.h>
13 #endif
14
15 #define LOG10_256 2.4082399653118
16 #define LOG_256 5.5451774444795
17
18 #if (_MSC_VER >= 700) /* for Fractint */
19 #pragma code_seg ("bigflt1_text") /* place following in an overlay */
20 #endif
21
22 /********************************************************************/
23 /* bf_hexdump() - for debugging, dumps to stdout */
24
25 void bf_hexdump(bf_t r)
26 {
27 int i;
28
29 for (i=0; i<bflength; i++)
30 printf("%02X ", *(r+i));
31 printf(" e %04hX ", (S16)big_access16(r+bflength));
32 printf("\n");
33 return;
34 }
35
36 /**********************************************************************/
37 /* strtobf() - converts a string into a bigfloat */
38 /* r - pointer to a bigfloat */
39 /* s - string in the floating point format [+-][dig].[dig]e[+-][dig]*/
40 /* note: the string may not be empty or have extra space. */
41 /* It may use scientific notation. */
42 /* USES: bftmp1 */
43
44 bf_t strtobf(bf_t r, char *s)
45 {
46 BYTE onesbyte;
47 int signflag=0;
48 char *l, *d, *e; /* pointer to s, ".", "[eE]" */
49 int powerten=0, keeplooping;
50
51 clear_bf(r);
52
53 if (s[0] == '+') /* for + sign */
54 {
55 s++;
56 }
57 else if (s[0] == '-') /* for neg sign */
58 {
59 signflag = 1;
60 s++;
61 }
62
63 d = strchr(s, '.');
64 e = strchr(s, 'e');
65 if (e == NULL)
66 e = strchr(s, 'E');
67 if (e != NULL)
68 {
69 powerten = atoi(e+1); /* read in the e (x10^) part */
70 l = e - 1; /* just before e */
71 }
72 else
73 l = s + strlen(s) - 1; /* last digit */
74
75 if (d != NULL) /* is there a decimal point? */
76 {
77 while (*l >= '0' && *l <= '9') /* while a digit */
78 {
79 onesbyte = (BYTE)(*(l--) - '0');
80 inttobf(bftmp1,onesbyte);
81 unsafe_add_a_bf(r, bftmp1);
82 div_a_bf_int(r, 10);
83 }
84
85 if (*(l--) == '.') /* the digit was found */
86 {
87 keeplooping = *l >= '0' && *l <= '9' && l>=s;
88 while (keeplooping) /* while a digit */
89 {
90 onesbyte = (BYTE)(*(l--) - '0');
91 inttobf(bftmp1,onesbyte);
92 unsafe_add_a_bf(r, bftmp1);
93 keeplooping = *l >= '0' && *l <= '9' && l>=s;
94 if (keeplooping)
95 {
96 div_a_bf_int(r, 10);
97 powerten++; /* increase the power of ten */
98 }
99 }
100 }
101 }
102 else
103 {
104 keeplooping = *l >= '0' && *l <= '9' && l>=s;
105 while (keeplooping) /* while a digit */
106 {
107 onesbyte = (BYTE)(*(l--) - '0');
108 inttobf(bftmp1,onesbyte);
109 unsafe_add_a_bf(r, bftmp1);
110 keeplooping = *l >= '0' && *l <= '9' && l>=s;
111 if (keeplooping)
112 {
113 div_a_bf_int(r, 10);
114 powerten++; /* increase the power of ten */
115 }
116 }
117 }
118
119 if (powerten > 0)
120 {
121 for (; powerten>0; powerten--)
122 {
123 mult_a_bf_int(r, 10);
124 }
125 }
126 else if (powerten < 0)
127 {
128 for (; powerten<0; powerten++)
129 {
130 div_a_bf_int(r, 10);
131 }
132 }
133 if (signflag)
134 neg_a_bf(r);
135
136 return r;
137 }
138
139 /********************************************************************/
140 /* strlen_needed() - returns string length needed to hold bigfloat */
141
142 int strlen_needed_bf()
143 {
144 int length;
145
146 /* first space for integer part */
147 length = 1;
148 length += decimals; /* decimal part */
149 length += 2; /* decimal point and sign */
150 length += 2; /* e and sign */
151 length += 4; /* exponent */
152 length += 4; /* null and a little extra for safety */
153 return(length);
154 }
155
156 /********************************************************************/
157 /* bftostr() - converts a bigfloat into a scientific notation string */
158 /* s - string, must be large enough to hold the number. */
159 /* dec - decimal places, 0 for max */
160 /* r - bigfloat */
161 /* will convert to a floating point notation */
162 /* SIDE-EFFECT: the bigfloat, r, is destroyed. */
163 /* Copy it first if necessary. */
164 /* USES: bftmp1 - bftmp2 */
165 /********************************************************************/
166
167 char *unsafe_bftostr(char *s, int dec, bf_t r)
168 {
169 LDBL value;
170 int power;
171
172 value = bftofloat(r);
173 if (value == 0.0)
174 {
175 strcpy(s, "0.0");
176 return s;
177 }
178
179 copy_bf(bftmp1, r);
180 unsafe_bftobf10(bf10tmp, dec, bftmp1);
181 power = (S16)big_access16(bf10tmp+dec+2); /* where the exponent is stored */
182 if (power > -4 && power < 6) /* tinker with this */
183 bf10tostr_f(s, dec, bf10tmp);
184 else
185 bf10tostr_e(s, dec, bf10tmp);
186 return s;
187 }
188
189
190 /********************************************************************/
191 /* the e version puts it in scientific notation, (like printf's %e) */
192 char *unsafe_bftostr_e(char *s, int dec, bf_t r)
193 {
194 LDBL value;
195
196 value = bftofloat(r);
197 if (value == 0.0)
198 {
199 strcpy(s, "0.0");
200 return s;
201 }
202
203 copy_bf(bftmp1, r);
204 unsafe_bftobf10(bf10tmp, dec, bftmp1);
205 bf10tostr_e(s, dec, bf10tmp);
206 return s;
207 }
208
209 /********************************************************************/
210 /* the f version puts it in decimal notation, (like printf's %f) */
211 char *unsafe_bftostr_f(char *s, int dec, bf_t r)
212 {
213 LDBL value;
214
215 value = bftofloat(r);
216 if (value == 0.0)
217 {
218 strcpy(s, "0.0");
219 return s;
220 }
221
222 copy_bf(bftmp1, r);
223 unsafe_bftobf10(bf10tmp, dec, bftmp1);
224 bf10tostr_f(s, dec, bf10tmp);
225 return s;
226 }
227
228 #if (_MSC_VER >= 700)
229 #pragma code_seg ( ) /* back to normal segment */
230 #endif
231
232 /*********************************************************************/
233 /* bn = floor(bf) */
234 /* Converts a bigfloat to a bignumber (integer) */
235 /* bflength must be at least bnlength+2 */
236 bn_t bftobn(bn_t n, bf_t f)
237 {
238 int fexp;
239 int movebytes;
240 BYTE hibyte;
241
242 fexp = (S16)big_access16(f+bflength);
243 if (fexp >= intlength)
244 { /* if it's too big, use max value */
245 max_bn(n);
246 if (is_bf_neg(f))
247 neg_a_bn(n);
248 return n;
249 }
250
251 if (-fexp > bnlength - intlength) /* too small, return zero */
252 {
253 clear_bn(n);
254 return n;
255 }
256
257 /* already checked for over/underflow, this should be ok */
258 movebytes = bnlength - intlength + fexp + 1;
259 _fmemcpy(n, f+bflength-movebytes-1, movebytes);
260 hibyte = *(f+bflength-1);
261 _fmemset(n+movebytes, hibyte, bnlength-movebytes); /* sign extends */
262 return n;
263 }
264
265 /*********************************************************************/
266 /* bf = bn */
267 /* Converts a bignumber (integer) to a bigfloat */
268 /* bflength must be at least bnlength+2 */
269 bf_t bntobf(bf_t f, bn_t n)
270 {
271 _fmemcpy(f+bflength-bnlength-1, n, bnlength);
272 _fmemset(f, 0, bflength - bnlength - 1);
273 *(f+bflength-1) = (BYTE)(is_bn_neg(n) ? 0xFF : 0x00); /* sign extend */
274 big_set16(f+bflength, (S16)(intlength - 1)); /* exp */
275 norm_bf(f);
276 return f;
277 }
278
279 /*********************************************************************/
280 /* b = l */
281 /* Converts a long to a bigfloat */
282 bf_t inttobf(bf_t r, long longval)
283 {
284 clear_bf(r);
285 big_set32(r+bflength-4, (S32)longval);
286 big_set16(r+bflength, (S16)2);
287 norm_bf(r);
288 return r;
289 }
290
291 /*********************************************************************/
292 /* l = floor(b), floor rounds down */
293 /* Converts a bigfloat to a long */
294 /* note: a bf value of 2.999... will be return a value of 2, not 3 */
295 long bftoint(bf_t f)
296 {
297 int fexp;
298 long longval;
299
300 fexp = (S16)big_access16(f+bflength);
301 if (fexp > 3)
302 {
303 longval = 0x7FFFFFFFL;
304 if (is_bf_neg(f))
305 longval = -longval;
306 return longval;
307 }
308 longval = big_access32(f+bflength-5);
309 longval >>= 8*(3-fexp);
310 return longval;
311 }
312
313 /********************************************************************/
314 /* sign(r) */
315 int sign_bf(bf_t n)
316 {
317 return is_bf_neg(n) ? -1 : is_bf_not_zero(n) ? 1 : 0;
318 }
319
320 /********************************************************************/
321 /* r = |n| */
322 bf_t abs_bf(bf_t r, bf_t n)
323 {
324 copy_bf(r,n);
325 if (is_bf_neg(r))
326 {
327 neg_a_bf(r);
328 }
329 return r;
330 }
331
332 /********************************************************************/
333 /* r = |r| */
334 bf_t abs_a_bf(bf_t r)
335 {
336 if (is_bf_neg(r))
337 neg_a_bf(r);
338 return r;
339 }
340
341 /********************************************************************/
342 /* r = 1/n */
343 /* uses bftmp1 - bftmp2 - global temp bigfloats */
344 /* SIDE-EFFECTS: */
345 /* n ends up as |n|/256^exp Make copy first if necessary. */
346 bf_t unsafe_inv_bf(bf_t r, bf_t n)
347 {
348 int signflag=0, i;
349 int fexp, rexp;
350 LDBL f;
351 bf_t orig_r, orig_n; /* orig_bftmp1 not needed here */
352 int orig_bflength,
353 orig_bnlength,
354 orig_padding,
355 orig_rlength,
356 orig_shiftfactor,
357 orig_rbflength;
358
359 /* use Newton's recursive method for zeroing in on 1/n : r=r(2-rn) */
360
361 if (is_bf_neg(n))
362 { /* will be a lot easier to deal with just positives */
363 signflag = 1;
364 neg_a_bf(n);
365 }
366
367 fexp = (S16)big_access16(n+bflength);
368 big_set16(n+bflength, (S16)0); /* put within LDBL range */
369
370 f = bftofloat(n);
371 if (f == 0) /* division by zero */
372 {
373 max_bf(r);
374 return r;
375 }
376 f = 1/f; /* approximate inverse */
377
378 /* With Newton's Method, there is no need to calculate all the digits */
379 /* every time. The precision approximately doubles each iteration. */
380 /* Save original values. */
381 orig_bflength = bflength;
382 orig_bnlength = bnlength;
383 orig_padding = padding;
384 orig_rlength = rlength;
385 orig_shiftfactor = shiftfactor;
386 orig_rbflength = rbflength;
387 orig_r = r;
388 orig_n = n;
389 /* orig_bftmp1 = bftmp1; */
390
391 /* calculate new starting values */
392 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
393 if (bnlength > orig_bnlength)
394 bnlength = orig_bnlength;
395 calc_lengths();
396
397 /* adjust pointers */
398 r = orig_r + orig_bflength - bflength;
399 n = orig_n + orig_bflength - bflength;
400 /* bftmp1 = orig_bftmp1 + orig_bflength - bflength; */
401
402 floattobf(r, f); /* start with approximate inverse */
403
404 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
405 {
406 /* adjust lengths */
407 bnlength <<= 1; /* double precision */
408 if (bnlength > orig_bnlength)
409 bnlength = orig_bnlength;
410 calc_lengths();
411 r = orig_r + orig_bflength - bflength;
412 n = orig_n + orig_bflength - bflength;
413 /* bftmp1 = orig_bftmp1 + orig_bflength - bflength; */
414
415 unsafe_mult_bf(bftmp1, r, n); /* bftmp1=rn */
416 inttobf(bftmp2, 1); /* will be used as 1.0 */
417
418 /* There seems to very little difficulty getting bftmp1 to be EXACTLY 1 */
419 if (bflength == orig_bflength && cmp_bf(bftmp1, bftmp2) == 0)
420 break;
421
422 inttobf(bftmp2, 2); /* will be used as 2.0 */
423 unsafe_sub_a_bf(bftmp2, bftmp1); /* bftmp2=2-rn */
424 unsafe_mult_bf(bftmp1, r, bftmp2); /* bftmp1=r(2-rn) */
425 copy_bf(r, bftmp1); /* r = bftmp1 */
426 }
427
428 /* restore original values */
429 bflength = orig_bflength;
430 bnlength = orig_bnlength;
431 padding = orig_padding;
432 rlength = orig_rlength;
433 shiftfactor = orig_shiftfactor;
434 rbflength = orig_rbflength;
435 r = orig_r;
436 n = orig_n;
437 /* bftmp1 = orig_bftmp1; */
438
439 if (signflag)
440 {
441 neg_a_bf(r);
442 }
443 rexp = (S16)big_access16(r+bflength);
444 rexp -= fexp;
445 big_set16(r+bflength, (S16)rexp); /* adjust result exponent */
446 return r;
447 }
448
449 /********************************************************************/
450 /* r = n1/n2 */
451 /* r - result of length bflength */
452 /* uses bftmp1 - bftmp2 - global temp bigfloats */
453 /* SIDE-EFFECTS: */
454 /* n1, n2 end up as |n1|/256^x, |n2|/256^x */
455 /* Make copies first if necessary. */
456 bf_t unsafe_div_bf(bf_t r, bf_t n1, bf_t n2)
457 {
458 int aexp, bexp, rexp;
459 LDBL a, b;
460
461 /* first, check for valid data */
462
463 aexp = (S16)big_access16(n1+bflength);
464 big_set16(n1+bflength, (S16)0); /* put within LDBL range */
465
466 a = bftofloat(n1);
467 if (a == 0) /* division into zero */
468 {
469 clear_bf(r); /* return 0 */
470 return r;
471 }
472
473 bexp = (S16)big_access16(n2+bflength);
474 big_set16(n2+bflength, (S16)0); /* put within LDBL range */
475
476 b = bftofloat(n2);
477 if (b == 0) /* division by zero */
478 {
479 max_bf(r);
480 return r;
481 }
482
483 unsafe_inv_bf(r, n2);
484 unsafe_mult_bf(bftmp1, n1, r);
485 copy_bf(r, bftmp1); /* r = bftmp1 */
486
487 rexp = (S16)big_access16(r+bflength);
488 rexp += aexp - bexp;
489 big_set16(r+bflength, (S16)rexp); /* adjust result exponent */
490
491 return r;
492 }
493
494 /********************************************************************/
495 /* sqrt(r) */
496 /* uses bftmp1 - bftmp3 - global temp bigfloats */
497 /* SIDE-EFFECTS: */
498 /* n ends up as |n| */
499 bf_t unsafe_sqrt_bf(bf_t r, bf_t n)
500 {
501 int i, comp, almost_match=0;
502 LDBL f;
503 bf_t orig_r, orig_n;
504 int orig_bflength,
505 orig_bnlength,
506 orig_padding,
507 orig_rlength,
508 orig_shiftfactor,
509 orig_rbflength;
510
511 /* use Newton's recursive method for zeroing in on sqrt(n): r=.5(r+n/r) */
512
513 if (is_bf_neg(n))
514 { /* sqrt of a neg, return 0 */
515 clear_bf(r);
516 return r;
517 }
518
519 f = bftofloat(n);
520 if (f == 0) /* division by zero will occur */
521 {
522 clear_bf(r); /* sqrt(0) = 0 */
523 return r;
524 }
525 f = sqrtl(f); /* approximate square root */
526 /* no need to check overflow */
527
528 /* With Newton's Method, there is no need to calculate all the digits */
529 /* every time. The precision approximately doubles each iteration. */
530 /* Save original values. */
531 orig_bflength = bflength;
532 orig_bnlength = bnlength;
533 orig_padding = padding;
534 orig_rlength = rlength;
535 orig_shiftfactor = shiftfactor;
536 orig_rbflength = rbflength;
537 orig_r = r;
538 orig_n = n;
539
540 /* calculate new starting values */
541 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
542 if (bnlength > orig_bnlength)
543 bnlength = orig_bnlength;
544 calc_lengths();
545
546 /* adjust pointers */
547 r = orig_r + orig_bflength - bflength;
548 n = orig_n + orig_bflength - bflength;
549
550 floattobf(r, f); /* start with approximate sqrt */
551
552 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
553 {
554 /* adjust lengths */
555 bnlength <<= 1; /* double precision */
556 if (bnlength > orig_bnlength)
557 bnlength = orig_bnlength;
558 calc_lengths();
559 r = orig_r + orig_bflength - bflength;
560 n = orig_n + orig_bflength - bflength;
561
562 unsafe_div_bf(bftmp3, n, r);
563 unsafe_add_a_bf(r, bftmp3);
564 half_a_bf(r);
565 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp3))) < 8 ) /* if match or almost match */
566 {
567 if (comp < 4 /* perfect or near perfect match */
568 || almost_match == 1) /* close enough for 2nd time */
569 break;
570 else /* this is the first time they almost matched */
571 almost_match++;
572 }
573 }
574
575 /* restore original values */
576 bflength = orig_bflength;
577 bnlength = orig_bnlength;
578 padding = orig_padding;
579 rlength = orig_rlength;
580 shiftfactor = orig_shiftfactor;
581 rbflength = orig_rbflength;
582 r = orig_r;
583 n = orig_n;
584
585 return r;
586 }
587
588 /********************************************************************/
589 /* exp(r) */
590 /* uses bftmp1, bftmp2, bftmp3 - global temp bigfloats */
591 bf_t exp_bf(bf_t r, bf_t n)
592 {
593 U16 fact=1;
594 S16 BIGDIST * testexp, BIGDIST * rexp;
595
596 testexp = (S16 BIGDIST *)(bftmp2+bflength);
597 rexp = (S16 BIGDIST *)(r+bflength);
598
599 if (is_bf_zero(n))
600 {
601 inttobf(r, 1);
602 return r;
603 }
604
605 /* use Taylor Series (very slow convergence) */
606 inttobf(r, 1); /* start with r=1.0 */
607 copy_bf(bftmp2, r);
608 for (;;)
609 {
610 copy_bf(bftmp1, n);
611 unsafe_mult_bf(bftmp3, bftmp2, bftmp1);
612 unsafe_div_bf_int(bftmp2, bftmp3, fact);
613 if (big_accessS16(testexp) < big_accessS16(rexp)-(bflength-2))
614 break; /* too small to register */
615 unsafe_add_a_bf(r, bftmp2);
616 fact++;
617 }
618
619 return r;
620 }
621
622 /********************************************************************/
623 /* ln(r) */
624 /* uses bftmp1 - bftmp6 - global temp bigfloats */
625 /* SIDE-EFFECTS: */
626 /* n ends up as |n| */
627 bf_t unsafe_ln_bf(bf_t r, bf_t n)
628 {
629 int i, comp, almost_match=0;
630 LDBL f;
631 bf_t orig_r, orig_n, orig_bftmp5;
632 int orig_bflength,
633 orig_bnlength,
634 orig_padding,
635 orig_rlength,
636 orig_shiftfactor,
637 orig_rbflength;
638
639 /* use Newton's recursive method for zeroing in on ln(n): r=r+n*exp(-r)-1 */
640
641 if (is_bf_neg(n) || is_bf_zero(n))
642 { /* error, return largest neg value */
643 max_bf(r);
644 neg_a_bf(r);
645 return r;
646 }
647
648 f = bftofloat(n);
649 f = logl(f); /* approximate ln(x) */
650 /* no need to check overflow */
651 /* appears to be ok, do ln */
652
653 /* With Newton's Method, there is no need to calculate all the digits */
654 /* every time. The precision approximately doubles each iteration. */
655 /* Save original values. */
656 orig_bflength = bflength;
657 orig_bnlength = bnlength;
658 orig_padding = padding;
659 orig_rlength = rlength;
660 orig_shiftfactor = shiftfactor;
661 orig_rbflength = rbflength;
662 orig_r = r;
663 orig_n = n;
664 orig_bftmp5 = bftmp5;
665
666 /* calculate new starting values */
667 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
668 if (bnlength > orig_bnlength)
669 bnlength = orig_bnlength;
670 calc_lengths();
671
672 /* adjust pointers */
673 r = orig_r + orig_bflength - bflength;
674 n = orig_n + orig_bflength - bflength;
675 bftmp5 = orig_bftmp5 + orig_bflength - bflength;
676
677 floattobf(r, f); /* start with approximate ln */
678 neg_a_bf(r); /* -r */
679 copy_bf(bftmp5, r); /* -r */
680
681 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
682 {
683 /* adjust lengths */
684 bnlength <<= 1; /* double precision */
685 if (bnlength > orig_bnlength)
686 bnlength = orig_bnlength;
687 calc_lengths();
688 r = orig_r + orig_bflength - bflength;
689 n = orig_n + orig_bflength - bflength;
690 bftmp5 = orig_bftmp5 + orig_bflength - bflength;
691
692 exp_bf(bftmp6, r); /* exp(-r) */
693 unsafe_mult_bf(bftmp2, bftmp6, n); /* n*exp(-r) */
694 inttobf(bftmp4, 1);
695 unsafe_sub_a_bf(bftmp2, bftmp4); /* n*exp(-r) - 1 */
696 unsafe_sub_a_bf(r, bftmp2); /* -r - (n*exp(-r) - 1) */
697 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp5))) < 8 ) /* if match or almost match */
698 {
699 if (comp < 4 /* perfect or near perfect match */
700 || almost_match == 1) /* close enough for 2nd time */
701 break;
702 else /* this is the first time they almost matched */
703 almost_match++;
704 }
705 copy_bf(bftmp5, r); /* -r */
706 }
707
708 /* restore original values */
709 bflength = orig_bflength;
710 bnlength = orig_bnlength;
711 padding = orig_padding;
712 rlength = orig_rlength;
713 shiftfactor = orig_shiftfactor;
714 rbflength = orig_rbflength;
715 r = orig_r;
716 n = orig_n;
717 bftmp5 = orig_bftmp5;
718
719 neg_a_bf(r); /* -(-r) */
720 return r;
721 }
722
723 /********************************************************************/
724 /* sincos_bf(r) */
725 /* uses bftmp1 - bftmp2 - global temp bigfloats */
726 /* SIDE-EFFECTS: */
727 /* n ends up as |n| mod (pi/4) */
728 bf_t unsafe_sincos_bf(bf_t s, bf_t c, bf_t n)
729 {
730 U16 fact=2;
731 int k=0, i, halves;
732 int signcos=0, signsin=0, switch_sincos=0;
733 int sin_done=0, cos_done=0;
734 S16 BIGDIST * testexp, BIGDIST * cexp, BIGDIST * sexp;
735
736 testexp = (S16 BIGDIST *)(bftmp1+bflength);
737 cexp = (S16 BIGDIST *)(c+bflength);
738 sexp = (S16 BIGDIST *)(s+bflength);
739
740 #ifndef CALCULATING_BIG_PI
741 /* assure range 0 <= x < pi/4 */
742
743 if (is_bf_zero(n))
744 {
745 clear_bf(s); /* sin(0) = 0 */
746 inttobf(c, 1); /* cos(0) = 1 */
747 return s;
748 }
749
750 if (is_bf_neg(n))
751 {
752 signsin = !signsin; /* sin(-x) = -sin(x), odd; cos(-x) = cos(x), even */
753 neg_a_bf(n);
754 }
755 /* n >= 0 */
756
757 double_bf(bftmp1, bf_pi); /* 2*pi */
758 /* this could be done with remainders, but it would probably be slower */
759 while (cmp_bf(n, bftmp1) >= 0) /* while n >= 2*pi */
760 {
761 copy_bf(bftmp2, bftmp1);
762 unsafe_sub_a_bf(n, bftmp2);
763 }
764 /* 0 <= n < 2*pi */
765
766 copy_bf(bftmp1, bf_pi); /* pi */
767 if (cmp_bf(n, bftmp1) >= 0) /* if n >= pi */
768 {
769 unsafe_sub_a_bf(n, bftmp1);
770 signsin = !signsin;
771 signcos = !signcos;
772 }
773 /* 0 <= n < pi */
774
775 half_bf(bftmp1, bf_pi); /* pi/2 */
776 if (cmp_bf(n, bftmp1) > 0) /* if n > pi/2 */
777 {
778 copy_bf(bftmp2, bf_pi);
779 unsafe_sub_bf(n, bftmp2, n);
780 signcos = !signcos;
781 }
782 /* 0 <= n < pi/2 */
783
784 half_bf(bftmp1, bf_pi); /* pi/2 */
785 half_a_bf(bftmp1); /* pi/4 */
786 if (cmp_bf(n, bftmp1) > 0) /* if n > pi/4 */
787 {
788 copy_bf(bftmp2, n);
789 half_bf(bftmp1, bf_pi); /* pi/2 */
790 unsafe_sub_bf(n, bftmp1, bftmp2); /* pi/2 - n */
791 switch_sincos = !switch_sincos;
792 }
793 /* 0 <= n < pi/4 */
794
795 /* this looks redundant, but n could now be zero when it wasn't before */
796 if (is_bf_zero(n))
797 {
798 clear_bf(s); /* sin(0) = 0 */
799 inttobf(c, 1); /* cos(0) = 1 */
800 return s;
801 }
802
803
804 /* at this point, the double angle trig identities could be used as many */
805 /* times as desired to reduce the range to pi/8, pi/16, etc... Each time */
806 /* the range is cut in half, the number of iterations required is reduced */
807 /* by "quite a bit." It's just a matter of testing to see what gives the */
808 /* optimal results. */
809 /* halves = bflength / 10; */ /* this is experimental */
810 halves = 1;
811 for (i = 0; i < halves; i++)
812 half_a_bf(n);
813 #endif
814
815 /* use Taylor Series (very slow convergence) */
816 copy_bf(s, n); /* start with s=n */
817 inttobf(c, 1); /* start with c=1 */
818 copy_bf(bftmp1, n); /* the current x^n/n! */
819 do
820 {
821 /* even terms for cosine */
822 copy_bf(bftmp2, bftmp1);
823 unsafe_mult_bf(bftmp1, bftmp2, n);
824 div_a_bf_int(bftmp1, fact++);
825 if (!cos_done)
826 {
827 cos_done = (big_accessS16(testexp) < big_accessS16(cexp)-(bflength-2)); /* too small to register */
828 if (!cos_done)
829 {
830 if (k) /* alternate between adding and subtracting */
831 unsafe_add_a_bf(c, bftmp1);
832 else
833 unsafe_sub_a_bf(c, bftmp1);
834 }
835 }
836
837 /* odd terms for sine */
838 copy_bf(bftmp2, bftmp1);
839 unsafe_mult_bf(bftmp1, bftmp2, n);
840 div_a_bf_int(bftmp1, fact++);
841 if (!sin_done)
842 {
843 sin_done = (big_accessS16(testexp) < big_accessS16(sexp)-(bflength-2)); /* too small to register */
844 if (!sin_done)
845 {
846 if (k) /* alternate between adding and subtracting */
847 unsafe_add_a_bf(s, bftmp1);
848 else
849 unsafe_sub_a_bf(s, bftmp1);
850 }
851 }
852 k = !k; /* toggle */
853 #ifdef CALCULATING_BIG_PI
854 printf("."); /* lets you know it's doing something */ 855 #endif
856 } while (!cos_done || !sin_done);
857
858 #ifndef CALCULATING_BIG_PI
859 /* now need to undo what was done by cutting angles in half */
860 for (i = 0; i < halves; i++)
861 {
862 unsafe_mult_bf(bftmp2, s, c); /* no need for safe mult */
863 double_bf(s, bftmp2); /* sin(2x) = 2*sin(x)*cos(x) */
864 unsafe_square_bf(bftmp2,c);
865 double_a_bf(bftmp2);
866 inttobf(bftmp1, 1);
867 unsafe_sub_bf(c, bftmp2, bftmp1); /* cos(2x) = 2*cos(x)*cos(x) - 1 */
868 }
869
870 if (switch_sincos)
871 {
872 copy_bf(bftmp1, s);
873 copy_bf(s, c);
874 copy_bf(c, bftmp1);
875 }
876 if (signsin)
877 neg_a_bf(s);
878 if (signcos)
879 neg_a_bf(c);
880 #endif
881
882 return s; /* return sine I guess */
883 }
884
885 /********************************************************************/
886 /* atan(r) */
887 /* uses bftmp1 - bftmp5 - global temp bigfloats */
888 /* SIDE-EFFECTS: */
889 /* n ends up as |n| or 1/|n| */
890 bf_t unsafe_atan_bf(bf_t r, bf_t n)
891 {
892 int i, comp, almost_match=0, signflag=0;
893 LDBL f;
894 bf_t orig_r, orig_n, orig_bf_pi, orig_bftmp3;
895 int orig_bflength,
896 orig_bnlength,
897 orig_padding,
898 orig_rlength,
899 orig_shiftfactor,
900 orig_rbflength;
901 int large_arg;
902
903
904 /* use Newton's recursive method for zeroing in on atan(n): r=r-cos(r)(sin(r)-n*cos(r)) */
905
906 if (is_bf_neg(n))
907 {
908 signflag = 1;
909 neg_a_bf(n);
910 }
911
912 /* If n is very large, atanl() won't give enough decimal places to be a */
913 /* good enough initial guess for Newton's Method. If it is larger than */
914 /* say, 1, atan(n) = pi/2 - acot(n) = pi/2 - atan(1/n). */
915
916 f = bftofloat(n);
917 large_arg = f > 1.0;
918 if (large_arg)
919 {
920 unsafe_inv_bf(bftmp3, n);
921 copy_bf(n, bftmp3);
922 f = bftofloat(n);
923 }
924
925 clear_bf(bftmp3); /* not really necessary, but makes things more consistent */
926
927 /* With Newton's Method, there is no need to calculate all the digits */
928 /* every time. The precision approximately doubles each iteration. */
929 /* Save original values. */
930 orig_bflength = bflength;
931 orig_bnlength = bnlength;
932 orig_padding = padding;
933 orig_rlength = rlength;
934 orig_shiftfactor = shiftfactor;
935 orig_rbflength = rbflength;
936 orig_bf_pi = bf_pi;
937 orig_r = r;
938 orig_n = n;
939 orig_bftmp3 = bftmp3;
940
941 /* calculate new starting values */
942 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
943 if (bnlength > orig_bnlength)
944 bnlength = orig_bnlength;
945 calc_lengths();
946
947 /* adjust pointers */
948 r = orig_r + orig_bflength - bflength;
949 n = orig_n + orig_bflength - bflength;
950 bf_pi = orig_bf_pi + orig_bflength - bflength;
951 bftmp3 = orig_bftmp3 + orig_bflength - bflength;
952
953 f = atanl(f); /* approximate arctangent */
954 /* no need to check overflow */
955
956 floattobf(r, f); /* start with approximate atan */
957 copy_bf(bftmp3, r);
958
959 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
960 {
961 /* adjust lengths */
962 bnlength <<= 1; /* double precision */
963 if (bnlength > orig_bnlength)
964 bnlength = orig_bnlength;
965 calc_lengths();
966 r = orig_r + orig_bflength - bflength;
967 n = orig_n + orig_bflength - bflength;
968 bf_pi = orig_bf_pi + orig_bflength - bflength;
969 bftmp3 = orig_bftmp3 + orig_bflength - bflength;
970
971 #ifdef CALCULATING_BIG_PI
973 #endif
974 unsafe_sincos_bf(bftmp4, bftmp5, bftmp3); /* sin(r), cos(r) */
975 copy_bf(bftmp3, r); /* restore bftmp3 from sincos_bf() */
976 copy_bf(bftmp1, bftmp5);
977 unsafe_mult_bf(bftmp2, n, bftmp1); /* n*cos(r) */
978 unsafe_sub_a_bf(bftmp4, bftmp2); /* sin(r) - n*cos(r) */
979 unsafe_mult_bf(bftmp1, bftmp5, bftmp4); /* cos(r) * (sin(r) - n*cos(r)) */
980 copy_bf(bftmp3, r);
981 unsafe_sub_a_bf(r, bftmp1); /* r - cos(r) * (sin(r) - n*cos(r)) */
982 #ifdef CALCULATING_BIG_PI
985 #endif
986 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp3))) < 8 ) /* if match or almost match */
987 {
988 #ifdef CALCULATING_BIG_PI
990 #endif
991 if (comp < 4 /* perfect or near perfect match */
992 || almost_match == 1) /* close enough for 2nd time */
993 break;
994 else /* this is the first time they almost matched */
995 almost_match++;
996 }
997
998 #ifdef CALCULATING_BIG_PI
1001 #endif
1002
1003 copy_bf(bftmp3, r); /* make a copy for later comparison */
1004 }
1005
1006 /* restore original values */
1007 bflength = orig_bflength;
1008 bnlength = orig_bnlength;
1009 padding = orig_padding;
1010 rlength = orig_rlength;
1011 shiftfactor = orig_shiftfactor;
1012 rbflength = orig_rbflength;
1013 bf_pi = orig_bf_pi;
1014 r = orig_r;
1015 n = orig_n;
1016 bftmp3 = orig_bftmp3;
1017
1018 if (large_arg)
1019 {
1020 half_bf(bftmp3, bf_pi); /* pi/2 */
1021 sub_a_bf(bftmp3, r); /* pi/2 - atan(1/n) */
1022 copy_bf(r, bftmp3);
1023 }
1024
1025 if (signflag)
1026 neg_a_bf(r);
1027 return r;
1028 }
1029
1030 /********************************************************************/
1031 /* atan2(r,ny,nx) */
1032 /* uses bftmp1 - bftmp6 - global temp bigfloats */
1033 bf_t unsafe_atan2_bf(bf_t r, bf_t ny, bf_t nx)
1034 {
1035 int signx, signy;
1036
1037 signx = sign_bf(nx);
1038 signy = sign_bf(ny);
1039
1040 if (signy == 0)
1041 {
1042 if (signx < 0)
1043 copy_bf(r, bf_pi); /* negative x axis, 180 deg */
1044 else /* signx >= 0 positive x axis, 0 */
1045 clear_bf(r);
1046 return(r);
1047 }
1048 if (signx == 0)
1049 {
1050 copy_bf(r, bf_pi); /* y axis */
1051 half_a_bf(r); /* +90 deg */
1052 if (signy < 0)
1053 neg_a_bf(r); /* -90 deg */
1054 return(r);
1055 }
1056
1057 if (signy < 0)
1058 neg_a_bf(ny);
1059 if (signx < 0)
1060 neg_a_bf(nx);
1061 unsafe_div_bf(bftmp6,ny,nx);
1062 unsafe_atan_bf(r, bftmp6);
1063 if (signx < 0)
1064 sub_bf(r,bf_pi,r);
1065 if(signy < 0)
1066 neg_a_bf(r);
1067 return(r);
1068 }
1069
1070 /**********************************************************************/
1071 /* The rest of the functions are "safe" versions of the routines that */
1072 /* have side effects which alter the parameters. */
1073 /* Most bf routines change values of parameters, not just the sign. */
1074 /**********************************************************************/
1075
1076 /**********************************************************************/
1077 bf_t add_bf(bf_t r, bf_t n1, bf_t n2)
1078 {
1079 copy_bf(bftmpcpy1, n1);
1080 copy_bf(bftmpcpy2, n2);
1081 unsafe_add_bf(r, bftmpcpy1, bftmpcpy2);
1082 return r;
1083 }
1084
1085 /**********************************************************************/
1086 bf_t add_a_bf(bf_t r, bf_t n)
1087 {
1088 copy_bf(bftmpcpy1, n);
1089 unsafe_add_a_bf(r, bftmpcpy1);
1090 return r;
1091 }
1092
1093 /**********************************************************************/
1094 bf_t sub_bf(bf_t r, bf_t n1, bf_t n2)
1095 {
1096 copy_bf(bftmpcpy1, n1);
1097 copy_bf(bftmpcpy2, n2);
1098 unsafe_sub_bf(r, bftmpcpy1, bftmpcpy2);
1099 return r;
1100 }
1101
1102 /**********************************************************************/
1103 bf_t sub_a_bf(bf_t r, bf_t n)
1104 {
1105 copy_bf(bftmpcpy1, n);
1106 unsafe_sub_a_bf(r, bftmpcpy1);
1107 return r;
1108 }
1109
1110 /**********************************************************************/
1111 /* mult and div only change sign */
1112 bf_t full_mult_bf(bf_t r, bf_t n1, bf_t n2)
1113 {
1114 copy_bf(bftmpcpy1, n1);
1115 copy_bf(bftmpcpy2, n2);
1116 unsafe_full_mult_bf(r, bftmpcpy1, bftmpcpy2);
1117 return r;
1118 }
1119
1120 /**********************************************************************/
1121 bf_t mult_bf(bf_t r, bf_t n1, bf_t n2)
1122 {
1123 copy_bf(bftmpcpy1, n1);
1124 copy_bf(bftmpcpy2, n2);
1125 unsafe_mult_bf(r, bftmpcpy1, bftmpcpy2);
1126 return r;
1127 }
1128
1129 /**********************************************************************/
1130 bf_t full_square_bf(bf_t r, bf_t n)
1131 {
1132 copy_bf(bftmpcpy1, n);
1133 unsafe_full_square_bf(r, bftmpcpy1);
1134 return r;
1135 }
1136
1137 /**********************************************************************/
1138 bf_t square_bf(bf_t r, bf_t n)
1139 {
1140 copy_bf(bftmpcpy1, n);
1141 unsafe_square_bf(r, bftmpcpy1);
1142 return r;
1143 }
1144
1145 /**********************************************************************/
1146 bf_t mult_bf_int(bf_t r, bf_t n, U16 u)
1147 {
1148 copy_bf(bftmpcpy1, n);
1149 unsafe_mult_bf_int(r, bftmpcpy1, u);
1150 return r;
1151 }
1152
1153 /**********************************************************************/
1154 bf_t div_bf_int(bf_t r, bf_t n, U16 u)
1155 {
1156 copy_bf(bftmpcpy1, n);
1157 unsafe_div_bf_int(r, bftmpcpy1, u);
1158 return r;
1159 }
1160
1161 /**********************************************************************/
1162 char *bftostr(char *s, int dec, bf_t r)
1163 {
1164 copy_bf(bftmpcpy1, r);
1165 unsafe_bftostr(s, dec, bftmpcpy1);
1166 return s;
1167 }
1168
1169 /**********************************************************************/
1170 char *bftostr_e(char *s, int dec, bf_t r)
1171 {
1172 copy_bf(bftmpcpy1, r);
1173 unsafe_bftostr_e(s, dec, bftmpcpy1);
1174 return s;
1175 }
1176
1177 /**********************************************************************/
1178 char *bftostr_f(char *s, int dec, bf_t r)
1179 {
1180 copy_bf(bftmpcpy1, r);
1181 unsafe_bftostr_f(s, dec, bftmpcpy1);
1182 return s;
1183 }
1184
1185 /**********************************************************************/
1186 bf_t inv_bf(bf_t r, bf_t n)
1187 {
1188 copy_bf(bftmpcpy1, n);
1189 unsafe_inv_bf(r, bftmpcpy1);
1190 return r;
1191 }
1192
1193 /**********************************************************************/
1194 bf_t div_bf(bf_t r, bf_t n1, bf_t n2)
1195 {
1196 copy_bf(bftmpcpy1, n1);
1197 copy_bf(bftmpcpy2, n2);
1198 unsafe_div_bf(r, bftmpcpy1, bftmpcpy2);
1199 return r;
1200 }
1201
1202 /**********************************************************************/
1203 bf_t sqrt_bf(bf_t r, bf_t n)
1204 {
1205 copy_bf(bftmpcpy1, n);
1206 unsafe_sqrt_bf(r, bftmpcpy1);
1207 return r;
1208 }
1209
1210 /**********************************************************************/
1211 bf_t ln_bf(bf_t r, bf_t n)
1212 {
1213 copy_bf(bftmpcpy1, n);
1214 unsafe_ln_bf(r, bftmpcpy1);
1215 return r;
1216 }
1217
1218 /**********************************************************************/
1219 bf_t sincos_bf(bf_t s, bf_t c, bf_t n)
1220 {
1221 copy_bf(bftmpcpy1, n);
1222 return unsafe_sincos_bf(s, c, bftmpcpy1);
1223 }
1224
1225 /**********************************************************************/
1226 bf_t atan_bf(bf_t r, bf_t n)
1227 {
1228 copy_bf(bftmpcpy1, n);
1229 unsafe_atan_bf(r, bftmpcpy1);
1230 return r;
1231 }
1232
1233 /**********************************************************************/
1234 bf_t atan2_bf(bf_t r, bf_t ny, bf_t nx)
1235 {
1236 copy_bf(bftmpcpy1, ny);
1237 copy_bf(bftmpcpy2, nx);
1238 unsafe_atan2_bf(r, bftmpcpy1, bftmpcpy2);
1239 return r;
1240 }
1241
1242 /**********************************************************************/
1243 int is_bf_zero(bf_t n)
1244 {
1245 return !is_bf_not_zero(n);
1246 }
1247
1248 /************************************************************************/
1249 /* convert_bf -- convert bigfloat numbers from old to new lengths */
1250 int convert_bf(bf_t new, bf_t old, int newbflength, int oldbflength)
1251 {
1252 int savebflength;
1253
1254 /* save lengths so not dependent on external environment */
1255 savebflength = bflength;
1256 bflength = newbflength;
1257 clear_bf(new);
1258 bflength = savebflength;
1259
1260 if(newbflength > oldbflength)
1261 _fmemcpy(new+newbflength-oldbflength,old,oldbflength+2);
1262 else
1263 _fmemcpy(new,old+oldbflength-newbflength,newbflength+2);
1264 return(0);
1265 }
1266
1267 /* The following used to be in bigfltc.c */
1268 /********************************************************************/
1269 /* normalize big float */
1270 bf_t norm_bf(bf_t r)
1271 {
1272 int scale;
1273 BYTE hi_byte;
1274 S16 BIGDIST *rexp;
1275
1276 rexp = (S16 BIGDIST *)(r+bflength);
1277
1278 /* check for overflow */
1279 hi_byte = r[bflength-1];
1280 if (hi_byte != 0x00 && hi_byte != 0xFF)
1281 {
1282 _fmemmove(r, r+1, bflength-1);
1283 r[bflength-1] = (BYTE)(hi_byte & 0x80 ? 0xFF : 0x00);
1284 big_setS16(rexp,big_accessS16(rexp)+(S16)1); /* exp */
1285 }
1286
1287 /* check for underflow */
1288 else
1289 {
1290 for (scale = 2; scale < bflength && r[bflength-scale] == hi_byte; scale++)
1291 ; /* do nothing */
1292 if (scale == bflength && hi_byte == 0) /* zero */
1293 big_setS16(rexp,0);
1294 else
1295 {
1296 scale -= 2;
1297 if (scale > 0) /* it did underflow */
1298 {
1299 _fmemmove(r+scale, r, bflength-scale-1);
1300 _fmemset(r, 0, scale);
1301 big_setS16(rexp,big_accessS16(rexp)-(S16)scale); /* exp */
1302 }
1303 }
1304 }
1305
1306 return r;
1307 }
1308
1309 /********************************************************************/
1310 /* normalize big float with forced sign */
1311 /* positive = 1, force to be positive */
1312 /* = 0, force to be negative */
1313 void norm_sign_bf(bf_t r, int positive)
1314 {
1315 norm_bf(r);
1316 r[bflength-1] = (BYTE)(positive ? 0x00 : 0xFF);
1317 }
1318 /******************************************************/
1319 /* adjust n1, n2 for before addition or subtraction */
1320 /* by forcing exp's to match. */
1321 /* returns the value of the adjusted exponents */
1322 S16 adjust_bf_add(bf_t n1, bf_t n2)
1323 {
1324 int scale, fill_byte;
1325 S16 rexp;
1326 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1327
1328 /* scale n1 or n2 */
1329 /* compare exp's */
1330 n1exp = (S16 BIGDIST *)(n1+bflength);
1331 n2exp = (S16 BIGDIST *)(n2+bflength);
1332 if (big_accessS16(n1exp) > big_accessS16(n2exp))
1333 { /* scale n2 */
1334 scale = big_accessS16(n1exp) - big_accessS16(n2exp); /* n1exp - n2exp */
1335 if (scale < bflength)
1336 {
1337 fill_byte = is_bf_neg(n2) ? 0xFF : 0x00;
1338 _fmemmove(n2, n2+scale, bflength-scale);
1339 _fmemset(n2+bflength-scale, fill_byte, scale);
1340 }
1341 else
1342 clear_bf(n2);
1343 big_setS16(n2exp,big_accessS16(n1exp)); /* *n2exp = *n1exp; set exp's = */
1344 rexp = big_accessS16(n2exp);
1345 }
1346 else if (big_accessS16(n1exp) < big_accessS16(n2exp))
1347 { /* scale n1 */
1348 scale = big_accessS16(n2exp) - big_accessS16(n1exp); /* n2exp - n1exp */
1349 if (scale < bflength)
1350 {
1351 fill_byte = is_bf_neg(n1) ? 0xFF : 0x00;
1352 _fmemmove(n1, n1+scale, bflength-scale);
1353 _fmemset(n1+bflength-scale, fill_byte, scale);
1354 }
1355 else
1356 clear_bf(n1);
1357 big_setS16(n1exp,big_accessS16(n2exp)); /* *n1exp = *n2exp; set exp's = */
1358 rexp = big_accessS16(n2exp);
1359 }
1360 else
1361 rexp = big_accessS16(n1exp);
1362 return rexp;
1363 }
1364
1365 /********************************************************************/
1366 /* r = max positive value */
1367 bf_t max_bf(bf_t r)
1368 {
1369 inttobf(r, 1);
1370 big_set16(r+bflength, (S16)(LDBL_MAX_EXP/8));
1371 return r;
1372 }
1373
1374 /****************************************************************************/
1375 /* n1 != n2 ? */
1376 /* RETURNS: */
1377 /* if n1 == n2 returns 0 */
1378 /* if n1 > n2 returns a positive (bytes left to go when mismatch occurred) */
1379 /* if n1 < n2 returns a negative (bytes left to go when mismatch occurred) */
1380
1381 int cmp_bf(bf_t n1, bf_t n2)
1382 {
1383 int i;
1384 int sign1, sign2;
1385 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1386 U16 value1, value2;
1387
1388
1389 /* compare signs */
1390 sign1 = sign_bf(n1);
1391 sign2 = sign_bf(n2);
1392 if (sign1 > sign2)
1393 return bflength;
1394 else if (sign1 < sign2)
1395 return -bflength;
1396 /* signs are the same */
1397
1398 /* compare exponents, using signed comparisons */
1399 n1exp = (S16 BIGDIST *)(n1+bflength);
1400 n2exp = (S16 BIGDIST *)(n2+bflength);
1401 if (big_accessS16(n1exp) > big_accessS16(n2exp))
1402 return sign1*(bflength);
1403 else if (big_accessS16(n1exp) < big_accessS16(n2exp))
1404 return -sign1*(bflength);
1405
1406 /* To get to this point, the signs must match */
1407 /* so unsigned comparison is ok. */
1408 /* two bytes at a time */
1409 for (i=bflength-2; i>=0; i-=2)
1410 {
1411 if ( (value1=big_access16(n1+i)) > (value2=big_access16(n2+i)) )
1412 { /* now determine which of the two bytes was different */
1413 if ( (value1&0xFF00) > (value2&0xFF00) ) /* compare just high bytes */
1414 return (i+2); /* high byte was different */
1415 else
1416 return (i+1); /* low byte was different */
1417 }
1418 else if (value1 < value2)
1419 { /* now determine which of the two bytes was different */
1420 if ( (value1&0xFF00) < (value2&0xFF00) ) /* compare just high bytes */
1421 return -(i+2); /* high byte was different */
1422 else
1423 return -(i+1); /* low byte was different */
1424 }
1425 }
1426 return 0;
1427 }
1428
1429 /********************************************************************/
1430 /* r < 0 ? */
1431 /* returns 1 if negative, 0 if positive or zero */
1432 int is_bf_neg(bf_t n)
1433 {
1434 return (S8)n[bflength-1] < 0;
1435 }
1436
1437 /********************************************************************/
1438 /* n != 0 ? */
1439 /* RETURNS: if n != 0 returns 1 */
1440 /* else returns 0 */
1441 int is_bf_not_zero(bf_t n)
1442 {
1443 int bnl;
1444 int retval;
1445
1446 bnl = bnlength;
1447 bnlength = bflength;
1448 retval = is_bn_not_zero(n);
1449 bnlength = bnl;
1450
1451 return retval;
1452 }
1453
1454 /********************************************************************/
1455 /* r = n1 + n2 */
1456 /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
1457 bf_t unsafe_add_bf(bf_t r, bf_t n1, bf_t n2)
1458 {
1459 int bnl;
1460 S16 BIGDIST *rexp;
1461
1462 if (is_bf_zero(n1))
1463 {
1464 copy_bf(r, n2);
1465 return r;
1466 }
1467 if (is_bf_zero(n2))
1468 {
1469 copy_bf(r, n1);
1470 return r;
1471 }
1472
1473 rexp = (S16 BIGDIST *)(r+bflength);
1474 big_setS16(rexp,adjust_bf_add(n1, n2));
1475
1476 bnl = bnlength;
1477 bnlength = bflength;
1478 add_bn(r, n1, n2);
1479 bnlength = bnl;
1480
1481 norm_bf(r);
1482 return r;
1483 }
1484
1485 /********************************************************************/
1486 /* r += n */
1487 bf_t unsafe_add_a_bf(bf_t r, bf_t n)
1488 {
1489 int bnl;
1490
1491 if (is_bf_zero(r))
1492 {
1493 copy_bf(r, n);
1494 return r;
1495 }
1496 if (is_bf_zero(n))
1497 {
1498 return r;
1499 }
1500
1501 adjust_bf_add(r, n);
1502
1503 bnl = bnlength;
1504 bnlength = bflength;
1505 add_a_bn(r, n);
1506 bnlength = bnl;
1507
1508 norm_bf(r);
1509
1510 return r;
1511 }
1512
1513 /********************************************************************/
1514 /* r = n1 - n2 */
1515 /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
1516 bf_t unsafe_sub_bf(bf_t r, bf_t n1, bf_t n2)
1517 {
1518 int bnl;
1519 S16 BIGDIST *rexp;
1520
1521 if (is_bf_zero(n1))
1522 {
1523 neg_bf(r, n2);
1524 return r;
1525 }
1526 if (is_bf_zero(n2))
1527 {
1528 copy_bf(r, n1);
1529 return r;
1530 }
1531
1532 rexp = (S16 BIGDIST *)(r+bflength);
1533 big_setS16(rexp,adjust_bf_add(n1, n2));
1534
1535 bnl = bnlength;
1536 bnlength = bflength;
1537 sub_bn(r, n1, n2);
1538 bnlength = bnl;
1539
1540 norm_bf(r);
1541 return r;
1542 }
1543
1544 /********************************************************************/
1545 /* r -= n */
1546 bf_t unsafe_sub_a_bf(bf_t r, bf_t n)
1547 {
1548 int bnl;
1549
1550 if (is_bf_zero(r))
1551 {
1552 neg_bf(r,n);
1553 return r;
1554 }
1555 if (is_bf_zero(n))
1556 {
1557 return r;
1558 }
1559
1560 adjust_bf_add(r, n);
1561
1562 bnl = bnlength;
1563 bnlength = bflength;
1564 sub_a_bn(r, n);
1565 bnlength = bnl;
1566
1567 norm_bf(r);
1568 return r;
1569 }
1570
1571 /********************************************************************/
1572 /* r = -n */
1573 bf_t neg_bf(bf_t r, bf_t n)
1574 {
1575 int bnl;
1576 S16 BIGDIST *rexp, BIGDIST *nexp;
1577
1578 rexp = (S16 BIGDIST *)(r+bflength);
1579 nexp = (S16 BIGDIST *)(n+bflength);
1580 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1581
1582 bnl = bnlength;
1583 bnlength = bflength;
1584 neg_bn(r, n);
1585 bnlength = bnl;
1586
1587 norm_bf(r);
1588 return r;
1589 }
1590
1591 /********************************************************************/
1592 /* r *= -1 */
1593 bf_t neg_a_bf(bf_t r)
1594 {
1595 int bnl;
1596
1597 bnl = bnlength;
1598 bnlength = bflength;
1599 neg_a_bn(r);
1600 bnlength = bnl;
1601
1602 norm_bf(r);
1603 return r;
1604 }
1605
1606 /********************************************************************/
1607 /* r = 2*n */
1608 bf_t double_bf(bf_t r, bf_t n)
1609 {
1610 int bnl;
1611 S16 BIGDIST *rexp, BIGDIST *nexp;
1612
1613 rexp = (S16 BIGDIST *)(r+bflength);
1614 nexp = (S16 BIGDIST *)(n+bflength);
1615 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1616
1617 bnl = bnlength;
1618 bnlength = bflength;
1619 double_bn(r, n);
1620 bnlength = bnl;
1621
1622 norm_bf(r);
1623 return r;
1624 }
1625
1626 /********************************************************************/
1627 /* r *= 2 */
1628 bf_t double_a_bf(bf_t r)
1629 {
1630 int bnl;
1631
1632 bnl = bnlength;
1633 bnlength = bflength;
1634 double_a_bn(r);
1635 bnlength = bnl;
1636
1637 norm_bf(r);
1638 return r;
1639 }
1640
1641 /********************************************************************/
1642 /* r = n/2 */
1643 bf_t half_bf(bf_t r, bf_t n)
1644 {
1645 int bnl;
1646 S16 BIGDIST *rexp, BIGDIST *nexp;
1647
1648 rexp = (S16 BIGDIST *)(r+bflength);
1649 nexp = (S16 BIGDIST *)(n+bflength);
1650 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1651
1652 bnl = bnlength;
1653 bnlength = bflength;
1654 half_bn(r, n);
1655 bnlength = bnl;
1656
1657 norm_bf(r);
1658 return r;
1659 }
1660
1661 /********************************************************************/
1662 /* r /= 2 */
1663 bf_t half_a_bf(bf_t r)
1664 {
1665 int bnl;
1666
1667 bnl = bnlength;
1668 bnlength = bflength;
1669 half_a_bn(r);
1670 bnlength = bnl;
1671
1672 norm_bf(r);
1673 return r;
1674 }
1675
1676 /************************************************************************/
1677 /* r = n1 * n2 */
1678 /* Note: r will be a double wide result, 2*bflength */
1679 /* n1 and n2 can be the same pointer */
1680 /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values */
1681 bf_t unsafe_full_mult_bf(bf_t r, bf_t n1, bf_t n2)
1682 {
1683 int bnl, dbfl;
1684 S16 BIGDIST *rexp, BIGDIST *n1exp, BIGDIST *n2exp;
1685
1686 if (is_bf_zero(n1) || is_bf_zero(n2) )
1687 {
1688 bflength <<= 1;
1689 clear_bf(r);
1690 bflength >>= 1;
1691 return r;
1692 }
1693
1694 dbfl = 2*bflength; /* double width bflength */
1695 rexp = (S16 BIGDIST *)(r+dbfl); /* note: 2*bflength */
1696 n1exp = (S16 BIGDIST *)(n1+bflength);
1697 n2exp = (S16 BIGDIST *)(n2+bflength);
1698 /* add exp's */
1699 big_setS16(rexp, (S16)(big_accessS16(n1exp) + big_accessS16(n2exp)) );
1700
1701 bnl = bnlength;
1702 bnlength = bflength;
1703 unsafe_full_mult_bn(r, n1, n2);
1704 bnlength = bnl;
1705
1706 /* handle normalizing full mult on individual basis */
1707
1708 return r;
1709 }
1710
1711 /************************************************************************/
1712 /* r = n1 * n2 calculating only the top rlength bytes */
1713 /* Note: r will be of length rlength */
1714 /* 2*bflength <= rlength < bflength */
1715 /* n1 and n2 can be the same pointer */
1716 /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values */
1717 bf_t unsafe_mult_bf(bf_t r, bf_t n1, bf_t n2)
1718 {
1719 int positive;
1720 int bnl, bfl, rl;
1721 int rexp;
1722 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1723
1724 if (is_bf_zero(n1) || is_bf_zero(n2) )
1725 {
1726 clear_bf(r);
1727 return r;
1728 }
1729
1730 n1exp = (S16 BIGDIST *)(n1+bflength);
1731 n2exp = (S16 BIGDIST *)(n2+bflength);
1732 /* add exp's */
1733 rexp = big_accessS16(n1exp) + big_accessS16(n2exp);
1734
1735 positive = is_bf_neg(n1) == is_bf_neg(n2); /* are they the same sign? */
1736
1737 bnl = bnlength;
1738 bnlength = bflength;
1739 rl = rlength;
1740 rlength = rbflength;
1741 unsafe_mult_bn(r, n1, n2);
1742 bnlength = bnl;
1743 rlength = rl;
1744
1745 bfl = bflength;
1746 bflength = rbflength;
1747 big_set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
1748 norm_sign_bf(r, positive);
1749 bflength = bfl;
1750 _fmemmove(r, r+padding, bflength+2); /* shift back */
1751
1752 return r;
1753 }
1754
1755 /************************************************************************/
1756 /* r = n^2 */
1757 /* because of the symmetry involved, n^2 is much faster than n*n */
1758 /* for a bignumber of length l */
1759 /* n*n takes l^2 multiplications */
1760 /* n^2 takes (l^2+l)/2 multiplications */
1761 /* which is about 1/2 n*n as l gets large */
1762 /* uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
1763 /* */
1764 /* SIDE-EFFECTS: n is changed to its absolute value */
1765 bf_t unsafe_full_square_bf(bf_t r, bf_t n)
1766 {
1767 int bnl, dbfl;
1768 S16 BIGDIST *rexp, BIGDIST *nexp;
1769
1770 if (is_bf_zero(n))
1771 {
1772 bflength <<= 1;
1773 clear_bf(r);
1774 bflength >>= 1;
1775 return r;
1776 }
1777
1778 dbfl = 2*bflength; /* double width bflength */
1779 rexp = (S16 BIGDIST *)(r+dbfl); /* note: 2*bflength */
1780 nexp = (S16 BIGDIST *)(n+bflength);
1781 big_setS16(rexp, 2 * big_accessS16(nexp));
1782
1783 bnl = bnlength;
1784 bnlength = bflength;
1785 unsafe_full_square_bn(r, n);
1786 bnlength = bnl;
1787
1788 /* handle normalizing full mult on individual basis */
1789
1790 return r;
1791 }
1792
1793
1794 /************************************************************************/
1795 /* r = n^2 */
1796 /* because of the symmetry involved, n^2 is much faster than n*n */
1797 /* for a bignumber of length l */
1798 /* n*n takes l^2 multiplications */
1799 /* n^2 takes (l^2+l)/2 multiplications */
1800 /* which is about 1/2 n*n as l gets large */
1801 /* uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
1802 /* */
1803 /* Note: r will be of length rlength */
1804 /* 2*bflength >= rlength > bflength */
1805 /* SIDE-EFFECTS: n is changed to its absolute value */
1806 bf_t unsafe_square_bf(bf_t r, bf_t n)
1807 {
1808 int bnl, bfl, rl;
1809 int rexp;
1810 S16 BIGDIST *nexp;
1811
1812 if (is_bf_zero(n))
1813 {
1814 clear_bf(r);
1815 return r;
1816 }
1817
1818 nexp = (S16 BIGDIST *)(n+bflength);
1819 rexp = (S16)(2 * big_accessS16(nexp));
1820
1821 bnl = bnlength;
1822 bnlength = bflength;
1823 rl = rlength;
1824 rlength = rbflength;
1825 unsafe_square_bn(r, n);
1826 bnlength = bnl;
1827 rlength = rl;
1828
1829 bfl = bflength;
1830 bflength = rbflength;
1831 big_set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
1832
1833 norm_sign_bf(r, 1);
1834 bflength = bfl;
1835 _fmemmove(r, r+padding, bflength+2); /* shift back */
1836
1837 return r;
1838 }
1839
1840 /********************************************************************/
1841 /* r = n * u where u is an unsigned integer */
1842 /* SIDE-EFFECTS: n can be "de-normalized" and lose precision */
1843 bf_t unsafe_mult_bf_int(bf_t r, bf_t n, U16 u)
1844 {
1845 int positive;
1846 int bnl;
1847 S16 BIGDIST *rexp, BIGDIST *nexp;
1848
1849 rexp = (S16 BIGDIST *)(r+bflength);
1850 nexp = (S16 BIGDIST *)(n+bflength);
1851 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1852
1853 positive = !is_bf_neg(n);
1854
1855 /*
1856 if u > 0x00FF, then the integer part of the mantissa will overflow the
1857 2 byte (16 bit) integer size. Therefore, make adjustment before
1858 multiplication is performed.
1859 */
1860 if (u > 0x00FF)
1861 { /* un-normalize n */
1862 _fmemmove(n, n+1, bflength-1); /* this sign extends as well */
1863 big_setS16(rexp,big_accessS16(rexp)+(S16)1);
1864 }
1865
1866 bnl = bnlength;
1867 bnlength = bflength;
1868 mult_bn_int(r, n, u);
1869 bnlength = bnl;
1870
1871 norm_sign_bf(r, positive);
1872 return r;
1873 }
1874
1875 /********************************************************************/
1876 /* r *= u where u is an unsigned integer */
1877 bf_t mult_a_bf_int(bf_t r, U16 u)
1878 {
1879 int positive;
1880 int bnl;
1881 S16 BIGDIST *rexp;
1882
1883 rexp = (S16 BIGDIST *)(r+bflength);
1884 positive = !is_bf_neg(r);
1885
1886 /*
1887 if u > 0x00FF, then the integer part of the mantissa will overflow the
1888 2 byte (16 bit) integer size. Therefore, make adjustment before
1889 multiplication is performed.
1890 */
1891 if (u > 0x00FF)
1892 { /* un-normalize n */
1893 _fmemmove(r, r+1, bflength-1); /* this sign extends as well */
1894 big_setS16(rexp,big_accessS16(rexp)+(S16)1);
1895 }
1896
1897 bnl = bnlength;
1898 bnlength = bflength;
1899 mult_a_bn_int(r, u);
1900 bnlength = bnl;
1901
1902 norm_sign_bf(r, positive);
1903 return r;
1904 }
1905
1906 /********************************************************************/
1907 /* r = n / u where u is an unsigned integer */
1908 bf_t unsafe_div_bf_int(bf_t r, bf_t n, U16 u)
1909 {
1910 int bnl;
1911 S16 BIGDIST *rexp, BIGDIST *nexp;
1912
1913 if (u == 0) /* division by zero */
1914 {
1915 max_bf(r);
1916 if (is_bf_neg(n))
1917 neg_a_bf(r);
1918 return r;
1919 }
1920
1921 rexp = (S16 BIGDIST *)(r+bflength);
1922 nexp = (S16 BIGDIST *)(n+bflength);
1923 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1924
1925 bnl = bnlength;
1926 bnlength = bflength;
1927 unsafe_div_bn_int(r, n, u);
1928 bnlength = bnl;
1929
1930 norm_bf(r);
1931 return r;
1932 }
1933
1934 /********************************************************************/
1935 /* r /= u where u is an unsigned integer */
1936 bf_t div_a_bf_int(bf_t r, U16 u)
1937 {
1938 int bnl;
1939
1940 if (u == 0) /* division by zero */
1941 {
1942 if (is_bf_neg(r))
1943 {
1944 max_bf(r);
1945 neg_a_bf(r);
1946 }
1947 else
1948 {
1949 max_bf(r);
1950 }
1951 return r;
1952 }
1953
1954 bnl = bnlength;
1955 bnlength = bflength;
1956 div_a_bn_int(r, u);
1957 bnlength = bnl;
1958
1959 norm_bf(r);
1960 return r;
1961 }
1962
1963 /********************************************************************/
1964 /* extracts the mantissa and exponent of f */
1965 /* finds m and n such that 1<=|m|<b and f = m*b^n */
1966 /* n is stored in *exp_ptr and m is returned, sort of like frexp() */
1967 LDBL extract_value(LDBL f, LDBL b, int *exp_ptr)
1968 {
1969 int n;
1970 LDBL af, ff, orig_b;
1971 LDBL value[15];
1972 unsigned powertwo;
1973
1974 if (b <= 0 || f == 0)
1975 {
1976 *exp_ptr = 0;
1977 return 0;
1978 }
1979
1980 orig_b = b;
1981 af = f >= 0 ? f: -f; /* abs value */
1982 ff = af > 1 ? af : 1/af;
1983 n = 0;
1984 powertwo = 1;
1985 while (b < ff)
1986 {
1987 value[n] = b;
1988 n++;
1989 powertwo <<= 1;
1990 b *= b;
1991 }
1992
1993 *exp_ptr = 0;
1994 for (; n > 0; n--)
1995 {
1996 powertwo >>= 1;
1997 if (value[n-1] < ff)
1998 {
1999 ff /= value[n-1];
2000 *exp_ptr += powertwo;
2001 }
2002 }
2003 if (f < 0)
2004 ff = -ff;
2005 if (af < 1)
2006 {
2007 ff = orig_b/ff;
2008 *exp_ptr = -*exp_ptr - 1;
2009 }
2010
2011 return ff;
2012 }
2013
2014 /********************************************************************/
2015 /* calculates and returns the value of f*b^n */
2016 /* sort of like ldexp() */
2017 LDBL scale_value( LDBL f, LDBL b , int n )
2018 {
2019 LDBL total=1;
2020 int an;
2021
2022 if (b == 0 || f == 0)
2023 return 0;
2024
2025 if (n == 0)
2026 return f;
2027
2028 an = abs(n);
2029
2030 while (an != 0)
2031 {
2032 if (an & 0x0001)
2033 total *= b;
2034 b *= b;
2035 an >>= 1;
2036 }
2037
2038 if (n > 0)
2039 f *= total;
2040 else /* n < 0 */
2041 f /= total;
2042 return f;
2043 }
2044
2045 /********************************************************************/
2046 /* extracts the mantissa and exponent of f */
2047 /* finds m and n such that 1<=|m|<10 and f = m*10^n */
2048 /* n is stored in *exp_ptr and m is returned, sort of like frexp() */
2049 LDBL extract_10(LDBL f, int *exp_ptr)
2050 {
2051 return extract_value(f, 10, exp_ptr);
2052 }
2053
2054 /********************************************************************/
2055 /* calculates and returns the value of f*10^n */
2056 /* sort of like ldexp() */
2057 LDBL scale_10( LDBL f, int n )
2058 {
2059 return scale_value( f, 10, n );
2060 }
2061
2062
2063
2064 /* big10flt.c - C routines for base 10 big floating point numbers */
2065
2066 /**********************************************************
2067 (Just when you thought it was safe to go back in the water.)
2068 Just when you thought you seen every type of format possible,
2069 16 bit integer, 32 bit integer, double, long double, mpmath,
2070 bn_t, bf_t, I now give you bf10_t (big float base 10)!
2071
2072 Why, because this is the only way (I can think of) to properly do a
2073 bftostr() without rounding errors. With out this, then
2074 -1.9999999999( > LDBL_DIG of 9's)9999999123456789...
2075 will round to -2.0. The good news is that we only need to do two
2076 mathematical operations: multiplication and division by integers
2077
2078 bf10_t format: (notice the position of the MSB and LSB)
2079
2080 MSB LSB
2081 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2082 n <><------------- dec --------------><> <->
2083 1 byte pad 1 byte rounding 2 byte exponent.
2084
2085 total length = dec + 4
2086
2087 ***********************************************************/
2088
2089 /**********************************************************************/
2090 /* unsafe_bftobf10() - converts a bigfloat into a bigfloat10 */
2091 /* n - pointer to a bigfloat */
2092 /* r - result array of BYTE big enough to hold the bf10_t number */
2093 /* dec - number of decimals, not including the one extra for rounding */
2094 /* SIDE-EFFECTS: n is changed to |n|. Make copy of n if necessary. */
2095
2096 bf10_t unsafe_bftobf10(bf10_t r, int dec, bf_t n)
2097 {
2098 int d;
2099 int power256;
2100 int p;
2101 int bnl;
2102 bf_t onesbyte;
2103 bf10_t power10;
2104
2105 if (is_bf_zero(n))
2106 { /* in scientific notation, the leading digit can't be zero */
2107 r[1] = (BYTE)0; /* unless the number is zero */
2108 return r;
2109 }
2110
2111 onesbyte = n + bflength - 1; /* really it's n+bflength-2 */
2112 power256 = (S16)big_access16(n + bflength) + 1; /* so adjust power256 by 1 */
2113
2114 if (dec == 0)
2115 dec = decimals;
2116 dec++; /* one extra byte for rounding */
2117 power10 = r + dec + 1;
2118
2119 if (is_bf_neg(n))
2120 {
2121 neg_a_bf(n);
2122 r[0] = 1; /* sign flag */
2123 }
2124 else
2125 r[0] = 0;
2126
2127 p = -1; /* multiply by 10 right away */
2128 bnl = bnlength;
2129 bnlength = bflength;
2130 for (d=1; d<=dec; d++)
2131 {
2132 /* pretend it's a bn_t instead of a bf_t */
2133 /* this leaves n un-normalized, which is what we want here */
2134 mult_a_bn_int(n, 10);
2135
2136 r[d] = *onesbyte;
2137 if (d == 1 && r[d] == 0)
2138 {
2139 d = 0; /* back up a digit */
2140 p--; /* and decrease by a factor of 10 */
2141 }
2142 *onesbyte = 0;
2143 }
2144 bnlength = bnl;
2145 big_set16(power10, (U16)p); /* save power of ten */
2146
2147 /* the digits are all read in, now scale it by 256^power256 */
2148 if (power256 > 0)
2149 for (d=0; d<power256; d++)
2150 mult_a_bf10_int(r, dec, 256);
2151
2152 else if (power256 < 0)
2153 for (d=0; d>power256; d--)
2154 div_a_bf10_int(r, dec, 256);
2155
2156 /* else power256 is zero, don't do anything */
2157
2158 /* round the last digit */
2159 if (r[dec] >= 5)
2160 {
2161 d = dec-1;
2162 while (d > 0) /* stop before you get to the sign flag */
2163 {
2164 r[d]++; /* round up */
2165 if (r[d] < 10)
2166 {
2167 d = -1; /* flag for below */
2168 break; /* finished rounding */
2169 }
2170 r[d] = 0;
2171 d--;
2172 }
2173 if (d == 0) /* rounding went back to the first digit and it overflowed */
2174 {
2175 r[1] = 0;
2176 _fmemmove(r+2, r+1, dec-1);
2177 r[1] = 1;
2178 p = (S16)big_access16(power10);
2179 big_set16(power10, (U16)(p+1));
2180 }
2181 }
2182 r[dec] = 0; /* truncate the rounded digit */
2183
2184 return r;
2185 }
2186
2187
2188 /**********************************************************************/
2189 /* mult_a_bf10_int() */
2190 /* r *= n */
2191 /* dec - number of decimals, including the one extra for rounding */
2192
2193 bf10_t mult_a_bf10_int(bf10_t r, int dec, U16 n)
2194 {
2195 int signflag;
2196 int d, p;
2197 unsigned value, overflow;
2198 bf10_t power10;
2199
2200 if (r[1] == 0 || n == 0)
2201 {
2202 r[1] = 0;
2203 return r;
2204 }
2205
2206 power10 = r + dec + 1;
2207 p = (S16)big_access16(power10);
2208
2209 signflag = r[0]; /* r[0] to be used as a padding */
2210 overflow = 0;
2211 for (d = dec; d>0; d--)
2212 {
2213 value = r[d] * n + overflow;
2214 r[d] = (BYTE)(value % 10);
2215 overflow = value / 10;
2216 }
2217 while (overflow)
2218 {
2219 p++;
2220 _fmemmove(r+2, r+1, dec-1);
2221 r[1] = (BYTE)(overflow % 10);
2222 overflow = overflow / 10;
2223 }
2224 big_set16(power10, (U16)p); /* save power of ten */
2225 r[0] = (BYTE)signflag; /* restore sign flag */
2226 return r;
2227 }
2228
2229 /**********************************************************************/
2230 /* div_a_bf10_int() */
2231 /* r /= n */
2232 /* dec - number of decimals, including the one extra for rounding */
2233
2234 bf10_t div_a_bf10_int (bf10_t r, int dec, U16 n)
2235 {
2236 int src, dest, p;
2237 unsigned value, remainder;
2238 bf10_t power10;
2239
2240 if (r[1] == 0 || n == 0)
2241 {
2242 r[1] = 0;
2243 return r;
2244 }
2245
2246 power10 = r + dec + 1;
2247 p = (S16)big_access16(power10);
2248
2249 remainder = 0;
2250 for (src=dest=1; src<=dec; dest++, src++)
2251 {
2252 value = 10*remainder + r[src];
2253 r[dest] = (BYTE)(value / n);
2254 remainder = value % n;
2255 if (dest == 1 && r[dest] == 0)
2256 {
2257 dest = 0; /* back up a digit */
2258 p--; /* and decrease by a factor of 10 */
2259 }
2260 }
2261 for (; dest<=dec; dest++)
2262 {
2263 value = 10*remainder;
2264 r[dest] = (BYTE)(value / n);
2265 remainder = value % n;
2266 if (dest == 1 && r[dest] == 0)
2267 {
2268 dest = 0; /* back up a digit */
2269 p--; /* and decrease by a factor of 10 */
2270 }
2271 }
2272
2273 big_set16(power10, (U16)p); /* save power of ten */
2274 return r;
2275 }
2276
2277
2278 /*************************************************************************/
2279 /* bf10tostr_e() */
2280 /* Takes a bf10 number and converts it to an ascii string, sci. notation */
2281 /* dec - number of decimals, not including the one extra for rounding */
2282
2283 char *bf10tostr_e(char *s, int dec, bf10_t n)
2284 {
2285 int d, p;
2286 bf10_t power10;
2287
2288 if (n[1] == 0)
2289 {
2290 strcpy(s, "0.0");
2291 return s;
2292 }
2293
2294 if (dec == 0)
2295 dec = decimals;
2296 dec++; /* one extra byte for rounding */
2297 power10 = n + dec + 1;
2298 p = (S16)big_access16(power10);
2299
2300 /* if p is negative, it is not necessary to show all the decimal places */
2301 if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
2302 {
2303 dec = dec + p;
2304 if (dec < 8) /* let's keep at least a few */
2305 dec = 8;
2306 }
2307
2308 if (n[0] == 1) /* sign flag */
2309 *(s++) = '-';
2310 *(s++) = (char)(n[1] + '0');
2311 *(s++) = '.';
2312 for (d=2; d<=dec; d++)
2313 {
2314 *(s++) = (char)(n[d] + '0');
2315 }
2316 /* clean up trailing 0's */
2317 while (*(s-1) == '0')
2318 s--;
2319 if (*(s-1) == '.') /* put at least one 0 after the decimal */
2320 *(s++) = '0';
2321 sprintf(s, "e%d", p);
2322 return s;
2323 }
2324
2325 /****************************************************************************/
2326 /* bf10tostr_f() */
2327 /* Takes a bf10 number and converts it to an ascii string, decimal notation */
2328
2329 char *bf10tostr_f(char *s, int dec, bf10_t n)
2330 {
2331 int d, p;
2332 bf10_t power10;
2333
2334 if (n[1] == 0)
2335 {
2336 strcpy(s, "0.0");
2337 return s;
2338 }
2339
2340 if (dec == 0)
2341 dec = decimals;
2342 dec++; /* one extra byte for rounding */
2343 power10 = n + dec + 1;
2344 p = (S16)big_access16(power10);
2345
2346 /* if p is negative, it is not necessary to show all the decimal places */
2347 if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
2348 {
2349 dec = dec + p;
2350 if (dec < 8) /* let's keep at least a few */
2351 dec = 8;
2352 }
2353
2354 if (n[0] == 1) /* sign flag */
2355 *(s++) = '-';
2356 if (p >= 0)
2357 {
2358 for (d=1; d<=p+1; d++)
2359 *(s++) = (char)(n[d] + '0');
2360 *(s++) = '.';
2361 for (; d<=dec; d++)
2362 *(s++) = (char)(n[d] + '0');
2363 }
2364 else
2365 {
2366 *(s++) = '0';
2367 *(s++) = '.';
2368 for (d=0; d>p+1; d--)
2369 *(s++) = '0';
2370 for (d=1; d<=dec; d++)
2371 *(s++) = (char)(n[d] + '0');
2372 }
2373
2374 /* clean up trailing 0's */
2375 while (*(s-1) == '0')
2376 s--;
2377 if (*(s-1) == '.') /* put at least one 0 after the decimal */
2378 *(s++) = '0';
2379 *s = '\0'; /* terminating nul */
2380 return s;
2381 }
2382