1 /* bignum.c - C routines for bignumbers */
2
3 /*
4 Wesley Loewer's Big Numbers. (C) 1994-95, Wesley B. Loewer
5
6 The bignumber format is simply a signed integer of variable length. The
7 bytes are stored in reverse order (least significant byte first, most
8 significant byte last). The sign bit is the highest bit of the most
9 significant byte. Negatives are stored in 2's complement form. The byte
10 length of the bignumbers must be a multiple of 4 for 386+ operations, and
11 a multiple of 2 for 8086/286 and non 80x86 machines.
12
13 Some of the arithmetic operations coded here may alter some of the
14 operands used. Therefore, take note of the SIDE-EFFECTS listed with each
15 procedure. If the value of an operand needs to be retained, just use
16 copy_bn() first. This was done for speed sake to avoid unnecessary
17 copying. If space is at such a premium that copying it would be
18 difficult, some of the operations only change the sign of the value. In
19 this case, the original could be optained by calling neg_a_bn().
20
21 Most of the bignumber routines operate on true integers. Some of the
22 procedures, however, are designed specifically for fixed decimal point
23 operations. The results and how the results are interpreted depend on
24 where the implied decimal point is located. The routines that depend on
25 where the decimal is located are: strtobn(), bntostr(), bntoint(), inttobn(),
26 bntofloat(), floattobn(), inv_bn(), div_bn(). The number of bytes
27 designated for the integer part may be 1, 2, or 4.
28
29
30 BIGNUMBER FORMAT:
31 The following is a discription of the bignumber format and associated
32 variables. The number is stored in reverse order (Least Significant Byte,
33 LSB, stored first in memory, Most Significant Byte, MSB, stored last).
34 Each '_' below represents a block of memory used for arithmetic (1 block =
35 4 bytes on 386+, 2 bytes on 286-). All lengths are given in bytes.
36
37 LSB MSB
38 _ _ _ _ _ _ _ _ _ _ _ _
39 n <----------- bnlength ----------->
40 intlength ---> <---
41
42 bnlength = the length in bytes of the bignumber
43 intlength = the number of bytes used to represent the integer part of
44 the bignumber. Possible values are 1, 2, or 4. This
45 determines the largest number that can be represented by
46 the bignumber.
47 intlength = 1, max value = 127.99...
48 intlength = 2, max value = 32,767.99...
49 intlength = 4, max value = 2,147,483,647.99...
50
51
52 FULL DOUBLE PRECISION MULTIPLICATION:
53 ( full_mult_bn(), full_square_bn() )
54
55 The product of two bignumbers, n1 and n2, will be a result, r, which is
56 a double wide bignumber. The integer part will also be twice as wide,
57 thereby eliminating the possiblity of overflowing the number.
58
59 LSB MSB
60 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
61 r <--------------------------- 2*bnlength ----------------------------->
62 2*intlength ---> <---
63
64 If this double wide bignumber, r, needs to be converted to a normal,
65 single width bignumber, this is easily done with pointer arithmetic. The
66 converted value starts at r+shiftfactor (where shiftfactor =
67 bnlength-intlength) and continues for bnlength bytes. The lower order
68 bytes and the upper integer part of the double wide number can then be
69 ignored.
70
71 LSB MSB
72 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
73 r <--------------------------- 2*bnlength ----------------------------->
74 2*intlength ---> <---
75 LSB MSB
76 r+shiftfactor <---------- bnlength ------------>
77 intlength ---> <---
78
79
80 PARTIAL PRECISION MULTIPLICATION:
81 ( mult_bn(), square_bn() )
82
83 In most cases, full double precision multiplication is not necessary. The
84 lower order bytes are usually thrown away anyway. The non-"full"
85 multiplication routines only calculate rlength bytes in the result. The
86 value of rlength must be in the range: 2*bnlength <= rlength < bnlength.
87 The amount by which rlength exceeds bnlength accounts for the extra bytes
88 that must be multiplied so that the first bnlength bytes are correct.
89 These extra bytes are refered to in the code as the "padding," that is:
90 rlength=bnlength+padding.
91
92 All three of the values, bnlength, rlength, and therefore padding, must be
93 multiples of the size of memory blocks being used for arithmetic (2 on
94 8086/286 and 4 on 386+). Typically, the padding is 2*blocksize. In the
95 case where bnlength=blocksize, padding can only be blocksize to keep
96 rlength from being too big.
97
98 The product of two bignumbers, n1 and n2, will then be a result, r, which
99 is of length rlength. The integer part will be twice as wide, thereby
100 eliminating the possiblity of overflowing the number.
101
102 LSB MSB
103 _ _ _ _ _ _ _ _ _ _ _ _ _ _
104 r <---- rlength = bnlength+padding ------>
105 2*intlength ---> <---
106
107 If r needs to be converted to a normal, single width bignumber, this is
108 easily done with pointer arithmetic. The converted value starts at
109 r+shiftfactor (where shiftfactor = padding-intlength) and continues for
110 bnlength bytes. The lower order bytes and the upper integer part of the
111 double wide number can then be ignored.
112
113 LSB MSB
114 _ _ _ _ _ _ _ _ _ _ _ _ _ _
115 r <---- rlength = bnlength+padding ------>
116 2*intlength ---> <---
117 LSB MSB
118 r+shiftfactor <---------- bnlength --------->
119 intlength ---> <---
120 */
121
122 /************************************************************************/
123 /* There are three parts to the bignum library: */
124 /* */
125 /* 1) bignum.c - initialization, general routines, routines that would */
126 /* not be speeded up much with assembler. */
127 /* */
128 /* 2) bignuma.asm - hand coded assembler routines. */
129 /* */
130 /* 3) bignumc.c - portable C versions of routines in bignuma.asm */
131 /* */
132 /************************************************************************/
133
134 #include <string.h>
135 /* see Fractint.c for a description of the "include" hierarchy */
136 #include "port.h"
137 #include "prototyp.h"
138 #include "big.h"
139 #ifndef BIG_ANSI_C
140 #include <malloc.h>
141 #endif
142
143 /*************************************************************************
144 * The original bignumber code was written specifically for a Little Endian
145 * system (80x86). The following is not particularly efficient, but was
146 * simple to incorporate. If speed with a Big Endian machine is critical,
147 * the bignumber format could be reversed.
148 **************************************************************************/
149 #ifdef ACCESS_BY_BYTE
big_access32(BYTE BIGDIST * addr)150 U32 big_access32(BYTE BIGDIST *addr)
151 {
152 return addr[0] | ((U32)addr[1]<<8) | ((U32)addr[2]<<16) | ((U32)addr[3]<<24);
153 }
154
big_access16(BYTE BIGDIST * addr)155 U16 big_access16(BYTE BIGDIST *addr)
156 {
157 return (U16)addr[0] | ((U16)addr[1]<<8);
158 }
159
big_accessS16(S16 BIGDIST * addr)160 S16 big_accessS16(S16 BIGDIST *addr)
161 {
162 return (S16)((BYTE *)addr)[0] | ((S16)((BYTE *)addr)[1]<<8);
163 }
164
big_set32(BYTE BIGDIST * addr,U32 val)165 U32 big_set32(BYTE BIGDIST *addr, U32 val)
166 {
167 addr[0] = (BYTE)(val&0xff);
168 addr[1] = (BYTE)((val>>8)&0xff);
169 addr[2] = (BYTE)((val>>16)&0xff);
170 addr[3] = (BYTE)((val>>24)&0xff);
171 return val;
172 }
173
big_set16(BYTE BIGDIST * addr,U16 val)174 U16 big_set16(BYTE BIGDIST *addr, U16 val)
175 {
176 addr[0] = (BYTE)(val&0xff);
177 addr[1] = (BYTE)((val>>8)&0xff);
178 return val;
179 }
180
big_setS16(S16 BIGDIST * addr,S16 val)181 S16 big_setS16(S16 BIGDIST *addr, S16 val)
182 {
183 ((BYTE *)addr)[0] = (BYTE)(val&0xff);
184 ((BYTE *)addr)[1] = (BYTE)((val>>8)&0xff);
185 return val;
186 }
187
188 #endif
189
190 #if (_MSC_VER >= 700)
191 #pragma code_seg ("bignum1_text") /* place following in an overlay */
192 #endif
193
194 /************************************************************************/
195 /* convert_bn -- convert bignum numbers from old to new lengths */
convert_bn(bn_t new,bn_t old,int newbnlength,int newintlength,int oldbnlength,int oldintlength)196 int convert_bn(bn_t new, bn_t old, int newbnlength, int newintlength,
197 int oldbnlength, int oldintlength)
198 {
199 int savebnlength, saveintlength;
200
201 /* save lengths so not dependent on external environment */
202 saveintlength = intlength;
203 savebnlength = bnlength;
204
205 intlength = newintlength;
206 bnlength = newbnlength;
207 clear_bn(new);
208
209 if(newbnlength - newintlength > oldbnlength - oldintlength)
210 {
211
212 /* This will keep the integer part from overflowing past the array. */
213 bnlength = oldbnlength - oldintlength + min(oldintlength, newintlength);
214
215 _fmemcpy(new+newbnlength-newintlength-oldbnlength+oldintlength,
216 old,bnlength);
217 }
218 else
219 {
220 bnlength = newbnlength - newintlength + min(oldintlength, newintlength);
221 _fmemcpy(new,old+oldbnlength-oldintlength-newbnlength+newintlength,
222 bnlength);
223 }
224 intlength = saveintlength;
225 bnlength = savebnlength;
226 return(0);
227 }
228
229 /********************************************************************/
230 /* bn_hexdump() - for debugging, dumps to stdout */
231
bn_hexdump(bn_t r)232 void bn_hexdump(bn_t r)
233 {
234 int i;
235
236 for (i=0; i<bnlength; i++)
237 printf("%02X ", r[i]);
238 printf("\n");
239 return;
240 }
241
242 /**********************************************************************/
243 /* strtobn() - converts a string into a bignumer */
244 /* r - pointer to a bignumber */
245 /* s - string in the floating point format [-][digits].[digits] */
246 /* note: the string may not be empty or have extra space and may */
247 /* not use scientific notation (2.3e4). */
248
strtobn(bn_t r,char * s)249 bn_t strtobn(bn_t r, char *s)
250 {
251 unsigned l;
252 bn_t onesbyte;
253 int signflag=0;
254 long longval;
255
256 clear_bn(r);
257 onesbyte = r + bnlength - intlength;
258
259 if (s[0] == '+') /* for + sign */
260 {
261 s++;
262 }
263 else if (s[0] == '-') /* for neg sign */
264 {
265 signflag = 1;
266 s++;
267 }
268
269 if (strchr(s, '.') != NULL) /* is there a decimal point? */
270 {
271 l = strlen(s) - 1; /* start with the last digit */
272 while (s[l] >= '0' && s[l] <= '9') /* while a digit */
273 {
274 *onesbyte = (BYTE)(s[l--] - '0');
275 div_a_bn_int(r, 10);
276 }
277
278 if (s[l] == '.')
279 {
280 longval = atol(s);
281 switch (intlength)
282 { /* only 1, 2, or 4 are allowed */
283 case 1:
284 *onesbyte = (BYTE)longval;
285 break;
286 case 2:
287 big_set16(onesbyte, (U16)longval);
288 break;
289 case 4:
290 big_set32(onesbyte, longval);
291 break;
292 }
293 }
294 }
295 else
296 {
297 longval = atol(s);
298 switch (intlength)
299 { /* only 1, 2, or 4 are allowed */
300 case 1:
301 *onesbyte = (BYTE)longval;
302 break;
303 case 2:
304 big_set16(onesbyte, (U16)longval);
305 break;
306 case 4:
307 big_set32(onesbyte, longval);
308 break;
309 }
310 }
311
312
313 if (signflag)
314 neg_a_bn(r);
315
316 return r;
317 }
318
319 /********************************************************************/
320 /* strlen_needed() - returns string length needed to hold bignumber */
321
strlen_needed()322 int strlen_needed()
323 {
324 int length = 3;
325
326 /* first space for integer part */
327 switch(intlength)
328 {
329 case 1:
330 length = 3; /* max 127 */
331 break;
332 case 2:
333 length = 5; /* max 32767 */
334 break;
335 case 4:
336 length = 10; /* max 2147483647 */
337 break;
338 }
339 length += decimals; /* decimal part */
340 length += 2; /* decimal point and sign */
341 length += 4; /* null and a little extra for safety */
342 return(length);
343 }
344
345 /********************************************************************/
346 /* bntostr() - converts a bignumber into a string */
347 /* s - string, must be large enough to hold the number. */
348 /* r - bignumber */
349 /* will covert to a floating point notation */
350 /* SIDE-EFFECT: the bignumber, r, is destroyed. */
351 /* Copy it first if necessary. */
352
unsafe_bntostr(char * s,int dec,bn_t r)353 char *unsafe_bntostr(char *s, int dec, bn_t r)
354 {
355 int l=0, d;
356 bn_t onesbyte;
357 long longval = 0;
358
359 if (dec == 0)
360 dec = decimals;
361 onesbyte = r + bnlength - intlength;
362
363 if (is_bn_neg(r))
364 {
365 neg_a_bn(r);
366 *(s++) = '-';
367 }
368 switch (intlength)
369 { /* only 1, 2, or 4 are allowed */
370 case 1:
371 longval = *onesbyte;
372 break;
373 case 2:
374 longval = big_access16(onesbyte);
375 break;
376 case 4:
377 longval = big_access32(onesbyte);
378 break;
379 }
380 ltoa(longval, s, 10);
381 l = strlen(s);
382 s[l++] = '.';
383 for (d=0; d < dec; d++)
384 {
385 *onesbyte = 0; /* clear out highest byte */
386 mult_a_bn_int(r, 10);
387 if (is_bn_zero(r))
388 break;
389 s[l++] = (BYTE)(*onesbyte + '0');
390 }
391 s[l] = '\0'; /* don't forget nul char */
392
393 return s;
394 }
395
396 #if (_MSC_VER >= 700)
397 #pragma code_seg ( ) /* back to normal segment */
398 #endif
399
400 /*********************************************************************/
401 /* b = l */
402 /* Converts a long to a bignumber */
inttobn(bn_t r,long longval)403 bn_t inttobn(bn_t r, long longval)
404 {
405 bn_t onesbyte;
406
407 clear_bn(r);
408 onesbyte = r + bnlength - intlength;
409 switch (intlength)
410 { /* only 1, 2, or 4 are allowed */
411 case 1:
412 *onesbyte = (BYTE)longval;
413 break;
414 case 2:
415 big_set16(onesbyte, (U16)longval);
416 break;
417 case 4:
418 big_set32(onesbyte, longval);
419 break;
420 }
421 return r;
422 }
423
424 /*********************************************************************/
425 /* l = floor(b), floor rounds down */
426 /* Converts the integer part a bignumber to a long */
bntoint(bn_t n)427 long bntoint(bn_t n)
428 {
429 bn_t onesbyte;
430 long longval = 0;
431
432 onesbyte = n + bnlength - intlength;
433 switch (intlength)
434 { /* only 1, 2, or 4 are allowed */
435 case 1:
436 longval = *onesbyte;
437 break;
438 case 2:
439 longval = big_access16(onesbyte);
440 break;
441 case 4:
442 longval = big_access32(onesbyte);
443 break;
444 }
445 return longval;
446 }
447
448 /*********************************************************************/
449 /* b = f */
450 /* Converts a double to a bignumber */
floattobn(bn_t r,LDBL f)451 bn_t floattobn(bn_t r, LDBL f)
452 {
453 #ifndef USE_BIGNUM_C_CODE
454 /* Only use this when using the ASM code as the C version of */
455 /* floattobf() calls floattobn(), an infinite recursive loop. */
456 floattobf(bftmp1, f);
457 return bftobn(r, bftmp1);
458
459 #else
460 bn_t onesbyte;
461 int i;
462 int signflag=0;
463
464 clear_bn(r);
465 onesbyte = r + bnlength - intlength;
466
467 if (f < 0)
468 {
469 signflag = 1;
470 f = -f;
471 }
472
473 switch (intlength)
474 { /* only 1, 2, or 4 are allowed */
475 case 1:
476 *onesbyte = (BYTE)f;
477 break;
478 case 2:
479 big_set16(onesbyte, (U16)f);
480 break;
481 case 4:
482 big_set32(onesbyte, (U32)f);
483 break;
484 }
485
486 f -= (long)f; /* keep only the decimal part */
487 for (i = bnlength - intlength - 1; i >= 0 && f != 0.0; i--)
488 {
489 f *= 256;
490 r[i] = (BYTE)f; /* keep use the integer part */
491 f -= (BYTE)f; /* now throw away the integer part */
492 }
493
494 if (signflag)
495 neg_a_bn(r);
496
497 return r;
498 #endif /* USE_BIGNUM_C_CODE */
499 }
500
501 /********************************************************************/
502 /* sign(r) */
sign_bn(bn_t n)503 int sign_bn(bn_t n)
504 {
505 return is_bn_neg(n) ? -1 : is_bn_not_zero(n) ? 1 : 0;
506 }
507
508 /********************************************************************/
509 /* r = |n| */
abs_bn(bn_t r,bn_t n)510 bn_t abs_bn(bn_t r, bn_t n)
511 {
512 copy_bn(r,n);
513 if (is_bn_neg(r))
514 neg_a_bn(r);
515 return r;
516 }
517
518 /********************************************************************/
519 /* r = |r| */
abs_a_bn(bn_t r)520 bn_t abs_a_bn(bn_t r)
521 {
522 if (is_bn_neg(r))
523 neg_a_bn(r);
524 return r;
525 }
526
527 /********************************************************************/
528 /* r = 1/n */
529 /* uses bntmp1 - bntmp3 - global temp bignumbers */
530 /* SIDE-EFFECTS: */
531 /* n ends up as |n| Make copy first if necessary. */
unsafe_inv_bn(bn_t r,bn_t n)532 bn_t unsafe_inv_bn(bn_t r, bn_t n)
533 {
534 int signflag=0, i;
535 long maxval;
536 LDBL f;
537 bn_t orig_r, orig_n; /* orig_bntmp1 not needed here */
538 int orig_bnlength,
539 orig_padding,
540 orig_rlength,
541 orig_shiftfactor;
542
543 /* use Newton's recursive method for zeroing in on 1/n : r=r(2-rn) */
544
545 if (is_bn_neg(n))
546 { /* will be a lot easier to deal with just positives */
547 signflag = 1;
548 neg_a_bn(n);
549 }
550
551 f = bntofloat(n);
552 if (f == 0) /* division by zero */
553 {
554 max_bn(r);
555 return r;
556 }
557 f = 1/f; /* approximate inverse */
558 maxval = (1L << ((intlength<<3)-1)) - 1;
559 if (f > maxval) /* check for overflow */
560 {
561 max_bn(r);
562 return r;
563 }
564 else if (f <= -maxval)
565 {
566 max_bn(r);
567 neg_a_bn(r);
568 return r;
569 }
570
571 /* With Newton's Method, there is no need to calculate all the digits */
572 /* every time. The precision approximately doubles each iteration. */
573 /* Save original values. */
574 orig_bnlength = bnlength;
575 orig_padding = padding;
576 orig_rlength = rlength;
577 orig_shiftfactor = shiftfactor;
578 orig_r = r;
579 orig_n = n;
580 /* orig_bntmp1 = bntmp1; */
581
582 /* calculate new starting values */
583 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
584 if (bnlength > orig_bnlength)
585 bnlength = orig_bnlength;
586 calc_lengths();
587
588 /* adjust pointers */
589 r = orig_r + orig_bnlength - bnlength;
590 n = orig_n + orig_bnlength - bnlength;
591 /* bntmp1 = orig_bntmp1 + orig_bnlength - bnlength; */
592
593 floattobn(r, f); /* start with approximate inverse */
594 clear_bn(bntmp2); /* will be used as 1.0 and 2.0 */
595
596 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
597 {
598 /* adjust lengths */
599 bnlength <<= 1; /* double precision */
600 if (bnlength > orig_bnlength)
601 bnlength = orig_bnlength;
602 calc_lengths();
603 r = orig_r + orig_bnlength - bnlength;
604 n = orig_n + orig_bnlength - bnlength;
605 /* bntmp1 = orig_bntmp1 + orig_bnlength - bnlength; */
606
607 unsafe_mult_bn(bntmp1, r, n); /* bntmp1=rn */
608 inttobn(bntmp2, 1); /* bntmp2 = 1.0 */
609 if (bnlength == orig_bnlength && cmp_bn(bntmp2, bntmp1+shiftfactor) == 0 ) /* if not different */
610 break; /* they must be the same */
611 inttobn(bntmp2, 2); /* bntmp2 = 2.0 */
612 sub_bn(bntmp3, bntmp2, bntmp1+shiftfactor); /* bntmp3=2-rn */
613 unsafe_mult_bn(bntmp1, r, bntmp3); /* bntmp1=r(2-rn) */
614 copy_bn(r, bntmp1+shiftfactor); /* r = bntmp1 */
615 }
616
617 /* restore original values */
618 bnlength = orig_bnlength;
619 padding = orig_padding;
620 rlength = orig_rlength;
621 shiftfactor = orig_shiftfactor;
622 r = orig_r;
623 n = orig_n;
624 /* bntmp1 = orig_bntmp1; */
625
626 if (signflag)
627 {
628 neg_a_bn(r);
629 }
630 return r;
631 }
632
633 /********************************************************************/
634 /* r = n1/n2 */
635 /* r - result of length bnlength */
636 /* uses bntmp1 - bntmp3 - global temp bignumbers */
637 /* SIDE-EFFECTS: */
638 /* n1, n2 can end up as GARBAGE */
639 /* Make copies first if necessary. */
unsafe_div_bn(bn_t r,bn_t n1,bn_t n2)640 bn_t unsafe_div_bn(bn_t r, bn_t n1, bn_t n2)
641 {
642 int scale1, scale2, scale, sign=0, i;
643 long maxval;
644 LDBL a, b, f;
645
646 /* first, check for valid data */
647 a = bntofloat(n1);
648 if (a == 0) /* division into zero */
649 {
650 clear_bn(r); /* return 0 */
651 return r;
652 }
653 b = bntofloat(n2);
654 if (b == 0) /* division by zero */
655 {
656 max_bn(r);
657 return r;
658 }
659 f = a/b; /* approximate quotient */
660 maxval = (1L << ((intlength<<3)-1)) - 1;
661 if (f > maxval) /* check for overflow */
662 {
663 max_bn(r);
664 return r;
665 }
666 else if (f <= -maxval)
667 {
668 max_bn(r);
669 neg_a_bn(r);
670 return r;
671 }
672 /* appears to be ok, do division */
673
674 if (is_bn_neg(n1))
675 {
676 neg_a_bn(n1);
677 sign = !sign;
678 }
679 if (is_bn_neg(n2))
680 {
681 neg_a_bn(n2);
682 sign = !sign;
683 }
684
685 /* scale n1 and n2 so: |n| >= 1/256 */
686 /* scale = (int)(log(1/fabs(a))/LOG_256) = LOG_256(1/|a|) */
687 i = bnlength-1;
688 while (i >= 0 && n1[i] == 0)
689 i--;
690 scale1 = bnlength - i - 2;
691 if (scale1 < 0)
692 scale1 = 0;
693 i = bnlength-1;
694 while (i >= 0 && n2[i] == 0)
695 i--;
696 scale2 = bnlength - i - 2;
697 if (scale2 < 0)
698 scale2 = 0;
699
700 /* shift n1, n2 */
701 /* important!, use memmove(), not memcpy() */
702 _fmemmove(n1+scale1, n1, bnlength-scale1); /* shift bytes over */
703 _fmemset(n1, 0, scale1); /* zero out the rest */
704 _fmemmove(n2+scale2, n2, bnlength-scale2); /* shift bytes over */
705 _fmemset(n2, 0, scale2); /* zero out the rest */
706
707 unsafe_inv_bn(r, n2);
708 unsafe_mult_bn(bntmp1, n1, r);
709 copy_bn(r, bntmp1+shiftfactor); /* r = bntmp1 */
710
711 if (scale1 != scale2)
712 {
713 /* Rescale r back to what it should be. Overflow has already been checked */
714 if (scale1 > scale2) /* answer is too big, adjust it */
715 {
716 scale = scale1-scale2;
717 _fmemmove(r, r+scale, bnlength-scale); /* shift bytes over */
718 _fmemset(r+bnlength-scale, 0, scale); /* zero out the rest */
719 }
720 else if (scale1 < scale2) /* answer is too small, adjust it */
721 {
722 scale = scale2-scale1;
723 _fmemmove(r+scale, r, bnlength-scale); /* shift bytes over */
724 _fmemset(r, 0, scale); /* zero out the rest */
725 }
726 /* else scale1 == scale2 */
727
728 }
729
730 if (sign)
731 neg_a_bn(r);
732
733 return r;
734 }
735
736 /********************************************************************/
737 /* sqrt(r) */
738 /* uses bntmp1 - bntmp6 - global temp bignumbers */
739 /* SIDE-EFFECTS: */
740 /* n ends up as |n| */
sqrt_bn(bn_t r,bn_t n)741 bn_t sqrt_bn(bn_t r, bn_t n)
742 {
743 int i, comp, almost_match=0;
744 LDBL f;
745 bn_t orig_r, orig_n;
746 int orig_bnlength,
747 orig_padding,
748 orig_rlength,
749 orig_shiftfactor;
750
751 /* use Newton's recursive method for zeroing in on sqrt(n): r=.5(r+n/r) */
752
753 if (is_bn_neg(n))
754 { /* sqrt of a neg, return 0 */
755 clear_bn(r);
756 return r;
757 }
758
759 f = bntofloat(n);
760 if (f == 0) /* division by zero will occur */
761 {
762 clear_bn(r); /* sqrt(0) = 0 */
763 return r;
764 }
765 f = sqrtl(f); /* approximate square root */
766 /* no need to check overflow */
767
768 /* With Newton's Method, there is no need to calculate all the digits */
769 /* every time. The precision approximately doubles each iteration. */
770 /* Save original values. */
771 orig_bnlength = bnlength;
772 orig_padding = padding;
773 orig_rlength = rlength;
774 orig_shiftfactor = shiftfactor;
775 orig_r = r;
776 orig_n = n;
777
778 /* calculate new starting values */
779 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
780 if (bnlength > orig_bnlength)
781 bnlength = orig_bnlength;
782 calc_lengths();
783
784 /* adjust pointers */
785 r = orig_r + orig_bnlength - bnlength;
786 n = orig_n + orig_bnlength - bnlength;
787
788 floattobn(r, f); /* start with approximate sqrt */
789 copy_bn(bntmp4, r);
790
791 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
792 {
793 /* adjust lengths */
794 bnlength <<= 1; /* double precision */
795 if (bnlength > orig_bnlength)
796 bnlength = orig_bnlength;
797 calc_lengths();
798 r = orig_r + orig_bnlength - bnlength;
799 n = orig_n + orig_bnlength - bnlength;
800
801 copy_bn(bntmp6, r);
802 copy_bn(bntmp5, n);
803 unsafe_div_bn(bntmp4, bntmp5, bntmp6);
804 add_a_bn(r, bntmp4);
805 half_a_bn(r);
806 if (bnlength == orig_bnlength && (comp=abs(cmp_bn(r, bntmp4))) < 8 ) /* if match or almost match */
807 {
808 if (comp < 4 /* perfect or near perfect match */
809 || almost_match == 1) /* close enough for 2nd time */
810 break;
811 else /* this is the first time they almost matched */
812 almost_match++;
813 }
814 }
815
816 /* restore original values */
817 bnlength = orig_bnlength;
818 padding = orig_padding;
819 rlength = orig_rlength;
820 shiftfactor = orig_shiftfactor;
821 r = orig_r;
822 n = orig_n;
823
824 return r;
825 }
826
827 /********************************************************************/
828 /* exp(r) */
829 /* uses bntmp1, bntmp2, bntmp3 - global temp bignumbers */
exp_bn(bn_t r,bn_t n)830 bn_t exp_bn(bn_t r, bn_t n)
831 {
832 U16 fact=1;
833
834 if (is_bn_zero(n))
835 {
836 inttobn(r, 1);
837 return r;
838 }
839
840 /* use Taylor Series (very slow convergence) */
841 inttobn(r, 1); /* start with r=1.0 */
842 copy_bn(bntmp2, r);
843 for (;;)
844 {
845 /* copy n, if n is negative, mult_bn() alters n */
846 unsafe_mult_bn(bntmp3, bntmp2, copy_bn(bntmp1, n));
847 copy_bn(bntmp2, bntmp3+shiftfactor);
848 div_a_bn_int(bntmp2, fact);
849 if (!is_bn_not_zero(bntmp2))
850 break; /* too small to register */
851 add_a_bn(r, bntmp2);
852 fact++;
853 }
854 return r;
855 }
856
857 /********************************************************************/
858 /* ln(r) */
859 /* uses bntmp1 - bntmp6 - global temp bignumbers */
860 /* SIDE-EFFECTS: */
861 /* n ends up as |n| */
unsafe_ln_bn(bn_t r,bn_t n)862 bn_t unsafe_ln_bn(bn_t r, bn_t n)
863 {
864 int i, comp, almost_match=0;
865 long maxval;
866 LDBL f;
867 bn_t orig_r, orig_n, orig_bntmp5, orig_bntmp4;
868 int orig_bnlength,
869 orig_padding,
870 orig_rlength,
871 orig_shiftfactor;
872
873 /* use Newton's recursive method for zeroing in on ln(n): r=r+n*exp(-r)-1 */
874
875 if (is_bn_neg(n) || is_bn_zero(n))
876 { /* error, return largest neg value */
877 max_bn(r);
878 neg_a_bn(r);
879 return r;
880 }
881
882 f = bntofloat(n);
883 f = logl(f); /* approximate ln(x) */
884 maxval = (1L << ((intlength<<3)-1)) - 1;
885 if (f > maxval) /* check for overflow */
886 {
887 max_bn(r);
888 return r;
889 }
890 else if (f <= -maxval)
891 {
892 max_bn(r);
893 neg_a_bn(r);
894 return r;
895 }
896 /* appears to be ok, do ln */
897
898 /* With Newton's Method, there is no need to calculate all the digits */
899 /* every time. The precision approximately doubles each iteration. */
900 /* Save original values. */
901 orig_bnlength = bnlength;
902 orig_padding = padding;
903 orig_rlength = rlength;
904 orig_shiftfactor = shiftfactor;
905 orig_r = r;
906 orig_n = n;
907 orig_bntmp5 = bntmp5;
908 orig_bntmp4 = bntmp4;
909
910 inttobn(bntmp4, 1); /* set before setting new values */
911
912 /* calculate new starting values */
913 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
914 if (bnlength > orig_bnlength)
915 bnlength = orig_bnlength;
916 calc_lengths();
917
918 /* adjust pointers */
919 r = orig_r + orig_bnlength - bnlength;
920 n = orig_n + orig_bnlength - bnlength;
921 bntmp5 = orig_bntmp5 + orig_bnlength - bnlength;
922 bntmp4 = orig_bntmp4 + orig_bnlength - bnlength;
923
924 floattobn(r, f); /* start with approximate ln */
925 neg_a_bn(r); /* -r */
926 copy_bn(bntmp5, r); /* -r */
927
928 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
929 {
930 /* adjust lengths */
931 bnlength <<= 1; /* double precision */
932 if (bnlength > orig_bnlength)
933 bnlength = orig_bnlength;
934 calc_lengths();
935 r = orig_r + orig_bnlength - bnlength;
936 n = orig_n + orig_bnlength - bnlength;
937 bntmp5 = orig_bntmp5 + orig_bnlength - bnlength;
938 bntmp4 = orig_bntmp4 + orig_bnlength - bnlength;
939 exp_bn(bntmp6, r); /* exp(-r) */
940 unsafe_mult_bn(bntmp2, bntmp6, n); /* n*exp(-r) */
941 sub_a_bn(bntmp2+shiftfactor, bntmp4); /* n*exp(-r) - 1 */
942 sub_a_bn(r, bntmp2+shiftfactor); /* -r - (n*exp(-r) - 1) */
943
944 if (bnlength == orig_bnlength && (comp=abs(cmp_bn(r, bntmp5))) < 8 ) /* if match or almost match */
945 {
946 if (comp < 4 /* perfect or near perfect match */
947 || almost_match == 1) /* close enough for 2nd time */
948 break;
949 else /* this is the first time they almost matched */
950 almost_match++;
951 }
952 copy_bn(bntmp5, r); /* -r */
953 }
954
955 /* restore original values */
956 bnlength = orig_bnlength;
957 padding = orig_padding;
958 rlength = orig_rlength;
959 shiftfactor = orig_shiftfactor;
960 r = orig_r;
961 n = orig_n;
962 bntmp5 = orig_bntmp5;
963 bntmp4 = orig_bntmp4;
964
965 neg_a_bn(r); /* -(-r) */
966 return r;
967 }
968
969 /********************************************************************/
970 /* sincos_bn(r) */
971 /* uses bntmp1 - bntmp2 - global temp bignumbers */
972 /* SIDE-EFFECTS: */
973 /* n ends up as |n| mod (pi/4) */
unsafe_sincos_bn(bn_t s,bn_t c,bn_t n)974 bn_t unsafe_sincos_bn(bn_t s, bn_t c, bn_t n)
975 {
976 U16 fact=2;
977 int k=0, i, halves;
978 int signcos=0, signsin=0, switch_sincos=0;
979
980 #ifndef CALCULATING_BIG_PI
981 /* assure range 0 <= x < pi/4 */
982
983 if (is_bn_zero(n))
984 {
985 clear_bn(s); /* sin(0) = 0 */
986 inttobn(c, 1); /* cos(0) = 1 */
987 return s;
988 }
989
990 if (is_bn_neg(n))
991 {
992 signsin = !signsin; /* sin(-x) = -sin(x), odd; cos(-x) = cos(x), even */
993 neg_a_bn(n);
994 }
995 /* n >= 0 */
996
997 double_bn(bntmp1, bn_pi); /* 2*pi */
998 /* this could be done with remainders, but it would probably be slower */
999 while (cmp_bn(n, bntmp1) >= 0) /* while n >= 2*pi */
1000 sub_a_bn(n, bntmp1);
1001 /* 0 <= n < 2*pi */
1002
1003 copy_bn(bntmp1, bn_pi); /* pi */
1004 if (cmp_bn(n, bntmp1) >= 0) /* if n >= pi */
1005 {
1006 sub_a_bn(n, bntmp1);
1007 signsin = !signsin;
1008 signcos = !signcos;
1009 }
1010 /* 0 <= n < pi */
1011
1012 half_bn(bntmp1, bn_pi); /* pi/2 */
1013 if (cmp_bn(n, bntmp1) > 0) /* if n > pi/2 */
1014 {
1015 sub_bn(n, bn_pi, n); /* pi - n */
1016 signcos = !signcos;
1017 }
1018 /* 0 <= n < pi/2 */
1019
1020 half_bn(bntmp1, bn_pi); /* pi/2 */
1021 half_a_bn(bntmp1); /* pi/4 */
1022 if (cmp_bn(n, bntmp1) > 0) /* if n > pi/4 */
1023 {
1024 half_bn(bntmp1, bn_pi); /* pi/2 */
1025 sub_bn(n, bntmp1, n); /* pi/2 - n */
1026 switch_sincos = !switch_sincos;
1027 }
1028 /* 0 <= n < pi/4 */
1029
1030 /* this looks redundant, but n could now be zero when it wasn't before */
1031 if (is_bn_zero(n))
1032 {
1033 clear_bn(s); /* sin(0) = 0 */
1034 inttobn(c, 1); /* cos(0) = 1 */
1035 return s;
1036 }
1037
1038
1039 /* at this point, the double angle trig identities could be used as many */
1040 /* times as desired to reduce the range to pi/8, pi/16, etc... Each time */
1041 /* the range is cut in half, the number of iterations required is reduced */
1042 /* by "quite a bit." It's just a matter of testing to see what gives the */
1043 /* optimal results. */
1044 /* halves = bnlength / 10; */ /* this is experimental */
1045 halves = 1;
1046 for (i = 0; i < halves; i++)
1047 half_a_bn(n);
1048 #endif
1049
1050 /* use Taylor Series (very slow convergence) */
1051 copy_bn(s, n); /* start with s=n */
1052 inttobn(c, 1); /* start with c=1 */
1053 copy_bn(bntmp1, n); /* the current x^n/n! */
1054
1055 for (;;)
1056 {
1057 /* even terms for cosine */
1058 unsafe_mult_bn(bntmp2, bntmp1, n);
1059 copy_bn(bntmp1, bntmp2+shiftfactor);
1060 div_a_bn_int(bntmp1, fact++);
1061 if (!is_bn_not_zero(bntmp1))
1062 break; /* too small to register */
1063 if (k) /* alternate between adding and subtracting */
1064 add_a_bn(c, bntmp1);
1065 else
1066 sub_a_bn(c, bntmp1);
1067
1068 /* odd terms for sine */
1069 unsafe_mult_bn(bntmp2, bntmp1, n);
1070 copy_bn(bntmp1, bntmp2+shiftfactor);
1071 div_a_bn_int(bntmp1, fact++);
1072 if (!is_bn_not_zero(bntmp1))
1073 break; /* too small to register */
1074 if (k) /* alternate between adding and subtracting */
1075 add_a_bn(s, bntmp1);
1076 else
1077 sub_a_bn(s, bntmp1);
1078 k = !k; /* toggle */
1079 #ifdef CALCULATING_BIG_PI
1080 printf("."); /* lets you know it's doing something */
1081 #endif
1082 }
1083
1084 #ifndef CALCULATING_BIG_PI
1085 /* now need to undo what was done by cutting angles in half */
1086 inttobn(bntmp1, 1);
1087 for (i = 0; i < halves; i++)
1088 {
1089 unsafe_mult_bn(bntmp2, s, c); /* no need for safe mult */
1090 double_bn(s, bntmp2+shiftfactor); /* sin(2x) = 2*sin(x)*cos(x) */
1091 unsafe_square_bn(bntmp2,c);
1092 double_a_bn(bntmp2+shiftfactor);
1093 sub_bn(c, bntmp2+shiftfactor, bntmp1); /* cos(2x) = 2*cos(x)*cos(x) - 1 */
1094 }
1095
1096 if (switch_sincos)
1097 {
1098 copy_bn(bntmp1, s);
1099 copy_bn(s, c);
1100 copy_bn(c, bntmp1);
1101 }
1102 if (signsin)
1103 neg_a_bn(s);
1104 if (signcos)
1105 neg_a_bn(c);
1106 #endif
1107
1108 return s; /* return sine I guess */
1109 }
1110
1111 /********************************************************************/
1112 /* atan(r) */
1113 /* uses bntmp1 - bntmp5 - global temp bignumbers */
1114 /* SIDE-EFFECTS: */
1115 /* n ends up as |n| or 1/|n| */
unsafe_atan_bn(bn_t r,bn_t n)1116 bn_t unsafe_atan_bn(bn_t r, bn_t n)
1117 {
1118 int i, comp, almost_match=0, signflag=0;
1119 LDBL f;
1120 bn_t orig_r, orig_n, orig_bn_pi, orig_bntmp3;
1121 int orig_bnlength,
1122 orig_padding,
1123 orig_rlength,
1124 orig_shiftfactor;
1125 int large_arg;
1126
1127 /* use Newton's recursive method for zeroing in on atan(n): r=r-cos(r)(sin(r)-n*cos(r)) */
1128
1129 if (is_bn_neg(n))
1130 {
1131 signflag = 1;
1132 neg_a_bn(n);
1133 }
1134
1135 /* If n is very large, atanl() won't give enough decimal places to be a */
1136 /* good enough initial guess for Newton's Method. If it is larger than */
1137 /* say, 1, atan(n) = pi/2 - acot(n) = pi/2 - atan(1/n). */
1138
1139 f = bntofloat(n);
1140 large_arg = f > 1.0;
1141 if (large_arg)
1142 {
1143 unsafe_inv_bn(bntmp4, n); /* need to use bntmp4 here since inv uses bntmp3 */
1144 copy_bn(n, bntmp4);
1145 f = bntofloat(n);
1146 }
1147
1148 clear_bn(bntmp3); /* not really necessary, but makes things more consistent */
1149
1150 /* With Newton's Method, there is no need to calculate all the digits */
1151 /* every time. The precision approximately doubles each iteration. */
1152 /* Save original values. */
1153 orig_bnlength = bnlength;
1154 orig_padding = padding;
1155 orig_rlength = rlength;
1156 orig_shiftfactor = shiftfactor;
1157 orig_bn_pi = bn_pi;
1158 orig_r = r;
1159 orig_n = n;
1160 orig_bntmp3 = bntmp3;
1161
1162 /* calculate new starting values */
1163 bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
1164 if (bnlength > orig_bnlength)
1165 bnlength = orig_bnlength;
1166 calc_lengths();
1167
1168 /* adjust pointers */
1169 r = orig_r + orig_bnlength - bnlength;
1170 n = orig_n + orig_bnlength - bnlength;
1171 bn_pi = orig_bn_pi + orig_bnlength - bnlength;
1172 bntmp3 = orig_bntmp3 + orig_bnlength - bnlength;
1173
1174 f = atanl(f); /* approximate arctangent */
1175 /* no need to check overflow */
1176
1177 floattobn(r, f); /* start with approximate atan */
1178 copy_bn(bntmp3, r);
1179
1180 for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
1181 {
1182 /* adjust lengths */
1183 bnlength <<= 1; /* double precision */
1184 if (bnlength > orig_bnlength)
1185 bnlength = orig_bnlength;
1186 calc_lengths();
1187 r = orig_r + orig_bnlength - bnlength;
1188 n = orig_n + orig_bnlength - bnlength;
1189 bn_pi = orig_bn_pi + orig_bnlength - bnlength;
1190 bntmp3 = orig_bntmp3 + orig_bnlength - bnlength;
1191
1192 #ifdef CALCULATING_BIG_PI
1193 printf("\natan() loop #%i, bnlength=%i\nsincos() loops\n", i, bnlength);
1194 #endif
1195 unsafe_sincos_bn(bntmp4, bntmp5, bntmp3); /* sin(r), cos(r) */
1196 copy_bn(bntmp3, r); /* restore bntmp3 from sincos_bn() */
1197 copy_bn(bntmp1, bntmp5);
1198 unsafe_mult_bn(bntmp2, n, bntmp1); /* n*cos(r) */
1199 sub_a_bn(bntmp4, bntmp2+shiftfactor); /* sin(r) - n*cos(r) */
1200 unsafe_mult_bn(bntmp1, bntmp5, bntmp4); /* cos(r) * (sin(r) - n*cos(r)) */
1201 sub_a_bn(r, bntmp1+shiftfactor); /* r - cos(r) * (sin(r) - n*cos(r)) */
1202
1203 #ifdef CALCULATING_BIG_PI
1204 putchar('\n');
1205 bn_hexdump(r);
1206 #endif
1207 if (bnlength == orig_bnlength && (comp=abs(cmp_bn(r, bntmp3))) < 8 ) /* if match or almost match */
1208 {
1209 #ifdef CALCULATING_BIG_PI
1210 printf("atan() loop comp=%i\n", comp);
1211 #endif
1212 if (comp < 4 /* perfect or near perfect match */
1213 || almost_match == 1) /* close enough for 2nd time */
1214 break;
1215 else /* this is the first time they almost matched */
1216 almost_match++;
1217 }
1218
1219 #ifdef CALCULATING_BIG_PI
1220 if (bnlength == orig_bnlength && comp >= 8)
1221 printf("atan() loop comp=%i\n", comp);
1222 #endif
1223
1224 /* the following gets set above and not altered before we use it! */
1225 /* copy_bn(bntmp3, r); */ /* make a copy for later comparison */
1226 }
1227
1228 /* restore original values */
1229 bnlength = orig_bnlength;
1230 padding = orig_padding;
1231 rlength = orig_rlength;
1232 shiftfactor = orig_shiftfactor;
1233 bn_pi = orig_bn_pi;
1234 r = orig_r;
1235 n = orig_n;
1236 bntmp3 = orig_bntmp3;
1237
1238 if (large_arg)
1239 {
1240 half_bn(bntmp3, bn_pi); /* pi/2 */
1241 sub_a_bn(bntmp3, r); /* pi/2 - atan(1/n) */
1242 copy_bn(r, bntmp3);
1243 }
1244
1245 if (signflag)
1246 neg_a_bn(r);
1247 return r;
1248 }
1249
1250 /********************************************************************/
1251 /* atan2(r,ny,nx) */
1252 /* uses bntmp1 - bntmp6 - global temp bigfloats */
unsafe_atan2_bn(bn_t r,bn_t ny,bn_t nx)1253 bn_t unsafe_atan2_bn(bn_t r, bn_t ny, bn_t nx)
1254 {
1255 int signx, signy;
1256
1257 signx = sign_bn(nx);
1258 signy = sign_bn(ny);
1259
1260 if (signy == 0)
1261 {
1262 if (signx < 0)
1263 copy_bn(r, bn_pi); /* negative x axis, 180 deg */
1264 else /* signx >= 0 positive x axis, 0 */
1265 clear_bn(r);
1266 return(r);
1267 }
1268 if (signx == 0)
1269 {
1270 copy_bn(r, bn_pi); /* y axis */
1271 half_a_bn(r); /* +90 deg */
1272 if (signy < 0)
1273 neg_a_bn(r); /* -90 deg */
1274 return(r);
1275 }
1276
1277 if (signy < 0)
1278 neg_a_bn(ny);
1279 if (signx < 0)
1280 neg_a_bn(nx);
1281 unsafe_div_bn(bntmp6,ny,nx);
1282 unsafe_atan_bn(r, bntmp6);
1283 if (signx < 0)
1284 sub_bn(r,bn_pi,r);
1285 if(signy < 0)
1286 neg_a_bn(r);
1287 return(r);
1288 }
1289
1290
1291 /**********************************************************************/
1292 /* The rest of the functions are "safe" versions of the routines that */
1293 /* have side effects which alter the parameters */
1294 /**********************************************************************/
1295
1296 /**********************************************************************/
full_mult_bn(bn_t r,bn_t n1,bn_t n2)1297 bn_t full_mult_bn(bn_t r, bn_t n1, bn_t n2)
1298 {
1299 int sign1, sign2;
1300
1301 sign1 = is_bn_neg(n1);
1302 sign2 = is_bn_neg(n2);
1303 unsafe_full_mult_bn(r, n1, n2);
1304 if (sign1)
1305 neg_a_bn(n1);
1306 if (sign2)
1307 neg_a_bn(n2);
1308 return r;
1309 }
1310
1311 /**********************************************************************/
mult_bn(bn_t r,bn_t n1,bn_t n2)1312 bn_t mult_bn(bn_t r, bn_t n1, bn_t n2)
1313 {
1314 int sign1, sign2;
1315
1316 /* TW ENDFIX */
1317 sign1 = is_bn_neg(n1);
1318 sign2 = is_bn_neg(n2);
1319 unsafe_mult_bn(r, n1, n2);
1320 if (sign1)
1321 neg_a_bn(n1);
1322 if (sign2)
1323 neg_a_bn(n2);
1324 return r;
1325 }
1326
1327 /**********************************************************************/
full_square_bn(bn_t r,bn_t n)1328 bn_t full_square_bn(bn_t r, bn_t n)
1329 {
1330 int sign;
1331
1332 sign = is_bn_neg(n);
1333 unsafe_full_square_bn(r, n);
1334 if (sign)
1335 neg_a_bn(n);
1336 return r;
1337 }
1338
1339 /**********************************************************************/
square_bn(bn_t r,bn_t n)1340 bn_t square_bn(bn_t r, bn_t n)
1341 {
1342 int sign;
1343
1344 sign = is_bn_neg(n);
1345 unsafe_square_bn(r, n);
1346 if (sign)
1347 neg_a_bn(n);
1348 return r;
1349 }
1350
1351 /**********************************************************************/
div_bn_int(bn_t r,bn_t n,U16 u)1352 bn_t div_bn_int(bn_t r, bn_t n, U16 u)
1353 {
1354 int sign;
1355
1356 sign = is_bn_neg(n);
1357 unsafe_div_bn_int(r, n, u);
1358 if (sign)
1359 neg_a_bn(n);
1360 return r;
1361 }
1362
1363 #if (_MSC_VER >= 700)
1364 #pragma code_seg ("bignum1_text") /* place following in an overlay */
1365 #endif
1366
1367 /**********************************************************************/
bntostr(char * s,int dec,bn_t r)1368 char *bntostr(char *s, int dec, bn_t r)
1369 {
1370 return unsafe_bntostr(s, dec, copy_bn(bntmpcpy2, r));
1371 }
1372
1373 #if (_MSC_VER >= 700)
1374 #pragma code_seg ( ) /* back to normal segment */
1375 #endif
1376
1377 /**********************************************************************/
inv_bn(bn_t r,bn_t n)1378 bn_t inv_bn(bn_t r, bn_t n)
1379 {
1380 int sign;
1381
1382 sign = is_bn_neg(n);
1383 unsafe_inv_bn(r, n);
1384 if (sign)
1385 neg_a_bn(n);
1386 return r;
1387 }
1388
1389 /**********************************************************************/
div_bn(bn_t r,bn_t n1,bn_t n2)1390 bn_t div_bn(bn_t r, bn_t n1, bn_t n2)
1391 {
1392 copy_bn(bntmpcpy1, n1);
1393 copy_bn(bntmpcpy2, n2);
1394 return unsafe_div_bn(r, bntmpcpy1, bntmpcpy2);
1395 }
1396
1397 /**********************************************************************/
ln_bn(bn_t r,bn_t n)1398 bn_t ln_bn(bn_t r, bn_t n)
1399 {
1400 #if 0
1401 int sign;
1402
1403 sign = is_bn_neg(n);
1404 unsafe_ln_bn(r, n);
1405 if (sign)
1406 neg_a_bn(n);
1407 #endif
1408 copy_bn(bntmpcpy1, n); /* allows r and n to overlap memory */
1409 unsafe_ln_bn(r, bntmpcpy1);
1410 return r;
1411 }
1412
1413 /**********************************************************************/
sincos_bn(bn_t s,bn_t c,bn_t n)1414 bn_t sincos_bn(bn_t s, bn_t c, bn_t n)
1415 {
1416 return unsafe_sincos_bn(s, c, copy_bn(bntmpcpy1, n));
1417 }
1418
1419 /**********************************************************************/
atan_bn(bn_t r,bn_t n)1420 bn_t atan_bn(bn_t r, bn_t n)
1421 {
1422 int sign;
1423
1424 sign = is_bn_neg(n);
1425 unsafe_atan_bn(r, n);
1426 if (sign)
1427 neg_a_bn(n);
1428 return r;
1429 }
1430
1431 /**********************************************************************/
atan2_bn(bn_t r,bn_t ny,bn_t nx)1432 bn_t atan2_bn(bn_t r, bn_t ny, bn_t nx)
1433 {
1434 copy_bn(bntmpcpy1, ny);
1435 copy_bn(bntmpcpy2, nx);
1436 unsafe_atan2_bn(r, bntmpcpy1, bntmpcpy2);
1437 return r;
1438 }
1439
1440
1441 /**********************************************************************/
1442 /* Tim's miscellaneous stuff follows */
1443
1444 /**********************************************************************/
is_bn_zero(bn_t n)1445 int is_bn_zero(bn_t n)
1446 {
1447 return !is_bn_not_zero(n);
1448 }
1449