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