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
bf_hexdump(bf_t r)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
strtobf(bf_t r,char * s)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
strlen_needed_bf()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
unsafe_bftostr(char * s,int dec,bf_t r)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) */
unsafe_bftostr_e(char * s,int dec,bf_t r)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) */
unsafe_bftostr_f(char * s,int dec,bf_t r)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 */
bftobn(bn_t n,bf_t f)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 */
bntobf(bf_t f,bn_t n)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 */
inttobf(bf_t r,long longval)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 */
bftoint(bf_t f)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) */
sign_bf(bf_t n)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| */
abs_bf(bf_t r,bf_t 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| */
abs_a_bf(bf_t 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. */
unsafe_inv_bf(bf_t r,bf_t n)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. */
unsafe_div_bf(bf_t r,bf_t n1,bf_t n2)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 - bftmp4 - global temp bigfloats */
497 /* SIDE-EFFECTS: */
498 /* n ends up as |n| */
unsafe_sqrt_bf(bf_t r,bf_t 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 copy_bf(bftmp4, r); /* r is used below, n is not */
563 unsafe_div_bf(bftmp3, n, bftmp4);
564 unsafe_add_a_bf(r, bftmp3);
565 half_a_bf(r);
566 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp3))) < 8 ) /* if match or almost match */
567 {
568 if (comp < 4 /* perfect or near perfect match */
569 || almost_match == 1) /* close enough for 2nd time */
570 break;
571 else /* this is the first time they almost matched */
572 almost_match++;
573 }
574 }
575
576 /* restore original values */
577 bflength = orig_bflength;
578 bnlength = orig_bnlength;
579 padding = orig_padding;
580 rlength = orig_rlength;
581 shiftfactor = orig_shiftfactor;
582 rbflength = orig_rbflength;
583 r = orig_r;
584 n = orig_n;
585
586 return r;
587 }
588
589 /********************************************************************/
590 /* exp(r) */
591 /* uses bftmp1, bftmp2, bftmp3 - global temp bigfloats */
exp_bf(bf_t r,bf_t n)592 bf_t exp_bf(bf_t r, bf_t n)
593 {
594 U16 fact=1;
595 S16 BIGDIST * testexp, BIGDIST * rexp;
596
597 testexp = (S16 BIGDIST *)(bftmp2+bflength);
598 rexp = (S16 BIGDIST *)(r+bflength);
599
600 if (is_bf_zero(n))
601 {
602 inttobf(r, 1);
603 return r;
604 }
605
606 /* use Taylor Series (very slow convergence) */
607 inttobf(r, 1); /* start with r=1.0 */
608 copy_bf(bftmp2, r);
609 for (;;)
610 {
611 copy_bf(bftmp1, n);
612 unsafe_mult_bf(bftmp3, bftmp2, bftmp1);
613 unsafe_div_bf_int(bftmp2, bftmp3, fact);
614 if (big_accessS16(testexp) < big_accessS16(rexp)-(bflength-2))
615 break; /* too small to register */
616 unsafe_add_a_bf(r, bftmp2);
617 fact++;
618 }
619
620 return r;
621 }
622
623 /********************************************************************/
624 /* ln(r) */
625 /* uses bftmp1 - bftmp6 - global temp bigfloats */
626 /* SIDE-EFFECTS: */
627 /* n ends up as |n| */
unsafe_ln_bf(bf_t r,bf_t n)628 bf_t unsafe_ln_bf(bf_t r, bf_t n)
629 {
630 int i, comp, almost_match=0;
631 LDBL f;
632 bf_t orig_r, orig_n, orig_bftmp5;
633 int orig_bflength,
634 orig_bnlength,
635 orig_padding,
636 orig_rlength,
637 orig_shiftfactor,
638 orig_rbflength;
639
640 /* use Newton's recursive method for zeroing in on ln(n): r=r+n*exp(-r)-1 */
641
642 if (is_bf_neg(n) || is_bf_zero(n))
643 { /* error, return largest neg value */
644 max_bf(r);
645 neg_a_bf(r);
646 return r;
647 }
648
649 f = bftofloat(n);
650 f = logl(f); /* approximate ln(x) */
651 /* no need to check overflow */
652 /* appears to be ok, do ln */
653
654 /* With Newton's Method, there is no need to calculate all the digits */
655 /* every time. The precision approximately doubles each iteration. */
656 /* Save original values. */
657 orig_bflength = bflength;
658 orig_bnlength = bnlength;
659 orig_padding = padding;
660 orig_rlength = rlength;
661 orig_shiftfactor = shiftfactor;
662 orig_rbflength = rbflength;
663 orig_r = r;
664 orig_n = n;
665 orig_bftmp5 = bftmp5;
666
667 /* calculate new starting values */
668 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
669 if (bnlength > orig_bnlength)
670 bnlength = orig_bnlength;
671 calc_lengths();
672
673 /* adjust pointers */
674 r = orig_r + orig_bflength - bflength;
675 n = orig_n + orig_bflength - bflength;
676 bftmp5 = orig_bftmp5 + orig_bflength - bflength;
677
678 floattobf(r, f); /* start with approximate ln */
679 neg_a_bf(r); /* -r */
680 copy_bf(bftmp5, r); /* -r */
681
682 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
683 {
684 /* adjust lengths */
685 bnlength <<= 1; /* double precision */
686 if (bnlength > orig_bnlength)
687 bnlength = orig_bnlength;
688 calc_lengths();
689 r = orig_r + orig_bflength - bflength;
690 n = orig_n + orig_bflength - bflength;
691 bftmp5 = orig_bftmp5 + orig_bflength - bflength;
692
693 exp_bf(bftmp6, r); /* exp(-r) */
694 unsafe_mult_bf(bftmp2, bftmp6, n); /* n*exp(-r) */
695 inttobf(bftmp4, 1);
696 unsafe_sub_a_bf(bftmp2, bftmp4); /* n*exp(-r) - 1 */
697 unsafe_sub_a_bf(r, bftmp2); /* -r - (n*exp(-r) - 1) */
698 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp5))) < 8 ) /* if match or almost match */
699 {
700 if (comp < 4 /* perfect or near perfect match */
701 || almost_match == 1) /* close enough for 2nd time */
702 break;
703 else /* this is the first time they almost matched */
704 almost_match++;
705 }
706 copy_bf(bftmp5, r); /* -r */
707 }
708
709 /* restore original values */
710 bflength = orig_bflength;
711 bnlength = orig_bnlength;
712 padding = orig_padding;
713 rlength = orig_rlength;
714 shiftfactor = orig_shiftfactor;
715 rbflength = orig_rbflength;
716 r = orig_r;
717 n = orig_n;
718 bftmp5 = orig_bftmp5;
719
720 neg_a_bf(r); /* -(-r) */
721 return r;
722 }
723
724 /********************************************************************/
725 /* sincos_bf(r) */
726 /* uses bftmp1 - bftmp2 - global temp bigfloats */
727 /* SIDE-EFFECTS: */
728 /* n ends up as |n| mod (pi/4) */
unsafe_sincos_bf(bf_t s,bf_t c,bf_t n)729 bf_t unsafe_sincos_bf(bf_t s, bf_t c, bf_t n)
730 {
731 U16 fact=2;
732 int k=0, i, halves;
733 int signcos=0, signsin=0, switch_sincos=0;
734 int sin_done=0, cos_done=0;
735 S16 BIGDIST * testexp, BIGDIST * cexp, BIGDIST * sexp;
736
737 testexp = (S16 BIGDIST *)(bftmp1+bflength);
738 cexp = (S16 BIGDIST *)(c+bflength);
739 sexp = (S16 BIGDIST *)(s+bflength);
740
741 #ifndef CALCULATING_BIG_PI
742 /* assure range 0 <= x < pi/4 */
743
744 if (is_bf_zero(n))
745 {
746 clear_bf(s); /* sin(0) = 0 */
747 inttobf(c, 1); /* cos(0) = 1 */
748 return s;
749 }
750
751 if (is_bf_neg(n))
752 {
753 signsin = !signsin; /* sin(-x) = -sin(x), odd; cos(-x) = cos(x), even */
754 neg_a_bf(n);
755 }
756 /* n >= 0 */
757
758 double_bf(bftmp1, bf_pi); /* 2*pi */
759 /* this could be done with remainders, but it would probably be slower */
760 while (cmp_bf(n, bftmp1) >= 0) /* while n >= 2*pi */
761 {
762 copy_bf(bftmp2, bftmp1);
763 unsafe_sub_a_bf(n, bftmp2);
764 }
765 /* 0 <= n < 2*pi */
766
767 copy_bf(bftmp1, bf_pi); /* pi */
768 if (cmp_bf(n, bftmp1) >= 0) /* if n >= pi */
769 {
770 unsafe_sub_a_bf(n, bftmp1);
771 signsin = !signsin;
772 signcos = !signcos;
773 }
774 /* 0 <= n < pi */
775
776 half_bf(bftmp1, bf_pi); /* pi/2 */
777 if (cmp_bf(n, bftmp1) > 0) /* if n > pi/2 */
778 {
779 copy_bf(bftmp2, bf_pi);
780 unsafe_sub_bf(n, bftmp2, n);
781 signcos = !signcos;
782 }
783 /* 0 <= n < pi/2 */
784
785 half_bf(bftmp1, bf_pi); /* pi/2 */
786 half_a_bf(bftmp1); /* pi/4 */
787 if (cmp_bf(n, bftmp1) > 0) /* if n > pi/4 */
788 {
789 copy_bf(bftmp2, n);
790 half_bf(bftmp1, bf_pi); /* pi/2 */
791 unsafe_sub_bf(n, bftmp1, bftmp2); /* pi/2 - n */
792 switch_sincos = !switch_sincos;
793 }
794 /* 0 <= n < pi/4 */
795
796 /* this looks redundant, but n could now be zero when it wasn't before */
797 if (is_bf_zero(n))
798 {
799 clear_bf(s); /* sin(0) = 0 */
800 inttobf(c, 1); /* cos(0) = 1 */
801 return s;
802 }
803
804
805 /* at this point, the double angle trig identities could be used as many */
806 /* times as desired to reduce the range to pi/8, pi/16, etc... Each time */
807 /* the range is cut in half, the number of iterations required is reduced */
808 /* by "quite a bit." It's just a matter of testing to see what gives the */
809 /* optimal results. */
810 /* halves = bflength / 10; */ /* this is experimental */
811 halves = 1;
812 for (i = 0; i < halves; i++)
813 half_a_bf(n);
814 #endif
815
816 /* use Taylor Series (very slow convergence) */
817 copy_bf(s, n); /* start with s=n */
818 inttobf(c, 1); /* start with c=1 */
819 copy_bf(bftmp1, n); /* the current x^n/n! */
820 do
821 {
822 /* even terms for cosine */
823 copy_bf(bftmp2, bftmp1);
824 unsafe_mult_bf(bftmp1, bftmp2, n);
825 div_a_bf_int(bftmp1, fact++);
826 if (!cos_done)
827 {
828 cos_done = (big_accessS16(testexp) < big_accessS16(cexp)-(bflength-2)); /* too small to register */
829 if (!cos_done)
830 {
831 if (k) /* alternate between adding and subtracting */
832 unsafe_add_a_bf(c, bftmp1);
833 else
834 unsafe_sub_a_bf(c, bftmp1);
835 }
836 }
837
838 /* odd terms for sine */
839 copy_bf(bftmp2, bftmp1);
840 unsafe_mult_bf(bftmp1, bftmp2, n);
841 div_a_bf_int(bftmp1, fact++);
842 if (!sin_done)
843 {
844 sin_done = (big_accessS16(testexp) < big_accessS16(sexp)-(bflength-2)); /* too small to register */
845 if (!sin_done)
846 {
847 if (k) /* alternate between adding and subtracting */
848 unsafe_add_a_bf(s, bftmp1);
849 else
850 unsafe_sub_a_bf(s, bftmp1);
851 }
852 }
853 k = !k; /* toggle */
854 #ifdef CALCULATING_BIG_PI
855 printf("."); /* lets you know it's doing something */
856 #endif
857 } while (!cos_done || !sin_done);
858
859 #ifndef CALCULATING_BIG_PI
860 /* now need to undo what was done by cutting angles in half */
861 for (i = 0; i < halves; i++)
862 {
863 unsafe_mult_bf(bftmp2, s, c); /* no need for safe mult */
864 double_bf(s, bftmp2); /* sin(2x) = 2*sin(x)*cos(x) */
865 unsafe_square_bf(bftmp2,c);
866 double_a_bf(bftmp2);
867 inttobf(bftmp1, 1);
868 unsafe_sub_bf(c, bftmp2, bftmp1); /* cos(2x) = 2*cos(x)*cos(x) - 1 */
869 }
870
871 if (switch_sincos)
872 {
873 copy_bf(bftmp1, s);
874 copy_bf(s, c);
875 copy_bf(c, bftmp1);
876 }
877 if (signsin)
878 neg_a_bf(s);
879 if (signcos)
880 neg_a_bf(c);
881 #endif
882
883 return s; /* return sine I guess */
884 }
885
886 /********************************************************************/
887 /* atan(r) */
888 /* uses bftmp1 - bftmp5 - global temp bigfloats */
889 /* SIDE-EFFECTS: */
890 /* n ends up as |n| or 1/|n| */
unsafe_atan_bf(bf_t r,bf_t n)891 bf_t unsafe_atan_bf(bf_t r, bf_t n)
892 {
893 int i, comp, almost_match=0, signflag=0;
894 LDBL f;
895 bf_t orig_r, orig_n, orig_bf_pi, orig_bftmp3;
896 int orig_bflength,
897 orig_bnlength,
898 orig_padding,
899 orig_rlength,
900 orig_shiftfactor,
901 orig_rbflength;
902 int large_arg;
903
904
905 /* use Newton's recursive method for zeroing in on atan(n): r=r-cos(r)(sin(r)-n*cos(r)) */
906
907 if (is_bf_neg(n))
908 {
909 signflag = 1;
910 neg_a_bf(n);
911 }
912
913 /* If n is very large, atanl() won't give enough decimal places to be a */
914 /* good enough initial guess for Newton's Method. If it is larger than */
915 /* say, 1, atan(n) = pi/2 - acot(n) = pi/2 - atan(1/n). */
916
917 f = bftofloat(n);
918 large_arg = f > 1.0;
919 if (large_arg)
920 {
921 unsafe_inv_bf(bftmp3, n);
922 copy_bf(n, bftmp3);
923 f = bftofloat(n);
924 }
925
926 clear_bf(bftmp3); /* not really necessary, but makes things more consistent */
927
928 /* With Newton's Method, there is no need to calculate all the digits */
929 /* every time. The precision approximately doubles each iteration. */
930 /* Save original values. */
931 orig_bflength = bflength;
932 orig_bnlength = bnlength;
933 orig_padding = padding;
934 orig_rlength = rlength;
935 orig_shiftfactor = shiftfactor;
936 orig_rbflength = rbflength;
937 orig_bf_pi = bf_pi;
938 orig_r = r;
939 orig_n = n;
940 orig_bftmp3 = bftmp3;
941
942 /* calculate new starting values */
943 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
944 if (bnlength > orig_bnlength)
945 bnlength = orig_bnlength;
946 calc_lengths();
947
948 /* adjust pointers */
949 r = orig_r + orig_bflength - bflength;
950 n = orig_n + orig_bflength - bflength;
951 bf_pi = orig_bf_pi + orig_bflength - bflength;
952 bftmp3 = orig_bftmp3 + orig_bflength - bflength;
953
954 f = atanl(f); /* approximate arctangent */
955 /* no need to check overflow */
956
957 floattobf(r, f); /* start with approximate atan */
958 copy_bf(bftmp3, r);
959
960 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
961 {
962 /* adjust lengths */
963 bnlength <<= 1; /* double precision */
964 if (bnlength > orig_bnlength)
965 bnlength = orig_bnlength;
966 calc_lengths();
967 r = orig_r + orig_bflength - bflength;
968 n = orig_n + orig_bflength - bflength;
969 bf_pi = orig_bf_pi + orig_bflength - bflength;
970 bftmp3 = orig_bftmp3 + orig_bflength - bflength;
971
972 #ifdef CALCULATING_BIG_PI
973 printf("\natan() loop #%i, bflength=%i\nsincos() loops\n", i, bflength);
974 #endif
975 unsafe_sincos_bf(bftmp4, bftmp5, bftmp3); /* sin(r), cos(r) */
976 copy_bf(bftmp3, r); /* restore bftmp3 from sincos_bf() */
977 copy_bf(bftmp1, bftmp5);
978 unsafe_mult_bf(bftmp2, n, bftmp1); /* n*cos(r) */
979 unsafe_sub_a_bf(bftmp4, bftmp2); /* sin(r) - n*cos(r) */
980 unsafe_mult_bf(bftmp1, bftmp5, bftmp4); /* cos(r) * (sin(r) - n*cos(r)) */
981 /* copy_bf(bftmp3, r); already done above! */
982 unsafe_sub_a_bf(r, bftmp1); /* r - cos(r) * (sin(r) - n*cos(r)) */
983 #ifdef CALCULATING_BIG_PI
984 putchar('\n');
985 bf_hexdump(r);
986 #endif
987 if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp3))) < 8 ) /* if match or almost match */
988 {
989 #ifdef CALCULATING_BIG_PI
990 printf("atan() loop comp=%i\n", comp);
991 #endif
992 if (comp < 4 /* perfect or near perfect match */
993 || almost_match == 1) /* close enough for 2nd time */
994 break;
995 else /* this is the first time they almost matched */
996 almost_match++;
997 }
998
999 #ifdef CALCULATING_BIG_PI
1000 if (bflength == orig_bflength && comp >= 8)
1001 printf("atan() loop comp=%i\n", comp);
1002 #endif
1003
1004 /* the following gets set above and not altered before we use it! */
1005 /* copy_bf(bftmp3, r); */ /* make a copy for later comparison */
1006 }
1007
1008 /* restore original values */
1009 bflength = orig_bflength;
1010 bnlength = orig_bnlength;
1011 padding = orig_padding;
1012 rlength = orig_rlength;
1013 shiftfactor = orig_shiftfactor;
1014 rbflength = orig_rbflength;
1015 bf_pi = orig_bf_pi;
1016 r = orig_r;
1017 n = orig_n;
1018 bftmp3 = orig_bftmp3;
1019
1020 if (large_arg)
1021 {
1022 half_bf(bftmp3, bf_pi); /* pi/2 */
1023 sub_a_bf(bftmp3, r); /* pi/2 - atan(1/n) */
1024 copy_bf(r, bftmp3);
1025 }
1026
1027 if (signflag)
1028 neg_a_bf(r);
1029 return r;
1030 }
1031
1032 /********************************************************************/
1033 /* atan2(r,ny,nx) */
1034 /* uses bftmp1 - bftmp6 - global temp bigfloats */
unsafe_atan2_bf(bf_t r,bf_t ny,bf_t nx)1035 bf_t unsafe_atan2_bf(bf_t r, bf_t ny, bf_t nx)
1036 {
1037 int signx, signy;
1038
1039 signx = sign_bf(nx);
1040 signy = sign_bf(ny);
1041
1042 if (signy == 0)
1043 {
1044 if (signx < 0)
1045 copy_bf(r, bf_pi); /* negative x axis, 180 deg */
1046 else /* signx >= 0 positive x axis, 0 */
1047 clear_bf(r);
1048 return(r);
1049 }
1050 if (signx == 0)
1051 {
1052 copy_bf(r, bf_pi); /* y axis */
1053 half_a_bf(r); /* +90 deg */
1054 if (signy < 0)
1055 neg_a_bf(r); /* -90 deg */
1056 return(r);
1057 }
1058
1059 if (signy < 0)
1060 neg_a_bf(ny);
1061 if (signx < 0)
1062 neg_a_bf(nx);
1063 unsafe_div_bf(bftmp6,ny,nx);
1064 unsafe_atan_bf(r, bftmp6);
1065 if (signx < 0)
1066 sub_bf(r,bf_pi,r);
1067 if(signy < 0)
1068 neg_a_bf(r);
1069 return(r);
1070 }
1071
1072 /**********************************************************************/
1073 /* The rest of the functions are "safe" versions of the routines that */
1074 /* have side effects which alter the parameters. */
1075 /* Most bf routines change values of parameters, not just the sign. */
1076 /**********************************************************************/
1077
1078 /**********************************************************************/
add_bf(bf_t r,bf_t n1,bf_t n2)1079 bf_t add_bf(bf_t r, bf_t n1, bf_t n2)
1080 {
1081 copy_bf(bftmpcpy1, n1);
1082 copy_bf(bftmpcpy2, n2);
1083 unsafe_add_bf(r, bftmpcpy1, bftmpcpy2);
1084 return r;
1085 }
1086
1087 /**********************************************************************/
add_a_bf(bf_t r,bf_t n)1088 bf_t add_a_bf(bf_t r, bf_t n)
1089 {
1090 copy_bf(bftmpcpy1, n);
1091 unsafe_add_a_bf(r, bftmpcpy1);
1092 return r;
1093 }
1094
1095 /**********************************************************************/
sub_bf(bf_t r,bf_t n1,bf_t n2)1096 bf_t sub_bf(bf_t r, bf_t n1, bf_t n2)
1097 {
1098 copy_bf(bftmpcpy1, n1);
1099 copy_bf(bftmpcpy2, n2);
1100 unsafe_sub_bf(r, bftmpcpy1, bftmpcpy2);
1101 return r;
1102 }
1103
1104 /**********************************************************************/
sub_a_bf(bf_t r,bf_t n)1105 bf_t sub_a_bf(bf_t r, bf_t n)
1106 {
1107 copy_bf(bftmpcpy1, n);
1108 unsafe_sub_a_bf(r, bftmpcpy1);
1109 return r;
1110 }
1111
1112 /**********************************************************************/
1113 /* mult and div only change sign */
full_mult_bf(bf_t r,bf_t n1,bf_t n2)1114 bf_t full_mult_bf(bf_t r, bf_t n1, bf_t n2)
1115 {
1116 copy_bf(bftmpcpy1, n1);
1117 copy_bf(bftmpcpy2, n2);
1118 unsafe_full_mult_bf(r, bftmpcpy1, bftmpcpy2);
1119 return r;
1120 }
1121
1122 /**********************************************************************/
mult_bf(bf_t r,bf_t n1,bf_t n2)1123 bf_t mult_bf(bf_t r, bf_t n1, bf_t n2)
1124 {
1125 copy_bf(bftmpcpy1, n1);
1126 copy_bf(bftmpcpy2, n2);
1127 unsafe_mult_bf(r, bftmpcpy1, bftmpcpy2);
1128 return r;
1129 }
1130
1131 /**********************************************************************/
full_square_bf(bf_t r,bf_t n)1132 bf_t full_square_bf(bf_t r, bf_t n)
1133 {
1134 copy_bf(bftmpcpy1, n);
1135 unsafe_full_square_bf(r, bftmpcpy1);
1136 return r;
1137 }
1138
1139 /**********************************************************************/
square_bf(bf_t r,bf_t n)1140 bf_t square_bf(bf_t r, bf_t n)
1141 {
1142 copy_bf(bftmpcpy1, n);
1143 unsafe_square_bf(r, bftmpcpy1);
1144 return r;
1145 }
1146
1147 /**********************************************************************/
mult_bf_int(bf_t r,bf_t n,U16 u)1148 bf_t mult_bf_int(bf_t r, bf_t n, U16 u)
1149 {
1150 copy_bf(bftmpcpy1, n);
1151 unsafe_mult_bf_int(r, bftmpcpy1, u);
1152 return r;
1153 }
1154
1155 /**********************************************************************/
div_bf_int(bf_t r,bf_t n,U16 u)1156 bf_t div_bf_int(bf_t r, bf_t n, U16 u)
1157 {
1158 copy_bf(bftmpcpy1, n);
1159 unsafe_div_bf_int(r, bftmpcpy1, u);
1160 return r;
1161 }
1162
1163 /**********************************************************************/
bftostr(char * s,int dec,bf_t r)1164 char *bftostr(char *s, int dec, bf_t r)
1165 {
1166 copy_bf(bftmpcpy1, r);
1167 unsafe_bftostr(s, dec, bftmpcpy1);
1168 return s;
1169 }
1170
1171 /**********************************************************************/
bftostr_e(char * s,int dec,bf_t r)1172 char *bftostr_e(char *s, int dec, bf_t r)
1173 {
1174 copy_bf(bftmpcpy1, r);
1175 unsafe_bftostr_e(s, dec, bftmpcpy1);
1176 return s;
1177 }
1178
1179 /**********************************************************************/
bftostr_f(char * s,int dec,bf_t r)1180 char *bftostr_f(char *s, int dec, bf_t r)
1181 {
1182 copy_bf(bftmpcpy1, r);
1183 unsafe_bftostr_f(s, dec, bftmpcpy1);
1184 return s;
1185 }
1186
1187 /**********************************************************************/
inv_bf(bf_t r,bf_t n)1188 bf_t inv_bf(bf_t r, bf_t n)
1189 {
1190 copy_bf(bftmpcpy1, n);
1191 unsafe_inv_bf(r, bftmpcpy1);
1192 return r;
1193 }
1194
1195 /**********************************************************************/
div_bf(bf_t r,bf_t n1,bf_t n2)1196 bf_t div_bf(bf_t r, bf_t n1, bf_t n2)
1197 {
1198 copy_bf(bftmpcpy1, n1);
1199 copy_bf(bftmpcpy2, n2);
1200 unsafe_div_bf(r, bftmpcpy1, bftmpcpy2);
1201 return r;
1202 }
1203
1204 /**********************************************************************/
sqrt_bf(bf_t r,bf_t n)1205 bf_t sqrt_bf(bf_t r, bf_t n)
1206 {
1207 copy_bf(bftmpcpy1, n);
1208 unsafe_sqrt_bf(r, bftmpcpy1);
1209 return r;
1210 }
1211
1212 /**********************************************************************/
ln_bf(bf_t r,bf_t n)1213 bf_t ln_bf(bf_t r, bf_t n)
1214 {
1215 copy_bf(bftmpcpy1, n);
1216 unsafe_ln_bf(r, bftmpcpy1);
1217 return r;
1218 }
1219
1220 /**********************************************************************/
sincos_bf(bf_t s,bf_t c,bf_t n)1221 bf_t sincos_bf(bf_t s, bf_t c, bf_t n)
1222 {
1223 copy_bf(bftmpcpy1, n);
1224 return unsafe_sincos_bf(s, c, bftmpcpy1);
1225 }
1226
1227 /**********************************************************************/
atan_bf(bf_t r,bf_t n)1228 bf_t atan_bf(bf_t r, bf_t n)
1229 {
1230 copy_bf(bftmpcpy1, n);
1231 unsafe_atan_bf(r, bftmpcpy1);
1232 return r;
1233 }
1234
1235 /**********************************************************************/
atan2_bf(bf_t r,bf_t ny,bf_t nx)1236 bf_t atan2_bf(bf_t r, bf_t ny, bf_t nx)
1237 {
1238 copy_bf(bftmpcpy1, ny);
1239 copy_bf(bftmpcpy2, nx);
1240 unsafe_atan2_bf(r, bftmpcpy1, bftmpcpy2);
1241 return r;
1242 }
1243
1244 /**********************************************************************/
is_bf_zero(bf_t n)1245 int is_bf_zero(bf_t n)
1246 {
1247 return !is_bf_not_zero(n);
1248 }
1249
1250 /************************************************************************/
1251 /* convert_bf -- convert bigfloat numbers from old to new lengths */
convert_bf(bf_t new,bf_t old,int newbflength,int oldbflength)1252 int convert_bf(bf_t new, bf_t old, int newbflength, int oldbflength)
1253 {
1254 int savebflength;
1255
1256 /* save lengths so not dependent on external environment */
1257 savebflength = bflength;
1258 bflength = newbflength;
1259 clear_bf(new);
1260 bflength = savebflength;
1261
1262 if(newbflength > oldbflength)
1263 _fmemcpy(new+newbflength-oldbflength,old,oldbflength+2);
1264 else
1265 _fmemcpy(new,old+oldbflength-newbflength,newbflength+2);
1266 return(0);
1267 }
1268
1269 /* The following used to be in bigfltc.c */
1270 /********************************************************************/
1271 /* normalize big float */
norm_bf(bf_t r)1272 bf_t norm_bf(bf_t r)
1273 {
1274 int scale;
1275 BYTE hi_byte;
1276 S16 BIGDIST *rexp;
1277
1278 rexp = (S16 BIGDIST *)(r+bflength);
1279
1280 /* check for overflow */
1281 hi_byte = r[bflength-1];
1282 if (hi_byte != 0x00 && hi_byte != 0xFF)
1283 {
1284 _fmemmove(r, r+1, bflength-1);
1285 r[bflength-1] = (BYTE)(hi_byte & 0x80 ? 0xFF : 0x00);
1286 big_setS16(rexp,big_accessS16(rexp)+(S16)1); /* exp */
1287 }
1288
1289 /* check for underflow */
1290 else
1291 {
1292 for (scale = 2; scale < bflength && r[bflength-scale] == hi_byte; scale++)
1293 ; /* do nothing */
1294 if (scale == bflength && hi_byte == 0) /* zero */
1295 big_setS16(rexp,0);
1296 else
1297 {
1298 scale -= 2;
1299 if (scale > 0) /* it did underflow */
1300 {
1301 _fmemmove(r+scale, r, bflength-scale-1);
1302 _fmemset(r, 0, scale);
1303 big_setS16(rexp,big_accessS16(rexp)-(S16)scale); /* exp */
1304 }
1305 }
1306 }
1307
1308 return r;
1309 }
1310
1311 /********************************************************************/
1312 /* normalize big float with forced sign */
1313 /* positive = 1, force to be positive */
1314 /* = 0, force to be negative */
norm_sign_bf(bf_t r,int positive)1315 void norm_sign_bf(bf_t r, int positive)
1316 {
1317 norm_bf(r);
1318 r[bflength-1] = (BYTE)(positive ? 0x00 : 0xFF);
1319 }
1320 /******************************************************/
1321 /* adjust n1, n2 for before addition or subtraction */
1322 /* by forcing exp's to match. */
1323 /* returns the value of the adjusted exponents */
adjust_bf_add(bf_t n1,bf_t n2)1324 S16 adjust_bf_add(bf_t n1, bf_t n2)
1325 {
1326 int scale, fill_byte;
1327 S16 rexp;
1328 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1329
1330 /* scale n1 or n2 */
1331 /* compare exp's */
1332 n1exp = (S16 BIGDIST *)(n1+bflength);
1333 n2exp = (S16 BIGDIST *)(n2+bflength);
1334 if (big_accessS16(n1exp) > big_accessS16(n2exp))
1335 { /* scale n2 */
1336 scale = big_accessS16(n1exp) - big_accessS16(n2exp); /* n1exp - n2exp */
1337 if (scale < bflength)
1338 {
1339 fill_byte = is_bf_neg(n2) ? 0xFF : 0x00;
1340 _fmemmove(n2, n2+scale, bflength-scale);
1341 _fmemset(n2+bflength-scale, fill_byte, scale);
1342 }
1343 else
1344 clear_bf(n2);
1345 big_setS16(n2exp,big_accessS16(n1exp)); /* *n2exp = *n1exp; set exp's = */
1346 rexp = big_accessS16(n2exp);
1347 }
1348 else if (big_accessS16(n1exp) < big_accessS16(n2exp))
1349 { /* scale n1 */
1350 scale = big_accessS16(n2exp) - big_accessS16(n1exp); /* n2exp - n1exp */
1351 if (scale < bflength)
1352 {
1353 fill_byte = is_bf_neg(n1) ? 0xFF : 0x00;
1354 _fmemmove(n1, n1+scale, bflength-scale);
1355 _fmemset(n1+bflength-scale, fill_byte, scale);
1356 }
1357 else
1358 clear_bf(n1);
1359 big_setS16(n1exp,big_accessS16(n2exp)); /* *n1exp = *n2exp; set exp's = */
1360 rexp = big_accessS16(n2exp);
1361 }
1362 else
1363 rexp = big_accessS16(n1exp);
1364 return rexp;
1365 }
1366
1367 /********************************************************************/
1368 /* r = max positive value */
max_bf(bf_t r)1369 bf_t max_bf(bf_t r)
1370 {
1371 inttobf(r, 1);
1372 big_set16(r+bflength, (S16)(LDBL_MAX_EXP/8));
1373 return r;
1374 }
1375
1376 /****************************************************************************/
1377 /* n1 != n2 ? */
1378 /* RETURNS: */
1379 /* if n1 == n2 returns 0 */
1380 /* if n1 > n2 returns a positive (bytes left to go when mismatch occurred) */
1381 /* if n1 < n2 returns a negative (bytes left to go when mismatch occurred) */
1382
cmp_bf(bf_t n1,bf_t n2)1383 int cmp_bf(bf_t n1, bf_t n2)
1384 {
1385 int i;
1386 int sign1, sign2;
1387 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1388 U16 value1, value2;
1389
1390
1391 /* compare signs */
1392 sign1 = sign_bf(n1);
1393 sign2 = sign_bf(n2);
1394 if (sign1 > sign2)
1395 return bflength;
1396 else if (sign1 < sign2)
1397 return -bflength;
1398 /* signs are the same */
1399
1400 /* compare exponents, using signed comparisons */
1401 n1exp = (S16 BIGDIST *)(n1+bflength);
1402 n2exp = (S16 BIGDIST *)(n2+bflength);
1403 if (big_accessS16(n1exp) > big_accessS16(n2exp))
1404 return sign1*(bflength);
1405 else if (big_accessS16(n1exp) < big_accessS16(n2exp))
1406 return -sign1*(bflength);
1407
1408 /* To get to this point, the signs must match */
1409 /* so unsigned comparison is ok. */
1410 /* two bytes at a time */
1411 for (i=bflength-2; i>=0; i-=2)
1412 {
1413 if ( (value1=big_access16(n1+i)) > (value2=big_access16(n2+i)) )
1414 { /* now determine which of the two bytes was different */
1415 if ( (value1&0xFF00) > (value2&0xFF00) ) /* compare just high bytes */
1416 return (i+2); /* high byte was different */
1417 else
1418 return (i+1); /* low byte was different */
1419 }
1420 else if (value1 < value2)
1421 { /* now determine which of the two bytes was different */
1422 if ( (value1&0xFF00) < (value2&0xFF00) ) /* compare just high bytes */
1423 return -(i+2); /* high byte was different */
1424 else
1425 return -(i+1); /* low byte was different */
1426 }
1427 }
1428 return 0;
1429 }
1430
1431 /********************************************************************/
1432 /* r < 0 ? */
1433 /* returns 1 if negative, 0 if positive or zero */
is_bf_neg(bf_t n)1434 int is_bf_neg(bf_t n)
1435 {
1436 return (S8)n[bflength-1] < 0;
1437 }
1438
1439 /********************************************************************/
1440 /* n != 0 ? */
1441 /* RETURNS: if n != 0 returns 1 */
1442 /* else returns 0 */
is_bf_not_zero(bf_t n)1443 int is_bf_not_zero(bf_t n)
1444 {
1445 int bnl;
1446 int retval;
1447
1448 bnl = bnlength;
1449 bnlength = bflength;
1450 retval = is_bn_not_zero(n);
1451 bnlength = bnl;
1452
1453 return retval;
1454 }
1455
1456 /********************************************************************/
1457 /* r = n1 + n2 */
1458 /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
unsafe_add_bf(bf_t r,bf_t n1,bf_t n2)1459 bf_t unsafe_add_bf(bf_t r, bf_t n1, bf_t n2)
1460 {
1461 int bnl;
1462 S16 BIGDIST *rexp;
1463
1464 if (is_bf_zero(n1))
1465 {
1466 copy_bf(r, n2);
1467 return r;
1468 }
1469 if (is_bf_zero(n2))
1470 {
1471 copy_bf(r, n1);
1472 return r;
1473 }
1474
1475 rexp = (S16 BIGDIST *)(r+bflength);
1476 big_setS16(rexp,adjust_bf_add(n1, n2));
1477
1478 bnl = bnlength;
1479 bnlength = bflength;
1480 add_bn(r, n1, n2);
1481 bnlength = bnl;
1482
1483 norm_bf(r);
1484 return r;
1485 }
1486
1487 /********************************************************************/
1488 /* r += n */
unsafe_add_a_bf(bf_t r,bf_t n)1489 bf_t unsafe_add_a_bf(bf_t r, bf_t n)
1490 {
1491 int bnl;
1492
1493 if (is_bf_zero(r))
1494 {
1495 copy_bf(r, n);
1496 return r;
1497 }
1498 if (is_bf_zero(n))
1499 {
1500 return r;
1501 }
1502
1503 adjust_bf_add(r, n);
1504
1505 bnl = bnlength;
1506 bnlength = bflength;
1507 add_a_bn(r, n);
1508 bnlength = bnl;
1509
1510 norm_bf(r);
1511
1512 return r;
1513 }
1514
1515 /********************************************************************/
1516 /* r = n1 - n2 */
1517 /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
unsafe_sub_bf(bf_t r,bf_t n1,bf_t n2)1518 bf_t unsafe_sub_bf(bf_t r, bf_t n1, bf_t n2)
1519 {
1520 int bnl;
1521 S16 BIGDIST *rexp;
1522
1523 if (is_bf_zero(n1))
1524 {
1525 neg_bf(r, n2);
1526 return r;
1527 }
1528 if (is_bf_zero(n2))
1529 {
1530 copy_bf(r, n1);
1531 return r;
1532 }
1533
1534 rexp = (S16 BIGDIST *)(r+bflength);
1535 big_setS16(rexp,adjust_bf_add(n1, n2));
1536
1537 bnl = bnlength;
1538 bnlength = bflength;
1539 sub_bn(r, n1, n2);
1540 bnlength = bnl;
1541
1542 norm_bf(r);
1543 return r;
1544 }
1545
1546 /********************************************************************/
1547 /* r -= n */
unsafe_sub_a_bf(bf_t r,bf_t n)1548 bf_t unsafe_sub_a_bf(bf_t r, bf_t n)
1549 {
1550 int bnl;
1551
1552 if (is_bf_zero(r))
1553 {
1554 neg_bf(r,n);
1555 return r;
1556 }
1557 if (is_bf_zero(n))
1558 {
1559 return r;
1560 }
1561
1562 adjust_bf_add(r, n);
1563
1564 bnl = bnlength;
1565 bnlength = bflength;
1566 sub_a_bn(r, n);
1567 bnlength = bnl;
1568
1569 norm_bf(r);
1570 return r;
1571 }
1572
1573 /********************************************************************/
1574 /* r = -n */
neg_bf(bf_t r,bf_t n)1575 bf_t neg_bf(bf_t r, bf_t n)
1576 {
1577 int bnl;
1578 S16 BIGDIST *rexp, BIGDIST *nexp;
1579
1580 rexp = (S16 BIGDIST *)(r+bflength);
1581 nexp = (S16 BIGDIST *)(n+bflength);
1582 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1583
1584 bnl = bnlength;
1585 bnlength = bflength;
1586 neg_bn(r, n);
1587 bnlength = bnl;
1588
1589 norm_bf(r);
1590 return r;
1591 }
1592
1593 /********************************************************************/
1594 /* r *= -1 */
neg_a_bf(bf_t r)1595 bf_t neg_a_bf(bf_t r)
1596 {
1597 int bnl;
1598
1599 bnl = bnlength;
1600 bnlength = bflength;
1601 neg_a_bn(r);
1602 bnlength = bnl;
1603
1604 norm_bf(r);
1605 return r;
1606 }
1607
1608 /********************************************************************/
1609 /* r = 2*n */
double_bf(bf_t r,bf_t n)1610 bf_t double_bf(bf_t r, bf_t n)
1611 {
1612 int bnl;
1613 S16 BIGDIST *rexp, BIGDIST *nexp;
1614
1615 rexp = (S16 BIGDIST *)(r+bflength);
1616 nexp = (S16 BIGDIST *)(n+bflength);
1617 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1618
1619 bnl = bnlength;
1620 bnlength = bflength;
1621 double_bn(r, n);
1622 bnlength = bnl;
1623
1624 norm_bf(r);
1625 return r;
1626 }
1627
1628 /********************************************************************/
1629 /* r *= 2 */
double_a_bf(bf_t r)1630 bf_t double_a_bf(bf_t r)
1631 {
1632 int bnl;
1633
1634 bnl = bnlength;
1635 bnlength = bflength;
1636 double_a_bn(r);
1637 bnlength = bnl;
1638
1639 norm_bf(r);
1640 return r;
1641 }
1642
1643 /********************************************************************/
1644 /* r = n/2 */
half_bf(bf_t r,bf_t n)1645 bf_t half_bf(bf_t r, bf_t n)
1646 {
1647 int bnl;
1648 S16 BIGDIST *rexp, BIGDIST *nexp;
1649
1650 rexp = (S16 BIGDIST *)(r+bflength);
1651 nexp = (S16 BIGDIST *)(n+bflength);
1652 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1653
1654 bnl = bnlength;
1655 bnlength = bflength;
1656 half_bn(r, n);
1657 bnlength = bnl;
1658
1659 norm_bf(r);
1660 return r;
1661 }
1662
1663 /********************************************************************/
1664 /* r /= 2 */
half_a_bf(bf_t r)1665 bf_t half_a_bf(bf_t r)
1666 {
1667 int bnl;
1668
1669 bnl = bnlength;
1670 bnlength = bflength;
1671 half_a_bn(r);
1672 bnlength = bnl;
1673
1674 norm_bf(r);
1675 return r;
1676 }
1677
1678 /************************************************************************/
1679 /* r = n1 * n2 */
1680 /* Note: r will be a double wide result, 2*bflength */
1681 /* n1 and n2 can be the same pointer */
1682 /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values */
unsafe_full_mult_bf(bf_t r,bf_t n1,bf_t n2)1683 bf_t unsafe_full_mult_bf(bf_t r, bf_t n1, bf_t n2)
1684 {
1685 int bnl, dbfl;
1686 S16 BIGDIST *rexp, BIGDIST *n1exp, BIGDIST *n2exp;
1687
1688 if (is_bf_zero(n1) || is_bf_zero(n2) )
1689 {
1690 bflength <<= 1;
1691 clear_bf(r);
1692 bflength >>= 1;
1693 return r;
1694 }
1695
1696 dbfl = 2*bflength; /* double width bflength */
1697 rexp = (S16 BIGDIST *)(r+dbfl); /* note: 2*bflength */
1698 n1exp = (S16 BIGDIST *)(n1+bflength);
1699 n2exp = (S16 BIGDIST *)(n2+bflength);
1700 /* add exp's */
1701 big_setS16(rexp, (S16)(big_accessS16(n1exp) + big_accessS16(n2exp)) );
1702
1703 bnl = bnlength;
1704 bnlength = bflength;
1705 unsafe_full_mult_bn(r, n1, n2);
1706 bnlength = bnl;
1707
1708 /* handle normalizing full mult on individual basis */
1709
1710 return r;
1711 }
1712
1713 /************************************************************************/
1714 /* r = n1 * n2 calculating only the top rlength bytes */
1715 /* Note: r will be of length rlength */
1716 /* 2*bflength <= rlength < bflength */
1717 /* n1 and n2 can be the same pointer */
1718 /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values */
unsafe_mult_bf(bf_t r,bf_t n1,bf_t n2)1719 bf_t unsafe_mult_bf(bf_t r, bf_t n1, bf_t n2)
1720 {
1721 int positive;
1722 int bnl, bfl, rl;
1723 int rexp;
1724 S16 BIGDIST *n1exp, BIGDIST *n2exp;
1725
1726 if (is_bf_zero(n1) || is_bf_zero(n2) )
1727 {
1728 clear_bf(r);
1729 return r;
1730 }
1731
1732 n1exp = (S16 BIGDIST *)(n1+bflength);
1733 n2exp = (S16 BIGDIST *)(n2+bflength);
1734 /* add exp's */
1735 rexp = big_accessS16(n1exp) + big_accessS16(n2exp);
1736
1737 positive = is_bf_neg(n1) == is_bf_neg(n2); /* are they the same sign? */
1738
1739 bnl = bnlength;
1740 bnlength = bflength;
1741 rl = rlength;
1742 rlength = rbflength;
1743 unsafe_mult_bn(r, n1, n2);
1744 bnlength = bnl;
1745 rlength = rl;
1746
1747 bfl = bflength;
1748 bflength = rbflength;
1749 big_set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
1750 norm_sign_bf(r, positive);
1751 bflength = bfl;
1752 _fmemmove(r, r+padding, bflength+2); /* shift back */
1753
1754 return r;
1755 }
1756
1757 /************************************************************************/
1758 /* r = n^2 */
1759 /* because of the symmetry involved, n^2 is much faster than n*n */
1760 /* for a bignumber of length l */
1761 /* n*n takes l^2 multiplications */
1762 /* n^2 takes (l^2+l)/2 multiplications */
1763 /* which is about 1/2 n*n as l gets large */
1764 /* uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
1765 /* */
1766 /* SIDE-EFFECTS: n is changed to its absolute value */
unsafe_full_square_bf(bf_t r,bf_t n)1767 bf_t unsafe_full_square_bf(bf_t r, bf_t n)
1768 {
1769 int bnl, dbfl;
1770 S16 BIGDIST *rexp, BIGDIST *nexp;
1771
1772 if (is_bf_zero(n))
1773 {
1774 bflength <<= 1;
1775 clear_bf(r);
1776 bflength >>= 1;
1777 return r;
1778 }
1779
1780 dbfl = 2*bflength; /* double width bflength */
1781 rexp = (S16 BIGDIST *)(r+dbfl); /* note: 2*bflength */
1782 nexp = (S16 BIGDIST *)(n+bflength);
1783 big_setS16(rexp, 2 * big_accessS16(nexp));
1784
1785 bnl = bnlength;
1786 bnlength = bflength;
1787 unsafe_full_square_bn(r, n);
1788 bnlength = bnl;
1789
1790 /* handle normalizing full mult on individual basis */
1791
1792 return r;
1793 }
1794
1795
1796 /************************************************************************/
1797 /* r = n^2 */
1798 /* because of the symmetry involved, n^2 is much faster than n*n */
1799 /* for a bignumber of length l */
1800 /* n*n takes l^2 multiplications */
1801 /* n^2 takes (l^2+l)/2 multiplications */
1802 /* which is about 1/2 n*n as l gets large */
1803 /* uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
1804 /* */
1805 /* Note: r will be of length rlength */
1806 /* 2*bflength >= rlength > bflength */
1807 /* SIDE-EFFECTS: n is changed to its absolute value */
unsafe_square_bf(bf_t r,bf_t n)1808 bf_t unsafe_square_bf(bf_t r, bf_t n)
1809 {
1810 int bnl, bfl, rl;
1811 int rexp;
1812 S16 BIGDIST *nexp;
1813
1814 if (is_bf_zero(n))
1815 {
1816 clear_bf(r);
1817 return r;
1818 }
1819
1820 nexp = (S16 BIGDIST *)(n+bflength);
1821 rexp = (S16)(2 * big_accessS16(nexp));
1822
1823 bnl = bnlength;
1824 bnlength = bflength;
1825 rl = rlength;
1826 rlength = rbflength;
1827 unsafe_square_bn(r, n);
1828 bnlength = bnl;
1829 rlength = rl;
1830
1831 bfl = bflength;
1832 bflength = rbflength;
1833 big_set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
1834
1835 norm_sign_bf(r, 1);
1836 bflength = bfl;
1837 _fmemmove(r, r+padding, bflength+2); /* shift back */
1838
1839 return r;
1840 }
1841
1842 /********************************************************************/
1843 /* r = n * u where u is an unsigned integer */
1844 /* SIDE-EFFECTS: n can be "de-normalized" and lose precision */
unsafe_mult_bf_int(bf_t r,bf_t n,U16 u)1845 bf_t unsafe_mult_bf_int(bf_t r, bf_t n, U16 u)
1846 {
1847 int positive;
1848 int bnl;
1849 S16 BIGDIST *rexp, BIGDIST *nexp;
1850
1851 rexp = (S16 BIGDIST *)(r+bflength);
1852 nexp = (S16 BIGDIST *)(n+bflength);
1853 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1854
1855 positive = !is_bf_neg(n);
1856
1857 /*
1858 if u > 0x00FF, then the integer part of the mantissa will overflow the
1859 2 byte (16 bit) integer size. Therefore, make adjustment before
1860 multiplication is performed.
1861 */
1862 if (u > 0x00FF)
1863 { /* un-normalize n */
1864 _fmemmove(n, n+1, bflength-1); /* this sign extends as well */
1865 big_setS16(rexp,big_accessS16(rexp)+(S16)1);
1866 }
1867
1868 bnl = bnlength;
1869 bnlength = bflength;
1870 mult_bn_int(r, n, u);
1871 bnlength = bnl;
1872
1873 norm_sign_bf(r, positive);
1874 return r;
1875 }
1876
1877 /********************************************************************/
1878 /* r *= u where u is an unsigned integer */
mult_a_bf_int(bf_t r,U16 u)1879 bf_t mult_a_bf_int(bf_t r, U16 u)
1880 {
1881 int positive;
1882 int bnl;
1883 S16 BIGDIST *rexp;
1884
1885 rexp = (S16 BIGDIST *)(r+bflength);
1886 positive = !is_bf_neg(r);
1887
1888 /*
1889 if u > 0x00FF, then the integer part of the mantissa will overflow the
1890 2 byte (16 bit) integer size. Therefore, make adjustment before
1891 multiplication is performed.
1892 */
1893 if (u > 0x00FF)
1894 { /* un-normalize n */
1895 _fmemmove(r, r+1, bflength-1); /* this sign extends as well */
1896 big_setS16(rexp,big_accessS16(rexp)+(S16)1);
1897 }
1898
1899 bnl = bnlength;
1900 bnlength = bflength;
1901 mult_a_bn_int(r, u);
1902 bnlength = bnl;
1903
1904 norm_sign_bf(r, positive);
1905 return r;
1906 }
1907
1908 /********************************************************************/
1909 /* r = n / u where u is an unsigned integer */
unsafe_div_bf_int(bf_t r,bf_t n,U16 u)1910 bf_t unsafe_div_bf_int(bf_t r, bf_t n, U16 u)
1911 {
1912 int bnl;
1913 S16 BIGDIST *rexp, BIGDIST *nexp;
1914
1915 if (u == 0) /* division by zero */
1916 {
1917 max_bf(r);
1918 if (is_bf_neg(n))
1919 neg_a_bf(r);
1920 return r;
1921 }
1922
1923 rexp = (S16 BIGDIST *)(r+bflength);
1924 nexp = (S16 BIGDIST *)(n+bflength);
1925 big_setS16(rexp, big_accessS16(nexp)); /* *rexp = *nexp; */
1926
1927 bnl = bnlength;
1928 bnlength = bflength;
1929 unsafe_div_bn_int(r, n, u);
1930 bnlength = bnl;
1931
1932 norm_bf(r);
1933 return r;
1934 }
1935
1936 /********************************************************************/
1937 /* r /= u where u is an unsigned integer */
div_a_bf_int(bf_t r,U16 u)1938 bf_t div_a_bf_int(bf_t r, U16 u)
1939 {
1940 int bnl;
1941
1942 if (u == 0) /* division by zero */
1943 {
1944 if (is_bf_neg(r))
1945 {
1946 max_bf(r);
1947 neg_a_bf(r);
1948 }
1949 else
1950 {
1951 max_bf(r);
1952 }
1953 return r;
1954 }
1955
1956 bnl = bnlength;
1957 bnlength = bflength;
1958 div_a_bn_int(r, u);
1959 bnlength = bnl;
1960
1961 norm_bf(r);
1962 return r;
1963 }
1964
1965 /********************************************************************/
1966 /* extracts the mantissa and exponent of f */
1967 /* finds m and n such that 1<=|m|<b and f = m*b^n */
1968 /* n is stored in *exp_ptr and m is returned, sort of like frexp() */
extract_value(LDBL f,LDBL b,int * exp_ptr)1969 LDBL extract_value(LDBL f, LDBL b, int *exp_ptr)
1970 {
1971 int n;
1972 LDBL af, ff, orig_b;
1973 LDBL value[15];
1974 unsigned powertwo;
1975
1976 if (b <= 0 || f == 0)
1977 {
1978 *exp_ptr = 0;
1979 return 0;
1980 }
1981
1982 orig_b = b;
1983 af = f >= 0 ? f: -f; /* abs value */
1984 ff = af > 1 ? af : 1/af;
1985 n = 0;
1986 powertwo = 1;
1987 while (b < ff)
1988 {
1989 value[n] = b;
1990 n++;
1991 powertwo <<= 1;
1992 b *= b;
1993 }
1994
1995 *exp_ptr = 0;
1996 for (; n > 0; n--)
1997 {
1998 powertwo >>= 1;
1999 if (value[n-1] < ff)
2000 {
2001 ff /= value[n-1];
2002 *exp_ptr += powertwo;
2003 }
2004 }
2005 if (f < 0)
2006 ff = -ff;
2007 if (af < 1)
2008 {
2009 ff = orig_b/ff;
2010 *exp_ptr = -*exp_ptr - 1;
2011 }
2012
2013 return ff;
2014 }
2015
2016 /********************************************************************/
2017 /* calculates and returns the value of f*b^n */
2018 /* sort of like ldexp() */
scale_value(LDBL f,LDBL b,int n)2019 LDBL scale_value( LDBL f, LDBL b , int n )
2020 {
2021 LDBL total=1;
2022 int an;
2023
2024 if (b == 0 || f == 0)
2025 return 0;
2026
2027 if (n == 0)
2028 return f;
2029
2030 an = abs(n);
2031
2032 while (an != 0)
2033 {
2034 if (an & 0x0001)
2035 total *= b;
2036 b *= b;
2037 an >>= 1;
2038 }
2039
2040 if (n > 0)
2041 f *= total;
2042 else /* n < 0 */
2043 f /= total;
2044 return f;
2045 }
2046
2047 /********************************************************************/
2048 /* extracts the mantissa and exponent of f */
2049 /* finds m and n such that 1<=|m|<10 and f = m*10^n */
2050 /* n is stored in *exp_ptr and m is returned, sort of like frexp() */
extract_10(LDBL f,int * exp_ptr)2051 LDBL extract_10(LDBL f, int *exp_ptr)
2052 {
2053 return extract_value(f, 10, exp_ptr);
2054 }
2055
2056 /********************************************************************/
2057 /* calculates and returns the value of f*10^n */
2058 /* sort of like ldexp() */
scale_10(LDBL f,int n)2059 LDBL scale_10( LDBL f, int n )
2060 {
2061 return scale_value( f, 10, n );
2062 }
2063
2064
2065
2066 /* big10flt.c - C routines for base 10 big floating point numbers */
2067
2068 /**********************************************************
2069 (Just when you thought it was safe to go back in the water.)
2070 Just when you thought you seen every type of format possible,
2071 16 bit integer, 32 bit integer, double, long double, mpmath,
2072 bn_t, bf_t, I now give you bf10_t (big float base 10)!
2073
2074 Why, because this is the only way (I can think of) to properly do a
2075 bftostr() without rounding errors. With out this, then
2076 -1.9999999999( > LDBL_DIG of 9's)9999999123456789...
2077 will round to -2.0. The good news is that we only need to do two
2078 mathematical operations: multiplication and division by integers
2079
2080 bf10_t format: (notice the position of the MSB and LSB)
2081
2082 MSB LSB
2083 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2084 n <><------------- dec --------------><> <->
2085 1 byte pad 1 byte rounding 2 byte exponent.
2086
2087 total length = dec + 4
2088
2089 ***********************************************************/
2090
2091 /**********************************************************************/
2092 /* unsafe_bftobf10() - converts a bigfloat into a bigfloat10 */
2093 /* n - pointer to a bigfloat */
2094 /* r - result array of BYTE big enough to hold the bf10_t number */
2095 /* dec - number of decimals, not including the one extra for rounding */
2096 /* SIDE-EFFECTS: n is changed to |n|. Make copy of n if necessary. */
2097
unsafe_bftobf10(bf10_t r,int dec,bf_t n)2098 bf10_t unsafe_bftobf10(bf10_t r, int dec, bf_t n)
2099 {
2100 int d;
2101 int power256;
2102 int p;
2103 int bnl;
2104 bf_t onesbyte;
2105 bf10_t power10;
2106
2107 if (is_bf_zero(n))
2108 { /* in scientific notation, the leading digit can't be zero */
2109 r[1] = (BYTE)0; /* unless the number is zero */
2110 return r;
2111 }
2112
2113 onesbyte = n + bflength - 1; /* really it's n+bflength-2 */
2114 power256 = (S16)big_access16(n + bflength) + 1; /* so adjust power256 by 1 */
2115
2116 if (dec == 0)
2117 dec = decimals;
2118 dec++; /* one extra byte for rounding */
2119 power10 = r + dec + 1;
2120
2121 if (is_bf_neg(n))
2122 {
2123 neg_a_bf(n);
2124 r[0] = 1; /* sign flag */
2125 }
2126 else
2127 r[0] = 0;
2128
2129 p = -1; /* multiply by 10 right away */
2130 bnl = bnlength;
2131 bnlength = bflength;
2132 for (d=1; d<=dec; d++)
2133 {
2134 /* pretend it's a bn_t instead of a bf_t */
2135 /* this leaves n un-normalized, which is what we want here */
2136 mult_a_bn_int(n, 10);
2137
2138 r[d] = *onesbyte;
2139 if (d == 1 && r[d] == 0)
2140 {
2141 d = 0; /* back up a digit */
2142 p--; /* and decrease by a factor of 10 */
2143 }
2144 *onesbyte = 0;
2145 }
2146 bnlength = bnl;
2147 big_set16(power10, (U16)p); /* save power of ten */
2148
2149 /* the digits are all read in, now scale it by 256^power256 */
2150 if (power256 > 0)
2151 for (d=0; d<power256; d++)
2152 mult_a_bf10_int(r, dec, 256);
2153
2154 else if (power256 < 0)
2155 for (d=0; d>power256; d--)
2156 div_a_bf10_int(r, dec, 256);
2157
2158 /* else power256 is zero, don't do anything */
2159
2160 /* round the last digit */
2161 if (r[dec] >= 5)
2162 {
2163 d = dec-1;
2164 while (d > 0) /* stop before you get to the sign flag */
2165 {
2166 r[d]++; /* round up */
2167 if (r[d] < 10)
2168 {
2169 d = -1; /* flag for below */
2170 break; /* finished rounding */
2171 }
2172 r[d] = 0;
2173 d--;
2174 }
2175 if (d == 0) /* rounding went back to the first digit and it overflowed */
2176 {
2177 r[1] = 0;
2178 _fmemmove(r+2, r+1, dec-1);
2179 r[1] = 1;
2180 p = (S16)big_access16(power10);
2181 big_set16(power10, (U16)(p+1));
2182 }
2183 }
2184 r[dec] = 0; /* truncate the rounded digit */
2185
2186 return r;
2187 }
2188
2189
2190 /**********************************************************************/
2191 /* mult_a_bf10_int() */
2192 /* r *= n */
2193 /* dec - number of decimals, including the one extra for rounding */
2194
mult_a_bf10_int(bf10_t r,int dec,U16 n)2195 bf10_t mult_a_bf10_int(bf10_t r, int dec, U16 n)
2196 {
2197 int signflag;
2198 int d, p;
2199 unsigned value, overflow;
2200 bf10_t power10;
2201
2202 if (r[1] == 0 || n == 0)
2203 {
2204 r[1] = 0;
2205 return r;
2206 }
2207
2208 power10 = r + dec + 1;
2209 p = (S16)big_access16(power10);
2210
2211 signflag = r[0]; /* r[0] to be used as a padding */
2212 overflow = 0;
2213 for (d = dec; d>0; d--)
2214 {
2215 value = r[d] * n + overflow;
2216 r[d] = (BYTE)(value % 10);
2217 overflow = value / 10;
2218 }
2219 while (overflow)
2220 {
2221 p++;
2222 _fmemmove(r+2, r+1, dec-1);
2223 r[1] = (BYTE)(overflow % 10);
2224 overflow = overflow / 10;
2225 }
2226 big_set16(power10, (U16)p); /* save power of ten */
2227 r[0] = (BYTE)signflag; /* restore sign flag */
2228 return r;
2229 }
2230
2231 /**********************************************************************/
2232 /* div_a_bf10_int() */
2233 /* r /= n */
2234 /* dec - number of decimals, including the one extra for rounding */
2235
div_a_bf10_int(bf10_t r,int dec,U16 n)2236 bf10_t div_a_bf10_int (bf10_t r, int dec, U16 n)
2237 {
2238 int src, dest, p;
2239 unsigned value, remainder;
2240 bf10_t power10;
2241
2242 if (r[1] == 0 || n == 0)
2243 {
2244 r[1] = 0;
2245 return r;
2246 }
2247
2248 power10 = r + dec + 1;
2249 p = (S16)big_access16(power10);
2250
2251 remainder = 0;
2252 for (src=dest=1; src<=dec; dest++, src++)
2253 {
2254 value = 10*remainder + r[src];
2255 r[dest] = (BYTE)(value / n);
2256 remainder = value % n;
2257 if (dest == 1 && r[dest] == 0)
2258 {
2259 dest = 0; /* back up a digit */
2260 p--; /* and decrease by a factor of 10 */
2261 }
2262 }
2263 for (; dest<=dec; dest++)
2264 {
2265 value = 10*remainder;
2266 r[dest] = (BYTE)(value / n);
2267 remainder = value % n;
2268 if (dest == 1 && r[dest] == 0)
2269 {
2270 dest = 0; /* back up a digit */
2271 p--; /* and decrease by a factor of 10 */
2272 }
2273 }
2274
2275 big_set16(power10, (U16)p); /* save power of ten */
2276 return r;
2277 }
2278
2279
2280 /*************************************************************************/
2281 /* bf10tostr_e() */
2282 /* Takes a bf10 number and converts it to an ascii string, sci. notation */
2283 /* dec - number of decimals, not including the one extra for rounding */
2284
bf10tostr_e(char * s,int dec,bf10_t n)2285 char *bf10tostr_e(char *s, int dec, bf10_t n)
2286 {
2287 int d, p;
2288 bf10_t power10;
2289
2290 if (n[1] == 0)
2291 {
2292 strcpy(s, "0.0");
2293 return s;
2294 }
2295
2296 if (dec == 0)
2297 dec = decimals;
2298 dec++; /* one extra byte for rounding */
2299 power10 = n + dec + 1;
2300 p = (S16)big_access16(power10);
2301
2302 /* if p is negative, it is not necessary to show all the decimal places */
2303 if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
2304 {
2305 dec = dec + p;
2306 if (dec < 8) /* let's keep at least a few */
2307 dec = 8;
2308 }
2309
2310 if (n[0] == 1) /* sign flag */
2311 *(s++) = '-';
2312 *(s++) = (char)(n[1] + '0');
2313 *(s++) = '.';
2314 for (d=2; d<=dec; d++)
2315 {
2316 *(s++) = (char)(n[d] + '0');
2317 }
2318 /* clean up trailing 0's */
2319 while (*(s-1) == '0')
2320 s--;
2321 if (*(s-1) == '.') /* put at least one 0 after the decimal */
2322 *(s++) = '0';
2323 sprintf(s, "e%d", p);
2324 return s;
2325 }
2326
2327 /****************************************************************************/
2328 /* bf10tostr_f() */
2329 /* Takes a bf10 number and converts it to an ascii string, decimal notation */
2330
bf10tostr_f(char * s,int dec,bf10_t n)2331 char *bf10tostr_f(char *s, int dec, bf10_t n)
2332 {
2333 int d, p;
2334 bf10_t power10;
2335
2336 if (n[1] == 0)
2337 {
2338 strcpy(s, "0.0");
2339 return s;
2340 }
2341
2342 if (dec == 0)
2343 dec = decimals;
2344 dec++; /* one extra byte for rounding */
2345 power10 = n + dec + 1;
2346 p = (S16)big_access16(power10);
2347
2348 /* if p is negative, it is not necessary to show all the decimal places */
2349 if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
2350 {
2351 dec = dec + p;
2352 if (dec < 8) /* let's keep at least a few */
2353 dec = 8;
2354 }
2355
2356 if (n[0] == 1) /* sign flag */
2357 *(s++) = '-';
2358 if (p >= 0)
2359 {
2360 for (d=1; d<=p+1; d++)
2361 *(s++) = (char)(n[d] + '0');
2362 *(s++) = '.';
2363 for (; d<=dec; d++)
2364 *(s++) = (char)(n[d] + '0');
2365 }
2366 else
2367 {
2368 *(s++) = '0';
2369 *(s++) = '.';
2370 for (d=0; d>p+1; d--)
2371 *(s++) = '0';
2372 for (d=1; d<=dec; d++)
2373 *(s++) = (char)(n[d] + '0');
2374 }
2375
2376 /* clean up trailing 0's */
2377 while (*(s-1) == '0')
2378 s--;
2379 if (*(s-1) == '.') /* put at least one 0 after the decimal */
2380 *(s++) = '0';
2381 *s = '\0'; /* terminating nul */
2382 return s;
2383 }
2384