contrib/bc/lib/number.c

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