"Fossies" - the Fresh Open Source Software Archive

Member "bc-1.06.95/lib/number.c" (5 Sep 2006, 41529 Bytes) of archive /linux/misc/old/bc-1.06.95.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    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       51 Franklin Street, Fifth Floor
   20       Boston, MA 02110-1301  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 #ifdef HAVE_STDLIB_H
   37 #include <stdlib.h>
   38 #endif
   39 #ifdef HAVE_STRING_H
   40 #include <string.h>
   41 #endif
   42 #include <ctype.h>
   43 
   44 /* Prototypes needed for external utility routines. */
   45 
   46 #define bc_rt_warn rt_warn
   47 #define bc_rt_error rt_error
   48 #define bc_out_of_memory out_of_memory
   49 
   50 _PROTOTYPE(void rt_warn, (const char *mesg ,...));
   51 _PROTOTYPE(void rt_error, (const char *mesg ,...));
   52 _PROTOTYPE(void out_of_memory, (void));
   53 
   54 /* Storage used for special numbers. */
   55 bc_num _zero_;
   56 bc_num _one_;
   57 bc_num _two_;
   58 
   59 static bc_num _bc_Free_list = NULL;
   60 
   61 /* new_num allocates a number and sets fields to known values. */
   62 
   63 bc_num
   64 bc_new_num (length, scale)
   65      int length, scale;
   66 {
   67   bc_num temp;
   68 
   69   if (_bc_Free_list != NULL) {
   70     temp = _bc_Free_list;
   71     _bc_Free_list = temp->n_next;
   72   } else {
   73     temp = (bc_num) malloc (sizeof(bc_struct));
   74     if (temp == NULL) bc_out_of_memory ();
   75   }
   76   temp->n_sign = PLUS;
   77   temp->n_len = length;
   78   temp->n_scale = scale;
   79   temp->n_refs = 1;
   80   temp->n_ptr = (char *) malloc (length+scale);
   81   if (temp->n_ptr == NULL) bc_out_of_memory();
   82   temp->n_value = temp->n_ptr;
   83   memset (temp->n_ptr, 0, length+scale);
   84   return temp;
   85 }
   86 
   87 /* "Frees" a bc_num NUM.  Actually decreases reference count and only
   88    frees the storage if reference count is zero. */
   89 
   90 void
   91 bc_free_num (num)
   92     bc_num *num;
   93 {
   94   if (*num == NULL) return;
   95   (*num)->n_refs--;
   96   if ((*num)->n_refs == 0) {
   97     if ((*num)->n_ptr)
   98       free ((*num)->n_ptr);
   99     (*num)->n_next = _bc_Free_list;
  100     _bc_Free_list = *num;
  101   }
  102   *num = NULL;
  103 }
  104 
  105 
  106 /* Intitialize the number package! */
  107 
  108 void
  109 bc_init_numbers ()
  110 {
  111   _zero_ = bc_new_num (1,0);
  112   _one_  = bc_new_num (1,0);
  113   _one_->n_value[0] = 1;
  114   _two_  = bc_new_num (1,0);
  115   _two_->n_value[0] = 2;
  116 }
  117 
  118 
  119 /* Make a copy of a number!  Just increments the reference count! */
  120 
  121 bc_num
  122 bc_copy_num (num)
  123      bc_num num;
  124 {
  125   num->n_refs++;
  126   return num;
  127 }
  128 
  129 
  130 /* Initialize a number NUM by making it a copy of zero. */
  131 
  132 void
  133 bc_init_num (num)
  134      bc_num *num;
  135 {
  136   *num = bc_copy_num (_zero_);
  137 }
  138 
  139 /* For many things, we may have leading zeros in a number NUM.
  140    _bc_rm_leading_zeros just moves the data "value" pointer to the
  141    correct place and adjusts the length. */
  142 
  143 static void
  144 _bc_rm_leading_zeros (num)
  145      bc_num num;
  146 {
  147   /* We can move n_value to point to the first non zero digit! */
  148   while (*num->n_value == 0 && num->n_len > 1) {
  149     num->n_value++;
  150     num->n_len--;
  151   }
  152 }
  153 
  154 
  155 /* Compare two bc numbers.  Return value is 0 if equal, -1 if N1 is less
  156    than N2 and +1 if N1 is greater than N2.  If USE_SIGN is false, just
  157    compare the magnitudes. */
  158 
  159 static int
  160 _bc_do_compare (n1, n2, use_sign, ignore_last)
  161      bc_num n1, n2;
  162      int use_sign;
  163      int ignore_last;
  164 {
  165   char *n1ptr, *n2ptr;
  166   int  count;
  167 
  168   /* First, compare signs. */
  169   if (use_sign && n1->n_sign != n2->n_sign)
  170     {
  171       if (n1->n_sign == PLUS)
  172     return (1); /* Positive N1 > Negative N2 */
  173       else
  174     return (-1);    /* Negative N1 < Positive N1 */
  175     }
  176 
  177   /* Now compare the magnitude. */
  178   if (n1->n_len != n2->n_len)
  179     {
  180       if (n1->n_len > n2->n_len)
  181     {
  182       /* Magnitude of n1 > n2. */
  183       if (!use_sign || n1->n_sign == PLUS)
  184         return (1);
  185       else
  186         return (-1);
  187     }
  188       else
  189     {
  190       /* Magnitude of n1 < n2. */
  191       if (!use_sign || n1->n_sign == PLUS)
  192         return (-1);
  193       else
  194         return (1);
  195     }
  196     }
  197 
  198   /* If we get here, they have the same number of integer digits.
  199      check the integer part and the equal length part of the fraction. */
  200   count = n1->n_len + MIN (n1->n_scale, n2->n_scale);
  201   n1ptr = n1->n_value;
  202   n2ptr = n2->n_value;
  203 
  204   while ((count > 0) && (*n1ptr == *n2ptr))
  205     {
  206       n1ptr++;
  207       n2ptr++;
  208       count--;
  209     }
  210   if (ignore_last && count == 1 && n1->n_scale == n2->n_scale)
  211     return (0);
  212   if (count != 0)
  213     {
  214       if (*n1ptr > *n2ptr)
  215     {
  216       /* Magnitude of n1 > n2. */
  217       if (!use_sign || n1->n_sign == PLUS)
  218         return (1);
  219       else
  220         return (-1);
  221     }
  222       else
  223     {
  224       /* Magnitude of n1 < n2. */
  225       if (!use_sign || n1->n_sign == PLUS)
  226         return (-1);
  227       else
  228         return (1);
  229     }
  230     }
  231 
  232   /* They are equal up to the last part of the equal part of the fraction. */
  233   if (n1->n_scale != n2->n_scale)
  234     {
  235       if (n1->n_scale > n2->n_scale)
  236     {
  237       for (count = n1->n_scale-n2->n_scale; count>0; count--)
  238         if (*n1ptr++ != 0)
  239           {
  240         /* Magnitude of n1 > n2. */
  241         if (!use_sign || n1->n_sign == PLUS)
  242           return (1);
  243         else
  244           return (-1);
  245           }
  246     }
  247       else
  248     {
  249       for (count = n2->n_scale-n1->n_scale; count>0; count--)
  250         if (*n2ptr++ != 0)
  251           {
  252         /* Magnitude of n1 < n2. */
  253         if (!use_sign || n1->n_sign == PLUS)
  254           return (-1);
  255         else
  256           return (1);
  257           }
  258     }
  259     }
  260 
  261   /* They must be equal! */
  262   return (0);
  263 }
  264 
  265 
  266 /* This is the "user callable" routine to compare numbers N1 and N2. */
  267 
  268 int
  269 bc_compare (n1, n2)
  270      bc_num n1, n2;
  271 {
  272   return _bc_do_compare (n1, n2, TRUE, FALSE);
  273 }
  274 
  275 /* In some places we need to check if the number is negative. */
  276 
  277 char
  278 bc_is_neg (num)
  279      bc_num num;
  280 {
  281   return num->n_sign == MINUS;
  282 }
  283 
  284 /* In some places we need to check if the number NUM is zero. */
  285 
  286 char
  287 bc_is_zero (num)
  288      bc_num num;
  289 {
  290   int  count;
  291   char *nptr;
  292 
  293   /* Quick check. */
  294   if (num == _zero_) return TRUE;
  295 
  296   /* Initialize */
  297   count = num->n_len + num->n_scale;
  298   nptr = num->n_value;
  299 
  300   /* The check */
  301   while ((count > 0) && (*nptr++ == 0)) count--;
  302 
  303   if (count != 0)
  304     return FALSE;
  305   else
  306     return TRUE;
  307 }
  308 
  309 /* In some places we need to check if the number NUM is almost zero.
  310    Specifically, all but the last digit is 0 and the last digit is 1.
  311    Last digit is defined by scale. */
  312 
  313 char
  314 bc_is_near_zero (num, scale)
  315      bc_num num;
  316      int scale;
  317 {
  318   int  count;
  319   char *nptr;
  320 
  321   /* Error checking */
  322   if (scale > num->n_scale)
  323     scale = num->n_scale;
  324 
  325   /* Initialize */
  326   count = num->n_len + scale;
  327   nptr = num->n_value;
  328 
  329   /* The check */
  330   while ((count > 0) && (*nptr++ == 0)) count--;
  331 
  332   if (count != 0 && (count != 1 || *--nptr != 1))
  333     return FALSE;
  334   else
  335     return TRUE;
  336 }
  337 
  338 
  339 /* Perform addition: N1 is added to N2 and the value is
  340    returned.  The signs of N1 and N2 are ignored.
  341    SCALE_MIN is to set the minimum scale of the result. */
  342 
  343 static bc_num
  344 _bc_do_add (n1, n2, scale_min)
  345      bc_num n1, n2;
  346      int scale_min;
  347 {
  348   bc_num sum;
  349   int sum_scale, sum_digits;
  350   char *n1ptr, *n2ptr, *sumptr;
  351   int carry, n1bytes, n2bytes;
  352   int count;
  353 
  354   /* Prepare sum. */
  355   sum_scale = MAX (n1->n_scale, n2->n_scale);
  356   sum_digits = MAX (n1->n_len, n2->n_len) + 1;
  357   sum = bc_new_num (sum_digits, MAX(sum_scale, scale_min));
  358 
  359   /* Zero extra digits made by scale_min. */
  360   if (scale_min > sum_scale)
  361     {
  362       sumptr = (char *) (sum->n_value + sum_scale + sum_digits);
  363       for (count = scale_min - sum_scale; count > 0; count--)
  364     *sumptr++ = 0;
  365     }
  366 
  367   /* Start with the fraction part.  Initialize the pointers. */
  368   n1bytes = n1->n_scale;
  369   n2bytes = n2->n_scale;
  370   n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1);
  371   n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1);
  372   sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1);
  373 
  374   /* Add the fraction part.  First copy the longer fraction.*/
  375   if (n1bytes != n2bytes)
  376     {
  377       if (n1bytes > n2bytes)
  378     while (n1bytes>n2bytes)
  379       { *sumptr-- = *n1ptr--; n1bytes--;}
  380       else
  381     while (n2bytes>n1bytes)
  382       { *sumptr-- = *n2ptr--; n2bytes--;}
  383     }
  384 
  385   /* Now add the remaining fraction part and equal size integer parts. */
  386   n1bytes += n1->n_len;
  387   n2bytes += n2->n_len;
  388   carry = 0;
  389   while ((n1bytes > 0) && (n2bytes > 0))
  390     {
  391       *sumptr = *n1ptr-- + *n2ptr-- + carry;
  392       if (*sumptr > (BASE-1))
  393     {
  394        carry = 1;
  395        *sumptr -= BASE;
  396     }
  397       else
  398     carry = 0;
  399       sumptr--;
  400       n1bytes--;
  401       n2bytes--;
  402     }
  403 
  404   /* Now add carry the longer integer part. */
  405   if (n1bytes == 0)
  406     { n1bytes = n2bytes; n1ptr = n2ptr; }
  407   while (n1bytes-- > 0)
  408     {
  409       *sumptr = *n1ptr-- + carry;
  410       if (*sumptr > (BASE-1))
  411     {
  412        carry = 1;
  413        *sumptr -= BASE;
  414      }
  415       else
  416     carry = 0;
  417       sumptr--;
  418     }
  419 
  420   /* Set final carry. */
  421   if (carry == 1)
  422     *sumptr += 1;
  423 
  424   /* Adjust sum and return. */
  425   _bc_rm_leading_zeros (sum);
  426   return sum;
  427 }
  428 
  429 
  430 /* Perform subtraction: N2 is subtracted from N1 and the value is
  431    returned.  The signs of N1 and N2 are ignored.  Also, N1 is
  432    assumed to be larger than N2.  SCALE_MIN is the minimum scale
  433    of the result. */
  434 
  435 static bc_num
  436 _bc_do_sub (n1, n2, scale_min)
  437      bc_num n1, n2;
  438      int scale_min;
  439 {
  440   bc_num diff;
  441   int diff_scale, diff_len;
  442   int min_scale, min_len;
  443   char *n1ptr, *n2ptr, *diffptr;
  444   int borrow, count, val;
  445 
  446   /* Allocate temporary storage. */
  447   diff_len = MAX (n1->n_len, n2->n_len);
  448   diff_scale = MAX (n1->n_scale, n2->n_scale);
  449   min_len = MIN  (n1->n_len, n2->n_len);
  450   min_scale = MIN (n1->n_scale, n2->n_scale);
  451   diff = bc_new_num (diff_len, MAX(diff_scale, scale_min));
  452 
  453   /* Zero extra digits made by scale_min. */
  454   if (scale_min > diff_scale)
  455     {
  456       diffptr = (char *) (diff->n_value + diff_len + diff_scale);
  457       for (count = scale_min - diff_scale; count > 0; count--)
  458     *diffptr++ = 0;
  459     }
  460 
  461   /* Initialize the subtract. */
  462   n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1);
  463   n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1);
  464   diffptr = (char *) (diff->n_value + diff_len + diff_scale -1);
  465 
  466   /* Subtract the numbers. */
  467   borrow = 0;
  468 
  469   /* Take care of the longer scaled number. */
  470   if (n1->n_scale != min_scale)
  471     {
  472       /* n1 has the longer scale */
  473       for (count = n1->n_scale - min_scale; count > 0; count--)
  474     *diffptr-- = *n1ptr--;
  475     }
  476   else
  477     {
  478       /* n2 has the longer scale */
  479       for (count = n2->n_scale - min_scale; count > 0; count--)
  480     {
  481       val = - *n2ptr-- - borrow;
  482       if (val < 0)
  483         {
  484           val += BASE;
  485           borrow = 1;
  486         }
  487       else
  488         borrow = 0;
  489       *diffptr-- = val;
  490     }
  491     }
  492 
  493   /* Now do the equal length scale and integer parts. */
  494 
  495   for (count = 0; count < min_len + min_scale; count++)
  496     {
  497       val = *n1ptr-- - *n2ptr-- - borrow;
  498       if (val < 0)
  499     {
  500       val += BASE;
  501       borrow = 1;
  502     }
  503       else
  504     borrow = 0;
  505       *diffptr-- = val;
  506     }
  507 
  508   /* If n1 has more digits then n2, we now do that subtract. */
  509   if (diff_len != min_len)
  510     {
  511       for (count = diff_len - min_len; count > 0; count--)
  512     {
  513       val = *n1ptr-- - borrow;
  514       if (val < 0)
  515         {
  516           val += BASE;
  517           borrow = 1;
  518         }
  519       else
  520         borrow = 0;
  521       *diffptr-- = val;
  522     }
  523     }
  524 
  525   /* Clean up and return. */
  526   _bc_rm_leading_zeros (diff);
  527   return diff;
  528 }
  529 
  530 
  531 /* Here is the full subtract routine that takes care of negative numbers.
  532    N2 is subtracted from N1 and the result placed in RESULT.  SCALE_MIN
  533    is the minimum scale for the result. */
  534 
  535 void
  536 bc_sub (n1, n2, result, scale_min)
  537      bc_num n1, n2, *result;
  538      int scale_min;
  539 {
  540   bc_num diff = NULL;
  541   int cmp_res;
  542   int res_scale;
  543 
  544   if (n1->n_sign != n2->n_sign)
  545     {
  546       diff = _bc_do_add (n1, n2, scale_min);
  547       diff->n_sign = n1->n_sign;
  548     }
  549   else
  550     {
  551       /* subtraction must be done. */
  552       /* Compare magnitudes. */
  553       cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);
  554       switch (cmp_res)
  555     {
  556     case -1:
  557       /* n1 is less than n2, subtract n1 from n2. */
  558       diff = _bc_do_sub (n2, n1, scale_min);
  559       diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS);
  560       break;
  561     case  0:
  562       /* They are equal! return zero! */
  563       res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
  564       diff = bc_new_num (1, res_scale);
  565       memset (diff->n_value, 0, res_scale+1);
  566       break;
  567     case  1:
  568       /* n2 is less than n1, subtract n2 from n1. */
  569       diff = _bc_do_sub (n1, n2, scale_min);
  570       diff->n_sign = n1->n_sign;
  571       break;
  572     }
  573     }
  574 
  575   /* Clean up and return. */
  576   bc_free_num (result);
  577   *result = diff;
  578 }
  579 
  580 
  581 /* Here is the full add routine that takes care of negative numbers.
  582    N1 is added to N2 and the result placed into RESULT.  SCALE_MIN
  583    is the minimum scale for the result. */
  584 
  585 void
  586 bc_add (n1, n2, result, scale_min)
  587      bc_num n1, n2, *result;
  588      int scale_min;
  589 {
  590   bc_num sum = NULL;
  591   int cmp_res;
  592   int res_scale;
  593 
  594   if (n1->n_sign == n2->n_sign)
  595     {
  596       sum = _bc_do_add (n1, n2, scale_min);
  597       sum->n_sign = n1->n_sign;
  598     }
  599   else
  600     {
  601       /* subtraction must be done. */
  602       cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);  /* Compare magnitudes. */
  603       switch (cmp_res)
  604     {
  605     case -1:
  606       /* n1 is less than n2, subtract n1 from n2. */
  607       sum = _bc_do_sub (n2, n1, scale_min);
  608       sum->n_sign = n2->n_sign;
  609       break;
  610     case  0:
  611       /* They are equal! return zero with the correct scale! */
  612       res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
  613       sum = bc_new_num (1, res_scale);
  614       memset (sum->n_value, 0, res_scale+1);
  615       break;
  616     case  1:
  617       /* n2 is less than n1, subtract n2 from n1. */
  618       sum = _bc_do_sub (n1, n2, scale_min);
  619       sum->n_sign = n1->n_sign;
  620     }
  621     }
  622 
  623   /* Clean up and return. */
  624   bc_free_num (result);
  625   *result = sum;
  626 }
  627 
  628 /* Recursive vs non-recursive multiply crossover ranges. */
  629 #if defined(MULDIGITS)
  630 #include "muldigits.h"
  631 #else
  632 #define MUL_BASE_DIGITS 80
  633 #endif
  634 
  635 int mul_base_digits = MUL_BASE_DIGITS;
  636 #define MUL_SMALL_DIGITS mul_base_digits/4
  637 
  638 /* Multiply utility routines */
  639 
  640 static bc_num
  641 new_sub_num (length, scale, value)
  642      int length, scale;
  643      char *value;
  644 {
  645   bc_num temp;
  646 
  647   if (_bc_Free_list != NULL) {
  648     temp = _bc_Free_list;
  649     _bc_Free_list = temp->n_next;
  650   } else {
  651     temp = (bc_num) malloc (sizeof(bc_struct));
  652     if (temp == NULL) bc_out_of_memory ();
  653   }
  654   temp->n_sign = PLUS;
  655   temp->n_len = length;
  656   temp->n_scale = scale;
  657   temp->n_refs = 1;
  658   temp->n_ptr = NULL;
  659   temp->n_value = value;
  660   return temp;
  661 }
  662 
  663 static void
  664 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod)
  665 {
  666   char *n1ptr, *n2ptr, *pvptr;
  667   char *n1end, *n2end;      /* To the end of n1 and n2. */
  668   int indx, sum, prodlen;
  669 
  670   prodlen = n1len+n2len+1;
  671 
  672   *prod = bc_new_num (prodlen, 0);
  673 
  674   n1end = (char *) (n1->n_value + n1len - 1);
  675   n2end = (char *) (n2->n_value + n2len - 1);
  676   pvptr = (char *) ((*prod)->n_value + prodlen - 1);
  677   sum = 0;
  678 
  679   /* Here is the loop... */
  680   for (indx = 0; indx < prodlen-1; indx++)
  681     {
  682       n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
  683       n2ptr = (char *) (n2end - MIN(indx, n2len-1));
  684       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
  685     sum += *n1ptr-- * *n2ptr++;
  686       *pvptr-- = sum % BASE;
  687       sum = sum / BASE;
  688     }
  689   *pvptr = sum;
  690 }
  691 
  692 
  693 /* A special adder/subtractor for the recursive divide and conquer
  694    multiply algorithm.  Note: if sub is called, accum must
  695    be larger that what is being subtracted.  Also, accum and val
  696    must have n_scale = 0.  (e.g. they must look like integers. *) */
  697 static void
  698 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
  699 {
  700   signed char *accp, *valp;
  701   int  count, carry;
  702 
  703   count = val->n_len;
  704   if (val->n_value[0] == 0)
  705     count--;
  706   assert (accum->n_len+accum->n_scale >= shift+count);
  707   
  708   /* Set up pointers and others */
  709   accp = (signed char *)(accum->n_value +
  710              accum->n_len + accum->n_scale - shift - 1);
  711   valp = (signed char *)(val->n_value + val->n_len - 1);
  712   carry = 0;
  713 
  714   if (sub) {
  715     /* Subtraction, carry is really borrow. */
  716     while (count--) {
  717       *accp -= *valp-- + carry;
  718       if (*accp < 0) {
  719     carry = 1;
  720         *accp-- += BASE;
  721       } else {
  722     carry = 0;
  723     accp--;
  724       }
  725     }
  726     while (carry) {
  727       *accp -= carry;
  728       if (*accp < 0)
  729     *accp-- += BASE;
  730       else
  731     carry = 0;
  732     }
  733   } else {
  734     /* Addition */
  735     while (count--) {
  736       *accp += *valp-- + carry;
  737       if (*accp > (BASE-1)) {
  738     carry = 1;
  739         *accp-- -= BASE;
  740       } else {
  741     carry = 0;
  742     accp--;
  743       }
  744     }
  745     while (carry) {
  746       *accp += carry;
  747       if (*accp > (BASE-1))
  748     *accp-- -= BASE;
  749       else
  750     carry = 0;
  751     }
  752   }
  753 }
  754 
  755 /* Recursive divide and conquer multiply algorithm.  
  756    Based on 
  757    Let u = u0 + u1*(b^n)
  758    Let v = v0 + v1*(b^n)
  759    Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
  760 
  761    B is the base of storage, number of digits in u1,u0 close to equal.
  762 */
  763 static void
  764 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod)
  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);
  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);
  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);
  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);
  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 
  861 void
  862 bc_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);
  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 
  896 static 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 
  936 int
  937 bc_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 
 1132 int
 1133 bc_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 
 1169 int
 1170 bc_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 
 1181 int
 1182 bc_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 
 1241 void
 1242 bc_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 
 1325 int
 1326 bc_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. */
 1418 typedef struct stk_rec {
 1419     long  digit;
 1420     struct stk_rec *next;
 1421 } stk_rec;
 1422 
 1423 /* The reference string for digits. */
 1424 static 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 
 1432 void
 1433 bc_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 
 1460 void
 1461 bc_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  ix, 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 (ix=num->n_len; ix>0; ix--)
 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 (ix=0; ix<num->n_scale; ix++)
 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 
 1589 long
 1590 bc_num2long (num)
 1591      bc_num num;
 1592 {
 1593   long val;
 1594   char *nptr;
 1595   int  i;
 1596 
 1597   /* Extract the int value, ignore the fraction. */
 1598   val = 0;
 1599   nptr = num->n_value;
 1600   for (i=num->n_len; (i>0) && (val<=(LONG_MAX/BASE)); i--)
 1601     val = val*BASE + *nptr++;
 1602 
 1603   /* Check for overflow.  If overflow, return zero. */
 1604   if (i>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 
 1617 void
 1618 bc_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 
 1660 char
 1661 *bc_num2str (num)
 1662       bc_num num;
 1663 {
 1664   char *str, *sptr;
 1665   char *nptr;
 1666   int  i, 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 (i=num->n_len; i>0; i--)
 1683     *sptr++ = BCD_CHAR(*nptr++);
 1684 
 1685   /* Now the fraction. */
 1686   if (num->n_scale > 0)
 1687     {
 1688       *sptr++ = '.';
 1689       for (i=0; i<num->n_scale; i++)
 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 
 1699 void
 1700 bc_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 /* Debugging routines */
 1770 
 1771 #ifdef DEBUG
 1772 
 1773 /* pn prints the number NUM in base 10. */
 1774 
 1775 static void
 1776 out_char (int c)
 1777 {
 1778   putchar(c);
 1779 }
 1780 
 1781 
 1782 void
 1783 pn (num)
 1784      bc_num num;
 1785 {
 1786   bc_out_num (num, 10, out_char, 0);
 1787   out_char ('\n');
 1788 }
 1789 
 1790 
 1791 /* pv prints a character array as if it was a string of bcd digits. */
 1792 void
 1793 pv (name, num, len)
 1794      char *name;
 1795      unsigned char *num;
 1796      int len;
 1797 {
 1798   int i;
 1799   printf ("%s=", name);
 1800   for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i]));
 1801   printf ("\n");
 1802 }
 1803 
 1804 #endif