* Add this nice filesystem testing tool that I've recently
[dragonfly.git] / contrib / bc / lib / number.c
1 /* number.c: Implements arbitrary precision numbers. */
2 /*
3     Copyright (C) 1991, 1992, 1993, 1994, 1997, 2000 Free Software Foundation, Inc.
4
5     This program is free software; you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation; either version 2 of the License , or
8     (at your option) any later version.
9
10     This program is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13     GNU General Public License for more details.
14
15     You should have received a copy of the GNU General Public License
16     along with this program; see the file COPYING.  If not, write to:
17
18       The Free Software Foundation, Inc.
19       59 Temple Place, Suite 330
20       Boston, MA 02111-1307 USA.
21
22
23     You may contact the author by:
24        e-mail:  philnelson@acm.org
25       us-mail:  Philip A. Nelson
26                 Computer Science Department, 9062
27                 Western Washington University
28                 Bellingham, WA 98226-9062
29        
30 *************************************************************************/
31
32 #include <stdio.h>
33 #include <config.h>
34 #include <number.h>
35 #include <assert.h>
36 #include <stdlib.h>
37 #include <ctype.h>/* Prototypes needed for external utility routines. */
38
39 #define bc_rt_warn rt_warn
40 #define bc_rt_error rt_error
41 #define bc_out_of_memory out_of_memory
42
43 _PROTOTYPE(void rt_warn, (char *mesg ,...));
44 _PROTOTYPE(void rt_error, (char *mesg ,...));
45 _PROTOTYPE(void out_of_memory, (void));
46
47 /* Storage used for special numbers. */
48 bc_num _zero_;
49 bc_num _one_;
50 bc_num _two_;
51
52 static bc_num _bc_Free_list = NULL;
53
54 /* new_num allocates a number and sets fields to known values. */
55
56 bc_num
57 bc_new_num (length, scale)
58      int length, scale;
59 {
60   bc_num temp;
61
62   if (_bc_Free_list != NULL) {
63     temp = _bc_Free_list;
64     _bc_Free_list = temp->n_next;
65   } else {
66     temp = (bc_num) malloc (sizeof(bc_struct));
67     if (temp == NULL) bc_out_of_memory ();
68   }
69   temp->n_sign = PLUS;
70   temp->n_len = length;
71   temp->n_scale = scale;
72   temp->n_refs = 1;
73   temp->n_ptr = (char *) malloc (length+scale);
74   if (temp->n_ptr == NULL) bc_out_of_memory();
75   temp->n_value = temp->n_ptr;
76   memset (temp->n_ptr, 0, length+scale);
77   return temp;
78 }
79
80 /* "Frees" a bc_num NUM.  Actually decreases reference count and only
81    frees the storage if reference count is zero. */
82
83 void
84 bc_free_num (num)
85     bc_num *num;
86 {
87   if (*num == NULL) return;
88   (*num)->n_refs--;
89   if ((*num)->n_refs == 0) {
90     if ((*num)->n_ptr)
91       free ((*num)->n_ptr);
92     (*num)->n_next = _bc_Free_list;
93     _bc_Free_list = *num;
94   }
95   *num = NULL;
96 }
97
98
99 /* Intitialize the number package! */
100
101 void
102 bc_init_numbers ()
103 {
104   _zero_ = bc_new_num (1,0);
105   _one_  = bc_new_num (1,0);
106   _one_->n_value[0] = 1;
107   _two_  = bc_new_num (1,0);
108   _two_->n_value[0] = 2;
109 }
110
111
112 /* Make a copy of a number!  Just increments the reference count! */
113
114 bc_num
115 bc_copy_num (num)
116      bc_num num;
117 {
118   num->n_refs++;
119   return num;
120 }
121
122
123 /* Initialize a number NUM by making it a copy of zero. */
124
125 void
126 bc_init_num (num)
127      bc_num *num;
128 {
129   *num = bc_copy_num (_zero_);
130 }
131
132 /* For many things, we may have leading zeros in a number NUM.
133    _bc_rm_leading_zeros just moves the data "value" pointer to the
134    correct place and adjusts the length. */
135
136 static void
137 _bc_rm_leading_zeros (num)
138      bc_num num;
139 {
140   /* We can move n_value to point to the first non zero digit! */
141   while (*num->n_value == 0 && num->n_len > 1) {
142     num->n_value++;
143     num->n_len--;
144   }
145 }
146
147
148 /* Compare two bc numbers.  Return value is 0 if equal, -1 if N1 is less
149    than N2 and +1 if N1 is greater than N2.  If USE_SIGN is false, just
150    compare the magnitudes. */
151
152 static int
153 _bc_do_compare (n1, n2, use_sign, ignore_last)
154      bc_num n1, n2;
155      int use_sign;
156      int ignore_last;
157 {
158   char *n1ptr, *n2ptr;
159   int  count;
160
161   /* First, compare signs. */
162   if (use_sign && n1->n_sign != n2->n_sign)
163     {
164       if (n1->n_sign == PLUS)
165         return (1);     /* Positive N1 > Negative N2 */
166       else
167         return (-1);    /* Negative N1 < Positive N1 */
168     }
169
170   /* Now compare the magnitude. */
171   if (n1->n_len != n2->n_len)
172     {
173       if (n1->n_len > n2->n_len)
174         {
175           /* Magnitude of n1 > n2. */
176           if (!use_sign || n1->n_sign == PLUS)
177             return (1);
178           else
179             return (-1);
180         }
181       else
182         {
183           /* Magnitude of n1 < n2. */
184           if (!use_sign || n1->n_sign == PLUS)
185             return (-1);
186           else
187             return (1);
188         }
189     }
190
191   /* If we get here, they have the same number of integer digits.
192      check the integer part and the equal length part of the fraction. */
193   count = n1->n_len + MIN (n1->n_scale, n2->n_scale);
194   n1ptr = n1->n_value;
195   n2ptr = n2->n_value;
196
197   while ((count > 0) && (*n1ptr == *n2ptr))
198     {
199       n1ptr++;
200       n2ptr++;
201       count--;
202     }
203   if (ignore_last && count == 1 && n1->n_scale == n2->n_scale)
204     return (0);
205   if (count != 0)
206     {
207       if (*n1ptr > *n2ptr)
208         {
209           /* Magnitude of n1 > n2. */
210           if (!use_sign || n1->n_sign == PLUS)
211             return (1);
212           else
213             return (-1);
214         }
215       else
216         {
217           /* Magnitude of n1 < n2. */
218           if (!use_sign || n1->n_sign == PLUS)
219             return (-1);
220           else
221             return (1);
222         }
223     }
224
225   /* They are equal up to the last part of the equal part of the fraction. */
226   if (n1->n_scale != n2->n_scale)
227     {
228       if (n1->n_scale > n2->n_scale)
229         {
230           for (count = n1->n_scale-n2->n_scale; count>0; count--)
231             if (*n1ptr++ != 0)
232               {
233                 /* Magnitude of n1 > n2. */
234                 if (!use_sign || n1->n_sign == PLUS)
235                   return (1);
236                 else
237                   return (-1);
238               }
239         }
240       else
241         {
242           for (count = n2->n_scale-n1->n_scale; count>0; count--)
243             if (*n2ptr++ != 0)
244               {
245                 /* Magnitude of n1 < n2. */
246                 if (!use_sign || n1->n_sign == PLUS)
247                   return (-1);
248                 else
249                   return (1);
250               }
251         }
252     }
253
254   /* They must be equal! */
255   return (0);
256 }
257
258
259 /* This is the "user callable" routine to compare numbers N1 and N2. */
260
261 int
262 bc_compare (n1, n2)
263      bc_num n1, n2;
264 {
265   return _bc_do_compare (n1, n2, TRUE, FALSE);
266 }
267
268 /* In some places we need to check if the number is negative. */
269
270 char
271 bc_is_neg (num)
272      bc_num num;
273 {
274   return num->n_sign == MINUS;
275 }
276
277 /* In some places we need to check if the number NUM is zero. */
278
279 char
280 bc_is_zero (num)
281      bc_num num;
282 {
283   int  count;
284   char *nptr;
285
286   /* Quick check. */
287   if (num == _zero_) return TRUE;
288
289   /* Initialize */
290   count = num->n_len + num->n_scale;
291   nptr = num->n_value;
292
293   /* The check */
294   while ((count > 0) && (*nptr++ == 0)) count--;
295
296   if (count != 0)
297     return FALSE;
298   else
299     return TRUE;
300 }
301
302 /* In some places we need to check if the number NUM is almost zero.
303    Specifically, all but the last digit is 0 and the last digit is 1.
304    Last digit is defined by scale. */
305
306 char
307 bc_is_near_zero (num, scale)
308      bc_num num;
309      int scale;
310 {
311   int  count;
312   char *nptr;
313
314   /* Error checking */
315   if (scale > num->n_scale)
316     scale = num->n_scale;
317
318   /* Initialize */
319   count = num->n_len + scale;
320   nptr = num->n_value;
321
322   /* The check */
323   while ((count > 0) && (*nptr++ == 0)) count--;
324
325   if (count != 0 && (count != 1 || *--nptr != 1))
326     return FALSE;
327   else
328     return TRUE;
329 }
330
331
332 /* Perform addition: N1 is added to N2 and the value is
333    returned.  The signs of N1 and N2 are ignored.
334    SCALE_MIN is to set the minimum scale of the result. */
335
336 static bc_num
337 _bc_do_add (n1, n2, scale_min)
338      bc_num n1, n2;
339      int scale_min;
340 {
341   bc_num sum;
342   int sum_scale, sum_digits;
343   char *n1ptr, *n2ptr, *sumptr;
344   int carry, n1bytes, n2bytes;
345   int count;
346
347   /* Prepare sum. */
348   sum_scale = MAX (n1->n_scale, n2->n_scale);
349   sum_digits = MAX (n1->n_len, n2->n_len) + 1;
350   sum = bc_new_num (sum_digits, MAX(sum_scale, scale_min));
351
352   /* Zero extra digits made by scale_min. */
353   if (scale_min > sum_scale)
354     {
355       sumptr = (char *) (sum->n_value + sum_scale + sum_digits);
356       for (count = scale_min - sum_scale; count > 0; count--)
357         *sumptr++ = 0;
358     }
359
360   /* Start with the fraction part.  Initialize the pointers. */
361   n1bytes = n1->n_scale;
362   n2bytes = n2->n_scale;
363   n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1);
364   n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1);
365   sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1);
366
367   /* Add the fraction part.  First copy the longer fraction.*/
368   if (n1bytes != n2bytes)
369     {
370       if (n1bytes > n2bytes)
371         while (n1bytes>n2bytes)
372           { *sumptr-- = *n1ptr--; n1bytes--;}
373       else
374         while (n2bytes>n1bytes)
375           { *sumptr-- = *n2ptr--; n2bytes--;}
376     }
377
378   /* Now add the remaining fraction part and equal size integer parts. */
379   n1bytes += n1->n_len;
380   n2bytes += n2->n_len;
381   carry = 0;
382   while ((n1bytes > 0) && (n2bytes > 0))
383     {
384       *sumptr = *n1ptr-- + *n2ptr-- + carry;
385       if (*sumptr > (BASE-1))
386         {
387            carry = 1;
388            *sumptr -= BASE;
389         }
390       else
391         carry = 0;
392       sumptr--;
393       n1bytes--;
394       n2bytes--;
395     }
396
397   /* Now add carry the longer integer part. */
398   if (n1bytes == 0)
399     { n1bytes = n2bytes; n1ptr = n2ptr; }
400   while (n1bytes-- > 0)
401     {
402       *sumptr = *n1ptr-- + carry;
403       if (*sumptr > (BASE-1))
404         {
405            carry = 1;
406            *sumptr -= BASE;
407          }
408       else
409         carry = 0;
410       sumptr--;
411     }
412
413   /* Set final carry. */
414   if (carry == 1)
415     *sumptr += 1;
416
417   /* Adjust sum and return. */
418   _bc_rm_leading_zeros (sum);
419   return sum;
420 }
421
422
423 /* Perform subtraction: N2 is subtracted from N1 and the value is
424    returned.  The signs of N1 and N2 are ignored.  Also, N1 is
425    assumed to be larger than N2.  SCALE_MIN is the minimum scale
426    of the result. */
427
428 static bc_num
429 _bc_do_sub (n1, n2, scale_min)
430      bc_num n1, n2;
431      int scale_min;
432 {
433   bc_num diff;
434   int diff_scale, diff_len;
435   int min_scale, min_len;
436   char *n1ptr, *n2ptr, *diffptr;
437   int borrow, count, val;
438
439   /* Allocate temporary storage. */
440   diff_len = MAX (n1->n_len, n2->n_len);
441   diff_scale = MAX (n1->n_scale, n2->n_scale);
442   min_len = MIN  (n1->n_len, n2->n_len);
443   min_scale = MIN (n1->n_scale, n2->n_scale);
444   diff = bc_new_num (diff_len, MAX(diff_scale, scale_min));
445
446   /* Zero extra digits made by scale_min. */
447   if (scale_min > diff_scale)
448     {
449       diffptr = (char *) (diff->n_value + diff_len + diff_scale);
450       for (count = scale_min - diff_scale; count > 0; count--)
451         *diffptr++ = 0;
452     }
453
454   /* Initialize the subtract. */
455   n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1);
456   n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1);
457   diffptr = (char *) (diff->n_value + diff_len + diff_scale -1);
458
459   /* Subtract the numbers. */
460   borrow = 0;
461
462   /* Take care of the longer scaled number. */
463   if (n1->n_scale != min_scale)
464     {
465       /* n1 has the longer scale */
466       for (count = n1->n_scale - min_scale; count > 0; count--)
467         *diffptr-- = *n1ptr--;
468     }
469   else
470     {
471       /* n2 has the longer scale */
472       for (count = n2->n_scale - min_scale; count > 0; count--)
473         {
474           val = - *n2ptr-- - borrow;
475           if (val < 0)
476             {
477               val += BASE;
478               borrow = 1;
479             }
480           else
481             borrow = 0;
482           *diffptr-- = val;
483         }
484     }
485
486   /* Now do the equal length scale and integer parts. */
487
488   for (count = 0; count < min_len + min_scale; count++)
489     {
490       val = *n1ptr-- - *n2ptr-- - borrow;
491       if (val < 0)
492         {
493           val += BASE;
494           borrow = 1;
495         }
496       else
497         borrow = 0;
498       *diffptr-- = val;
499     }
500
501   /* If n1 has more digits then n2, we now do that subtract. */
502   if (diff_len != min_len)
503     {
504       for (count = diff_len - min_len; count > 0; count--)
505         {
506           val = *n1ptr-- - borrow;
507           if (val < 0)
508             {
509               val += BASE;
510               borrow = 1;
511             }
512           else
513             borrow = 0;
514           *diffptr-- = val;
515         }
516     }
517
518   /* Clean up and return. */
519   _bc_rm_leading_zeros (diff);
520   return diff;
521 }
522
523
524 /* Here is the full subtract routine that takes care of negative numbers.
525    N2 is subtracted from N1 and the result placed in RESULT.  SCALE_MIN
526    is the minimum scale for the result. */
527
528 void
529 bc_sub (n1, n2, result, scale_min)
530      bc_num n1, n2, *result;
531      int scale_min;
532 {
533   bc_num diff = NULL;
534   int cmp_res;
535   int res_scale;
536
537   if (n1->n_sign != n2->n_sign)
538     {
539       diff = _bc_do_add (n1, n2, scale_min);
540       diff->n_sign = n1->n_sign;
541     }
542   else
543     {
544       /* subtraction must be done. */
545       /* Compare magnitudes. */
546       cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);
547       switch (cmp_res)
548         {
549         case -1:
550           /* n1 is less than n2, subtract n1 from n2. */
551           diff = _bc_do_sub (n2, n1, scale_min);
552           diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS);
553           break;
554         case  0:
555           /* They are equal! return zero! */
556           res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
557           diff = bc_new_num (1, res_scale);
558           memset (diff->n_value, 0, res_scale+1);
559           break;
560         case  1:
561           /* n2 is less than n1, subtract n2 from n1. */
562           diff = _bc_do_sub (n1, n2, scale_min);
563           diff->n_sign = n1->n_sign;
564           break;
565         }
566     }
567
568   /* Clean up and return. */
569   bc_free_num (result);
570   *result = diff;
571 }
572
573
574 /* Here is the full add routine that takes care of negative numbers.
575    N1 is added to N2 and the result placed into RESULT.  SCALE_MIN
576    is the minimum scale for the result. */
577
578 void
579 bc_add (n1, n2, result, scale_min)
580      bc_num n1, n2, *result;
581      int scale_min;
582 {
583   bc_num sum = NULL;
584   int cmp_res;
585   int res_scale;
586
587   if (n1->n_sign == n2->n_sign)
588     {
589       sum = _bc_do_add (n1, n2, scale_min);
590       sum->n_sign = n1->n_sign;
591     }
592   else
593     {
594       /* subtraction must be done. */
595       cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);  /* Compare magnitudes. */
596       switch (cmp_res)
597         {
598         case -1:
599           /* n1 is less than n2, subtract n1 from n2. */
600           sum = _bc_do_sub (n2, n1, scale_min);
601           sum->n_sign = n2->n_sign;
602           break;
603         case  0:
604           /* They are equal! return zero with the correct scale! */
605           res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
606           sum = bc_new_num (1, res_scale);
607           memset (sum->n_value, 0, res_scale+1);
608           break;
609         case  1:
610           /* n2 is less than n1, subtract n2 from n1. */
611           sum = _bc_do_sub (n1, n2, scale_min);
612           sum->n_sign = n1->n_sign;
613         }
614     }
615
616   /* Clean up and return. */
617   bc_free_num (result);
618   *result = sum;
619 }
620
621 /* Recursive vs non-recursive multiply crossover ranges. */
622 #if defined(MULDIGITS)
623 #include "muldigits.h"
624 #else
625 #define MUL_BASE_DIGITS 80
626 #endif
627
628 int mul_base_digits = MUL_BASE_DIGITS;
629 #define MUL_SMALL_DIGITS mul_base_digits/4
630
631 /* Multiply utility routines */
632
633 static bc_num
634 new_sub_num (length, scale, value)
635      int length, scale;
636      char *value;
637 {
638   bc_num temp;
639
640   if (_bc_Free_list != NULL) {
641     temp = _bc_Free_list;
642     _bc_Free_list = temp->n_next;
643   } else {
644     temp = (bc_num) malloc (sizeof(bc_struct));
645     if (temp == NULL) bc_out_of_memory ();
646   }
647   temp->n_sign = PLUS;
648   temp->n_len = length;
649   temp->n_scale = scale;
650   temp->n_refs = 1;
651   temp->n_ptr = NULL;
652   temp->n_value = value;
653   return temp;
654 }
655
656 static void
657 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
658               int full_scale)
659 {
660   char *n1ptr, *n2ptr, *pvptr;
661   char *n1end, *n2end;          /* To the end of n1 and n2. */
662   int indx, sum, prodlen;
663
664   prodlen = n1len+n2len+1;
665
666   *prod = bc_new_num (prodlen, 0);
667
668   n1end = (char *) (n1->n_value + n1len - 1);
669   n2end = (char *) (n2->n_value + n2len - 1);
670   pvptr = (char *) ((*prod)->n_value + prodlen - 1);
671   sum = 0;
672
673   /* Here is the loop... */
674   for (indx = 0; indx < prodlen-1; indx++)
675     {
676       n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
677       n2ptr = (char *) (n2end - MIN(indx, n2len-1));
678       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
679         sum += *n1ptr-- * *n2ptr++;
680       *pvptr-- = sum % BASE;
681       sum = sum / BASE;
682     }
683   *pvptr = sum;
684 }
685
686
687 /* A special adder/subtractor for the recursive divide and conquer
688    multiply algorithm.  Note: if sub is called, accum must
689    be larger that what is being subtracted.  Also, accum and val
690    must have n_scale = 0.  (e.g. they must look like integers. *) */
691 static void
692 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
693 {
694   signed char *accp, *valp;
695   int  count, carry;
696
697   count = val->n_len;
698   if (val->n_value[0] == 0)
699     count--;
700   assert (accum->n_len+accum->n_scale >= shift+count);
701   
702   /* Set up pointers and others */
703   accp = (signed char *)(accum->n_value +
704                          accum->n_len + accum->n_scale - shift - 1);
705   valp = (signed char *)(val->n_value + val->n_len - 1);
706   carry = 0;
707
708   if (sub) {
709     /* Subtraction, carry is really borrow. */
710     while (count--) {
711       *accp -= *valp-- + carry;
712       if (*accp < 0) {
713         carry = 1;
714         *accp-- += BASE;
715       } else {
716         carry = 0;
717         accp--;
718       }
719     }
720     while (carry) {
721       *accp -= carry;
722       if (*accp < 0)
723         *accp-- += BASE;
724       else
725         carry = 0;
726     }
727   } else {
728     /* Addition */
729     while (count--) {
730       *accp += *valp-- + carry;
731       if (*accp > (BASE-1)) {
732         carry = 1;
733         *accp-- -= BASE;
734       } else {
735         carry = 0;
736         accp--;
737       }
738     }
739     while (carry) {
740       *accp += carry;
741       if (*accp > (BASE-1))
742         *accp-- -= BASE;
743       else
744         carry = 0;
745     }
746   }
747 }
748
749 /* Recursive divide and conquer multiply algorithm.  
750    Based on 
751    Let u = u0 + u1*(b^n)
752    Let v = v0 + v1*(b^n)
753    Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
754
755    B is the base of storage, number of digits in u1,u0 close to equal.
756 */
757 static void
758 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
759              int full_scale)
760
761   bc_num u0, u1, v0, v1;
762   int u0len, v0len;
763   bc_num m1, m2, m3, d1, d2;
764   int n, prodlen, m1zero;
765   int d1len, d2len;
766
767   /* Base case? */
768   if ((ulen+vlen) < mul_base_digits
769       || ulen < MUL_SMALL_DIGITS
770       || vlen < MUL_SMALL_DIGITS ) {
771     _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
772     return;
773   }
774
775   /* Calculate n -- the u and v split point in digits. */
776   n = (MAX(ulen, vlen)+1) / 2;
777
778   /* Split u and v. */
779   if (ulen < n) {
780     u1 = bc_copy_num (_zero_);
781     u0 = new_sub_num (ulen,0, u->n_value);
782   } else {
783     u1 = new_sub_num (ulen-n, 0, u->n_value);
784     u0 = new_sub_num (n, 0, u->n_value+ulen-n);
785   }
786   if (vlen < n) {
787     v1 = bc_copy_num (_zero_);
788     v0 = new_sub_num (vlen,0, v->n_value);
789   } else {
790     v1 = new_sub_num (vlen-n, 0, v->n_value);
791     v0 = new_sub_num (n, 0, v->n_value+vlen-n);
792     }
793   _bc_rm_leading_zeros (u1);
794   _bc_rm_leading_zeros (u0);
795   u0len = u0->n_len;
796   _bc_rm_leading_zeros (v1);
797   _bc_rm_leading_zeros (v0);
798   v0len = v0->n_len;
799
800   m1zero = bc_is_zero(u1) || bc_is_zero(v1);
801
802   /* Calculate sub results ... */
803
804   bc_init_num(&d1);
805   bc_init_num(&d2);
806   bc_sub (u1, u0, &d1, 0);
807   d1len = d1->n_len;
808   bc_sub (v0, v1, &d2, 0);
809   d2len = d2->n_len;
810
811
812   /* Do recursive multiplies and shifted adds. */
813   if (m1zero)
814     m1 = bc_copy_num (_zero_);
815   else
816     _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0);
817
818   if (bc_is_zero(d1) || bc_is_zero(d2))
819     m2 = bc_copy_num (_zero_);
820   else
821     _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0);
822
823   if (bc_is_zero(u0) || bc_is_zero(v0))
824     m3 = bc_copy_num (_zero_);
825   else
826     _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0);
827
828   /* Initialize product */
829   prodlen = ulen+vlen+1;
830   *prod = bc_new_num(prodlen, 0);
831
832   if (!m1zero) {
833     _bc_shift_addsub (*prod, m1, 2*n, 0);
834     _bc_shift_addsub (*prod, m1, n, 0);
835   }
836   _bc_shift_addsub (*prod, m3, n, 0);
837   _bc_shift_addsub (*prod, m3, 0, 0);
838   _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
839
840   /* Now clean up! */
841   bc_free_num (&u1);
842   bc_free_num (&u0);
843   bc_free_num (&v1);
844   bc_free_num (&m1);
845   bc_free_num (&v0);
846   bc_free_num (&m2);
847   bc_free_num (&m3);
848   bc_free_num (&d1);
849   bc_free_num (&d2);
850 }
851
852 /* The multiply routine.  N2 times N1 is put int PROD with the scale of
853    the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
854    */
855
856 void
857 bc_multiply (n1, n2, prod, scale)
858      bc_num n1, n2, *prod;
859      int scale;
860 {
861   bc_num pval; 
862   int len1, len2;
863   int full_scale, prod_scale;
864
865   /* Initialize things. */
866   len1 = n1->n_len + n1->n_scale;
867   len2 = n2->n_len + n2->n_scale;
868   full_scale = n1->n_scale + n2->n_scale;
869   prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
870
871   /* Do the multiply */
872   _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale);
873
874   /* Assign to prod and clean up the number. */
875   pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
876   pval->n_value = pval->n_ptr;
877   pval->n_len = len2 + len1 + 1 - full_scale;
878   pval->n_scale = prod_scale;
879   _bc_rm_leading_zeros (pval);
880   if (bc_is_zero (pval))
881     pval->n_sign = PLUS;
882   bc_free_num (prod);
883   *prod = pval;
884 }
885
886 /* Some utility routines for the divide:  First a one digit multiply.
887    NUM (with SIZE digits) is multiplied by DIGIT and the result is
888    placed into RESULT.  It is written so that NUM and RESULT can be
889    the same pointers.  */
890
891 static void
892 _one_mult (num, size, digit, result)
893      unsigned char *num;
894      int size, digit;
895      unsigned char *result;
896 {
897   int carry, value;
898   unsigned char *nptr, *rptr;
899
900   if (digit == 0)
901     memset (result, 0, size);
902   else
903     {
904       if (digit == 1)
905         memcpy (result, num, size);
906       else
907         {
908           /* Initialize */
909           nptr = (unsigned char *) (num+size-1);
910           rptr = (unsigned char *) (result+size-1);
911           carry = 0;
912
913           while (size-- > 0)
914             {
915               value = *nptr-- * digit + carry;
916               *rptr-- = value % BASE;
917               carry = value / BASE;
918             }
919
920           if (carry != 0) *rptr = carry;
921         }
922     }
923 }
924
925
926 /* The full division routine. This computes N1 / N2.  It returns
927    0 if the division is ok and the result is in QUOT.  The number of
928    digits after the decimal point is SCALE. It returns -1 if division
929    by zero is tried.  The algorithm is found in Knuth Vol 2. p237. */
930
931 int
932 bc_divide (n1, n2, quot, scale)
933      bc_num n1, n2, *quot;
934      int scale;
935 {
936   bc_num qval;
937   unsigned char *num1, *num2;
938   unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
939   int  scale1, val;
940   unsigned int  len1, len2, scale2, qdigits, extra, count;
941   unsigned int  qdig, qguess, borrow, carry;
942   unsigned char *mval;
943   char zero;
944   unsigned int  norm;
945
946   /* Test for divide by zero. */
947   if (bc_is_zero (n2)) return -1;
948
949   /* Test for divide by 1.  If it is we must truncate. */
950   if (n2->n_scale == 0)
951     {
952       if (n2->n_len == 1 && *n2->n_value == 1)
953         {
954           qval = bc_new_num (n1->n_len, scale);
955           qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
956           memset (&qval->n_value[n1->n_len],0,scale);
957           memcpy (qval->n_value, n1->n_value,
958                   n1->n_len + MIN(n1->n_scale,scale));
959           bc_free_num (quot);
960           *quot = qval;
961         }
962     }
963
964   /* Set up the divide.  Move the decimal point on n1 by n2's scale.
965      Remember, zeros on the end of num2 are wasted effort for dividing. */
966   scale2 = n2->n_scale;
967   n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
968   while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
969
970   len1 = n1->n_len + scale2;
971   scale1 = n1->n_scale - scale2;
972   if (scale1 < scale)
973     extra = scale - scale1;
974   else
975     extra = 0;
976   num1 = (unsigned char *) malloc (n1->n_len+n1->n_scale+extra+2);
977   if (num1 == NULL) bc_out_of_memory();
978   memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
979   memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
980
981   len2 = n2->n_len + scale2;
982   num2 = (unsigned char *) malloc (len2+1);
983   if (num2 == NULL) bc_out_of_memory();
984   memcpy (num2, n2->n_value, len2);
985   *(num2+len2) = 0;
986   n2ptr = num2;
987   while (*n2ptr == 0)
988     {
989       n2ptr++;
990       len2--;
991     }
992
993   /* Calculate the number of quotient digits. */
994   if (len2 > len1+scale)
995     {
996       qdigits = scale+1;
997       zero = TRUE;
998     }
999   else
1000     {
1001       zero = FALSE;
1002       if (len2>len1)
1003         qdigits = scale+1;      /* One for the zero integer part. */
1004       else
1005         qdigits = len1-len2+scale+1;
1006     }
1007
1008   /* Allocate and zero the storage for the quotient. */
1009   qval = bc_new_num (qdigits-scale,scale);
1010   memset (qval->n_value, 0, qdigits);
1011
1012   /* Allocate storage for the temporary storage mval. */
1013   mval = (unsigned char *) malloc (len2+1);
1014   if (mval == NULL) bc_out_of_memory ();
1015
1016   /* Now for the full divide algorithm. */
1017   if (!zero)
1018     {
1019       /* Normalize */
1020       norm =  10 / ((int)*n2ptr + 1);
1021       if (norm != 1)
1022         {
1023           _one_mult (num1, len1+scale1+extra+1, norm, num1);
1024           _one_mult (n2ptr, len2, norm, n2ptr);
1025         }
1026
1027       /* Initialize divide loop. */
1028       qdig = 0;
1029       if (len2 > len1)
1030         qptr = (unsigned char *) qval->n_value+len2-len1;
1031       else
1032         qptr = (unsigned char *) qval->n_value;
1033
1034       /* Loop */
1035       while (qdig <= len1+scale-len2)
1036         {
1037           /* Calculate the quotient digit guess. */
1038           if (*n2ptr == num1[qdig])
1039             qguess = 9;
1040           else
1041             qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
1042
1043           /* Test qguess. */
1044           if (n2ptr[1]*qguess >
1045               (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
1046                + num1[qdig+2])
1047             {
1048               qguess--;
1049               /* And again. */
1050               if (n2ptr[1]*qguess >
1051                   (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
1052                   + num1[qdig+2])
1053                 qguess--;
1054             }
1055
1056           /* Multiply and subtract. */
1057           borrow = 0;
1058           if (qguess != 0)
1059             {
1060               *mval = 0;
1061               _one_mult (n2ptr, len2, qguess, mval+1);
1062               ptr1 = (unsigned char *) num1+qdig+len2;
1063               ptr2 = (unsigned char *) mval+len2;
1064               for (count = 0; count < len2+1; count++)
1065                 {
1066                   val = (int) *ptr1 - (int) *ptr2-- - borrow;
1067                   if (val < 0)
1068                     {
1069                       val += 10;
1070                       borrow = 1;
1071                     }
1072                   else
1073                     borrow = 0;
1074                   *ptr1-- = val;
1075                 }
1076             }
1077
1078           /* Test for negative result. */
1079           if (borrow == 1)
1080             {
1081               qguess--;
1082               ptr1 = (unsigned char *) num1+qdig+len2;
1083               ptr2 = (unsigned char *) n2ptr+len2-1;
1084               carry = 0;
1085               for (count = 0; count < len2; count++)
1086                 {
1087                   val = (int) *ptr1 + (int) *ptr2-- + carry;
1088                   if (val > 9)
1089                     {
1090                       val -= 10;
1091                       carry = 1;
1092                     }
1093                   else
1094                     carry = 0;
1095                   *ptr1-- = val;
1096                 }
1097               if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
1098             }
1099
1100           /* We now know the quotient digit. */
1101           *qptr++ =  qguess;
1102           qdig++;
1103         }
1104     }
1105
1106   /* Clean up and return the number. */
1107   qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
1108   if (bc_is_zero (qval)) qval->n_sign = PLUS;
1109   _bc_rm_leading_zeros (qval);
1110   bc_free_num (quot);
1111   *quot = qval;
1112
1113   /* Clean up temporary storage. */
1114   free (mval);
1115   free (num1);
1116   free (num2);
1117
1118   return 0;     /* Everything is OK. */
1119 }
1120
1121
1122 /* Division *and* modulo for numbers.  This computes both NUM1 / NUM2 and
1123    NUM1 % NUM2  and puts the results in QUOT and REM, except that if QUOT
1124    is NULL then that store will be omitted.
1125  */
1126
1127 int
1128 bc_divmod (num1, num2, quot, rem, scale)
1129      bc_num num1, num2, *quot, *rem;
1130      int scale;
1131 {
1132   bc_num quotient = NULL;
1133   bc_num temp;
1134   int rscale;
1135
1136   /* Check for correct numbers. */
1137   if (bc_is_zero (num2)) return -1;
1138
1139   /* Calculate final scale. */
1140   rscale = MAX (num1->n_scale, num2->n_scale+scale);
1141   bc_init_num(&temp);
1142
1143   /* Calculate it. */
1144   bc_divide (num1, num2, &temp, scale);
1145   if (quot)
1146     quotient = bc_copy_num (temp);
1147   bc_multiply (temp, num2, &temp, rscale);
1148   bc_sub (num1, temp, rem, rscale);
1149   bc_free_num (&temp);
1150
1151   if (quot)
1152     {
1153       bc_free_num (quot);
1154       *quot = quotient;
1155     }
1156
1157   return 0;     /* Everything is OK. */
1158 }
1159
1160
1161 /* Modulo for numbers.  This computes NUM1 % NUM2  and puts the
1162    result in RESULT.   */
1163
1164 int
1165 bc_modulo (num1, num2, result, scale)
1166      bc_num num1, num2, *result;
1167      int scale;
1168 {
1169   return bc_divmod (num1, num2, NULL, result, scale);
1170 }
1171
1172 /* Raise BASE to the EXPO power, reduced modulo MOD.  The result is
1173    placed in RESULT.  If a EXPO is not an integer,
1174    only the integer part is used.  */
1175
1176 int
1177 bc_raisemod (base, expo, mod, result, scale)
1178      bc_num base, expo, mod, *result;
1179      int scale;
1180 {
1181   bc_num power, exponent, parity, temp;
1182   int rscale;
1183
1184   /* Check for correct numbers. */
1185   if (bc_is_zero(mod)) return -1;
1186   if (bc_is_neg(expo)) return -1;
1187
1188   /* Set initial values.  */
1189   power = bc_copy_num (base);
1190   exponent = bc_copy_num (expo);
1191   temp = bc_copy_num (_one_);
1192   bc_init_num(&parity);
1193
1194   /* Check the base for scale digits. */
1195   if (base->n_scale != 0)
1196       bc_rt_warn ("non-zero scale in base");
1197
1198   /* Check the exponent for scale digits. */
1199   if (exponent->n_scale != 0)
1200     {
1201       bc_rt_warn ("non-zero scale in exponent");
1202       bc_divide (exponent, _one_, &exponent, 0); /*truncate */
1203     }
1204
1205   /* Check the modulus for scale digits. */
1206   if (mod->n_scale != 0)
1207       bc_rt_warn ("non-zero scale in modulus");
1208
1209   /* Do the calculation. */
1210   rscale = MAX(scale, base->n_scale);
1211   while ( !bc_is_zero(exponent) )
1212     {
1213       (void) bc_divmod (exponent, _two_, &exponent, &parity, 0);
1214       if ( !bc_is_zero(parity) )
1215         {
1216           bc_multiply (temp, power, &temp, rscale);
1217           (void) bc_modulo (temp, mod, &temp, scale);
1218         }
1219
1220       bc_multiply (power, power, &power, rscale);
1221       (void) bc_modulo (power, mod, &power, scale);
1222     }
1223
1224   /* Assign the value. */
1225   bc_free_num (&power);
1226   bc_free_num (&exponent);
1227   bc_free_num (result);
1228   *result = temp;
1229   return 0;     /* Everything is OK. */
1230 }
1231
1232 /* Raise NUM1 to the NUM2 power.  The result is placed in RESULT.
1233    Maximum exponent is LONG_MAX.  If a NUM2 is not an integer,
1234    only the integer part is used.  */
1235
1236 void
1237 bc_raise (num1, num2, result, scale)
1238      bc_num num1, num2, *result;
1239      int scale;
1240 {
1241    bc_num temp, power;
1242    long exponent;
1243    int rscale;
1244    int pwrscale;
1245    int calcscale;
1246    char neg;
1247
1248    /* Check the exponent for scale digits and convert to a long. */
1249    if (num2->n_scale != 0)
1250      bc_rt_warn ("non-zero scale in exponent");
1251    exponent = bc_num2long (num2);
1252    if (exponent == 0 && (num2->n_len > 1 || num2->n_value[0] != 0))
1253        bc_rt_error ("exponent too large in raise");
1254
1255    /* Special case if exponent is a zero. */
1256    if (exponent == 0)
1257      {
1258        bc_free_num (result);
1259        *result = bc_copy_num (_one_);
1260        return;
1261      }
1262
1263    /* Other initializations. */
1264    if (exponent < 0)
1265      {
1266        neg = TRUE;
1267        exponent = -exponent;
1268        rscale = scale;
1269      }
1270    else
1271      {
1272        neg = FALSE;
1273        rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale));
1274      }
1275
1276    /* Set initial value of temp.  */
1277    power = bc_copy_num (num1);
1278    pwrscale = num1->n_scale;
1279    while ((exponent & 1) == 0)
1280      {
1281        pwrscale = 2*pwrscale;
1282        bc_multiply (power, power, &power, pwrscale);
1283        exponent = exponent >> 1;
1284      }
1285    temp = bc_copy_num (power);
1286    calcscale = pwrscale;
1287    exponent = exponent >> 1;
1288
1289    /* Do the calculation. */
1290    while (exponent > 0)
1291      {
1292        pwrscale = 2*pwrscale;
1293        bc_multiply (power, power, &power, pwrscale);
1294        if ((exponent & 1) == 1) {
1295          calcscale = pwrscale + calcscale;
1296          bc_multiply (temp, power, &temp, calcscale);
1297        }
1298        exponent = exponent >> 1;
1299      }
1300
1301    /* Assign the value. */
1302    if (neg)
1303      {
1304        bc_divide (_one_, temp, result, rscale);
1305        bc_free_num (&temp);
1306      }
1307    else
1308      {
1309        bc_free_num (result);
1310        *result = temp;
1311        if ((*result)->n_scale > rscale)
1312          (*result)->n_scale = rscale;
1313      }
1314    bc_free_num (&power);
1315 }
1316
1317 /* Take the square root NUM and return it in NUM with SCALE digits
1318    after the decimal place. */
1319
1320 int
1321 bc_sqrt (num, scale)
1322      bc_num *num;
1323      int scale;
1324 {
1325   int rscale, cmp_res, done;
1326   int cscale;
1327   bc_num guess, guess1, point5, diff;
1328
1329   /* Initial checks. */
1330   cmp_res = bc_compare (*num, _zero_);
1331   if (cmp_res < 0)
1332     return 0;           /* error */
1333   else
1334     {
1335       if (cmp_res == 0)
1336         {
1337           bc_free_num (num);
1338           *num = bc_copy_num (_zero_);
1339           return 1;
1340         }
1341     }
1342   cmp_res = bc_compare (*num, _one_);
1343   if (cmp_res == 0)
1344     {
1345       bc_free_num (num);
1346       *num = bc_copy_num (_one_);
1347       return 1;
1348     }
1349
1350   /* Initialize the variables. */
1351   rscale = MAX (scale, (*num)->n_scale);
1352   bc_init_num(&guess);
1353   bc_init_num(&guess1);
1354   bc_init_num(&diff);
1355   point5 = bc_new_num (1,1);
1356   point5->n_value[1] = 5;
1357
1358
1359   /* Calculate the initial guess. */
1360   if (cmp_res < 0)
1361     {
1362       /* The number is between 0 and 1.  Guess should start at 1. */
1363       guess = bc_copy_num (_one_);
1364       cscale = (*num)->n_scale;
1365     }
1366   else
1367     {
1368       /* The number is greater than 1.  Guess should start at 10^(exp/2). */
1369       bc_int2num (&guess,10);
1370
1371       bc_int2num (&guess1,(*num)->n_len);
1372       bc_multiply (guess1, point5, &guess1, 0);
1373       guess1->n_scale = 0;
1374       bc_raise (guess, guess1, &guess, 0);
1375       bc_free_num (&guess1);
1376       cscale = 3;
1377     }
1378
1379   /* Find the square root using Newton's algorithm. */
1380   done = FALSE;
1381   while (!done)
1382     {
1383       bc_free_num (&guess1);
1384       guess1 = bc_copy_num (guess);
1385       bc_divide (*num, guess, &guess, cscale);
1386       bc_add (guess, guess1, &guess, 0);
1387       bc_multiply (guess, point5, &guess, cscale);
1388       bc_sub (guess, guess1, &diff, cscale+1);
1389       if (bc_is_near_zero (diff, cscale))
1390         {
1391           if (cscale < rscale+1)
1392             cscale = MIN (cscale*3, rscale+1);
1393           else
1394             done = TRUE;
1395         }
1396     }
1397
1398   /* Assign the number and clean up. */
1399   bc_free_num (num);
1400   bc_divide (guess,_one_,num,rscale);
1401   bc_free_num (&guess);
1402   bc_free_num (&guess1);
1403   bc_free_num (&point5);
1404   bc_free_num (&diff);
1405   return 1;
1406 }
1407
1408
1409 /* The following routines provide output for bcd numbers package
1410    using the rules of POSIX bc for output. */
1411
1412 /* This structure is used for saving digits in the conversion process. */
1413 typedef struct stk_rec {
1414         long  digit;
1415         struct stk_rec *next;
1416 } stk_rec;
1417
1418 /* The reference string for digits. */
1419 static char ref_str[] = "0123456789ABCDEF";
1420
1421
1422 /* A special output routine for "multi-character digits."  Exactly
1423    SIZE characters must be output for the value VAL.  If SPACE is
1424    non-zero, we must output one space before the number.  OUT_CHAR
1425    is the actual routine for writing the characters. */
1426
1427 void
1428 bc_out_long (val, size, space, out_char)
1429      long val;
1430      int size, space;
1431 #ifdef __STDC__
1432      void (*out_char)(int);
1433 #else
1434      void (*out_char)();
1435 #endif
1436 {
1437   char digits[40];
1438   int len, ix;
1439
1440   if (space) (*out_char) (' ');
1441   sprintf (digits, "%ld", val);
1442   len = strlen (digits);
1443   while (size > len)
1444     {
1445       (*out_char) ('0');
1446       size--;
1447     }
1448   for (ix=0; ix < len; ix++)
1449     (*out_char) (digits[ix]);
1450 }
1451
1452 /* Output of a bcd number.  NUM is written in base O_BASE using OUT_CHAR
1453    as the routine to do the actual output of the characters. */
1454
1455 void
1456 bc_out_num (num, o_base, out_char, leading_zero)
1457      bc_num num;
1458      int o_base;
1459 #ifdef __STDC__
1460      void (*out_char)(int);
1461 #else
1462      void (*out_char)();
1463 #endif
1464      int leading_zero;
1465 {
1466   char *nptr;
1467   int  index, fdigit, pre_space;
1468   stk_rec *digits, *temp;
1469   bc_num int_part, frac_part, base, cur_dig, t_num, max_o_digit;
1470
1471   /* The negative sign if needed. */
1472   if (num->n_sign == MINUS) (*out_char) ('-');
1473
1474   /* Output the number. */
1475   if (bc_is_zero (num))
1476     (*out_char) ('0');
1477   else
1478     if (o_base == 10)
1479       {
1480         /* The number is in base 10, do it the fast way. */
1481         nptr = num->n_value;
1482         if (num->n_len > 1 || *nptr != 0)
1483           for (index=num->n_len; index>0; index--)
1484             (*out_char) (BCD_CHAR(*nptr++));
1485         else
1486           nptr++;
1487
1488         if (leading_zero && bc_is_zero (num))
1489           (*out_char) ('0');
1490
1491         /* Now the fraction. */
1492         if (num->n_scale > 0)
1493           {
1494             (*out_char) ('.');
1495             for (index=0; index<num->n_scale; index++)
1496               (*out_char) (BCD_CHAR(*nptr++));
1497           }
1498       }
1499     else
1500       {
1501         /* special case ... */
1502         if (leading_zero && bc_is_zero (num))
1503           (*out_char) ('0');
1504
1505         /* The number is some other base. */
1506         digits = NULL;
1507         bc_init_num (&int_part);
1508         bc_divide (num, _one_, &int_part, 0);
1509         bc_init_num (&frac_part);
1510         bc_init_num (&cur_dig);
1511         bc_init_num (&base);
1512         bc_sub (num, int_part, &frac_part, 0);
1513         /* Make the INT_PART and FRAC_PART positive. */
1514         int_part->n_sign = PLUS;
1515         frac_part->n_sign = PLUS;
1516         bc_int2num (&base, o_base);
1517         bc_init_num (&max_o_digit);
1518         bc_int2num (&max_o_digit, o_base-1);
1519
1520
1521         /* Get the digits of the integer part and push them on a stack. */
1522         while (!bc_is_zero (int_part))
1523           {
1524             bc_modulo (int_part, base, &cur_dig, 0);
1525             temp = (stk_rec *) malloc (sizeof(stk_rec));
1526             if (temp == NULL) bc_out_of_memory();
1527             temp->digit = bc_num2long (cur_dig);
1528             temp->next = digits;
1529             digits = temp;
1530             bc_divide (int_part, base, &int_part, 0);
1531           }
1532
1533         /* Print the digits on the stack. */
1534         if (digits != NULL)
1535           {
1536             /* Output the digits. */
1537             while (digits != NULL)
1538               {
1539                 temp = digits;
1540                 digits = digits->next;
1541                 if (o_base <= 16)
1542                   (*out_char) (ref_str[ (int) temp->digit]);
1543                 else
1544                   bc_out_long (temp->digit, max_o_digit->n_len, 1, out_char);
1545                 free (temp);
1546               }
1547           }
1548
1549         /* Get and print the digits of the fraction part. */
1550         if (num->n_scale > 0)
1551           {
1552             (*out_char) ('.');
1553             pre_space = 0;
1554             t_num = bc_copy_num (_one_);
1555             while (t_num->n_len <= num->n_scale) {
1556               bc_multiply (frac_part, base, &frac_part, num->n_scale);
1557               fdigit = bc_num2long (frac_part);
1558               bc_int2num (&int_part, fdigit);
1559               bc_sub (frac_part, int_part, &frac_part, 0);
1560               if (o_base <= 16)
1561                 (*out_char) (ref_str[fdigit]);
1562               else {
1563                 bc_out_long (fdigit, max_o_digit->n_len, pre_space, out_char);
1564                 pre_space = 1;
1565               }
1566               bc_multiply (t_num, base, &t_num, 0);
1567             }
1568             bc_free_num (&t_num);
1569           }
1570
1571         /* Clean up. */
1572         bc_free_num (&int_part);
1573         bc_free_num (&frac_part);
1574         bc_free_num (&base);
1575         bc_free_num (&cur_dig);
1576         bc_free_num (&max_o_digit);
1577       }
1578 }
1579 /* Convert a number NUM to a long.  The function returns only the integer
1580    part of the number.  For numbers that are too large to represent as
1581    a long, this function returns a zero.  This can be detected by checking
1582    the NUM for zero after having a zero returned. */
1583
1584 long
1585 bc_num2long (num)
1586      bc_num num;
1587 {
1588   long val;
1589   char *nptr;
1590   int  index;
1591
1592   /* Extract the int value, ignore the fraction. */
1593   val = 0;
1594   nptr = num->n_value;
1595   for (index=num->n_len; (index>0) && (val<=(LONG_MAX/BASE)); index--)
1596     val = val*BASE + *nptr++;
1597
1598   /* Check for overflow.  If overflow, return zero. */
1599   if (index>0) val = 0;
1600   if (val < 0) val = 0;
1601
1602   /* Return the value. */
1603   if (num->n_sign == PLUS)
1604     return (val);
1605   else
1606     return (-val);
1607 }
1608
1609
1610 /* Convert an integer VAL to a bc number NUM. */
1611
1612 void
1613 bc_int2num (num, val)
1614      bc_num *num;
1615      int val;
1616 {
1617   char buffer[30];
1618   char *bptr, *vptr;
1619   int  ix = 1;
1620   char neg = 0;
1621
1622   /* Sign. */
1623   if (val < 0)
1624     {
1625       neg = 1;
1626       val = -val;
1627     }
1628
1629   /* Get things going. */
1630   bptr = buffer;
1631   *bptr++ = val % BASE;
1632   val = val / BASE;
1633
1634   /* Extract remaining digits. */
1635   while (val != 0)
1636     {
1637       *bptr++ = val % BASE;
1638       val = val / BASE;
1639       ix++;             /* Count the digits. */
1640     }
1641
1642   /* Make the number. */
1643   bc_free_num (num);
1644   *num = bc_new_num (ix, 0);
1645   if (neg) (*num)->n_sign = MINUS;
1646
1647   /* Assign the digits. */
1648   vptr = (*num)->n_value;
1649   while (ix-- > 0)
1650     *vptr++ = *--bptr;
1651 }
1652
1653 /* Convert a numbers to a string.  Base 10 only.*/
1654
1655 char
1656 *num2str (num)
1657       bc_num num;
1658 {
1659   char *str, *sptr;
1660   char *nptr;
1661   int  index, signch;
1662
1663   /* Allocate the string memory. */
1664   signch = ( num->n_sign == PLUS ? 0 : 1 );  /* Number of sign chars. */
1665   if (num->n_scale > 0)
1666     str = (char *) malloc (num->n_len + num->n_scale + 2 + signch);
1667   else
1668     str = (char *) malloc (num->n_len + 1 + signch);
1669   if (str == NULL) bc_out_of_memory();
1670
1671   /* The negative sign if needed. */
1672   sptr = str;
1673   if (signch) *sptr++ = '-';
1674
1675   /* Load the whole number. */
1676   nptr = num->n_value;
1677   for (index=num->n_len; index>0; index--)
1678     *sptr++ = BCD_CHAR(*nptr++);
1679
1680   /* Now the fraction. */
1681   if (num->n_scale > 0)
1682     {
1683       *sptr++ = '.';
1684       for (index=0; index<num->n_scale; index++)
1685         *sptr++ = BCD_CHAR(*nptr++);
1686     }
1687
1688   /* Terminate the string and return it! */
1689   *sptr = '\0';
1690   return (str);
1691 }
1692 /* Convert strings to bc numbers.  Base 10 only.*/
1693
1694 void
1695 bc_str2num (num, str, scale)
1696      bc_num *num;
1697      char *str;
1698      int scale;
1699 {
1700   int digits, strscale;
1701   char *ptr, *nptr;
1702   char zero_int;
1703
1704   /* Prepare num. */
1705   bc_free_num (num);
1706
1707   /* Check for valid number and count digits. */
1708   ptr = str;
1709   digits = 0;
1710   strscale = 0;
1711   zero_int = FALSE;
1712   if ( (*ptr == '+') || (*ptr == '-'))  ptr++;  /* Sign */
1713   while (*ptr == '0') ptr++;                    /* Skip leading zeros. */
1714   while (isdigit((int)*ptr)) ptr++, digits++;   /* digits */
1715   if (*ptr == '.') ptr++;                       /* decimal point */
1716   while (isdigit((int)*ptr)) ptr++, strscale++; /* digits */
1717   if ((*ptr != '\0') || (digits+strscale == 0))
1718     {
1719       *num = bc_copy_num (_zero_);
1720       return;
1721     }
1722
1723   /* Adjust numbers and allocate storage and initialize fields. */
1724   strscale = MIN(strscale, scale);
1725   if (digits == 0)
1726     {
1727       zero_int = TRUE;
1728       digits = 1;
1729     }
1730   *num = bc_new_num (digits, strscale);
1731
1732   /* Build the whole number. */
1733   ptr = str;
1734   if (*ptr == '-')
1735     {
1736       (*num)->n_sign = MINUS;
1737       ptr++;
1738     }
1739   else
1740     {
1741       (*num)->n_sign = PLUS;
1742       if (*ptr == '+') ptr++;
1743     }
1744   while (*ptr == '0') ptr++;                    /* Skip leading zeros. */
1745   nptr = (*num)->n_value;
1746   if (zero_int)
1747     {
1748       *nptr++ = 0;
1749       digits = 0;
1750     }
1751   for (;digits > 0; digits--)
1752     *nptr++ = CH_VAL(*ptr++);
1753
1754
1755   /* Build the fractional part. */
1756   if (strscale > 0)
1757     {
1758       ptr++;  /* skip the decimal point! */
1759       for (;strscale > 0; strscale--)
1760         *nptr++ = CH_VAL(*ptr++);
1761     }
1762 }
1763
1764 /* pn prints the number NUM in base 10. */
1765
1766 static void
1767 out_char (int c)
1768 {
1769   putchar(c);
1770 }
1771
1772
1773 void
1774 pn (num)
1775      bc_num num;
1776 {
1777   bc_out_num (num, 10, out_char, 0);
1778   out_char ('\n');
1779 }
1780
1781
1782 /* pv prints a character array as if it was a string of bcd digits. */
1783 void
1784 pv (name, num, len)
1785      char *name;
1786      unsigned char *num;
1787      int len;
1788 {
1789   int i;
1790   printf ("%s=", name);
1791   for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i]));
1792   printf ("\n");
1793 }