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. */
48bc_num _zero_;
49bc_num _one_;
50bc_num _two_;
51
52static bc_num _bc_Free_list = NULL;
53
54/* new_num allocates a number and sets fields to known values. */
55
56bc_num
57bc_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
83void
84bc_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
101void
102bc_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
114bc_num
115bc_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
125void
126bc_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
136static 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
152static 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
261int
262bc_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
270char
271bc_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
279char
280bc_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
306char
307bc_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
336static 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
428static 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
528void
529bc_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
578void
579bc_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
628int mul_base_digits = MUL_BASE_DIGITS;
629#define MUL_SMALL_DIGITS mul_base_digits/4
630
631/* Multiply utility routines */
632
633static bc_num
634new_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
656static 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. *) */
691static 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*/
757static 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
856void
857bc_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
891static 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
931int
932bc_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
1127int
1128bc_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
1164int
1165bc_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
1176int
1177bc_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
1236void
1237bc_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
1320int
1321bc_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. */
1413typedef struct stk_rec {
1414	long  digit;
1415	struct stk_rec *next;
1416} stk_rec;
1417
1418/* The reference string for digits. */
1419static 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
1427void
1428bc_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
1455void
1456bc_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
1584long
1585bc_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
1612void
1613bc_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
1655char
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
1694void
1695bc_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
1766static void
1767out_char (int c)
1768{
1769  putchar(c);
1770}
1771
1772
1773void
1774pn (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. */
1783void
1784pv (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}
1794