gnu/usr.bin/bc/number.c

   1 /* number.c: Implements arbitrary precision numbers. */
   2
   3 /*  This file is part of bc written for MINIX.
   4     Copyright (C) 1991, 1992 Free Software Foundation, Inc.
   5
   6     This program is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License , or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program; see the file COPYING.  If not, write to
  18     the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  19
  20     You may contact the author by:
  21        e-mail:  phil@cs.wwu.edu
  22       us-mail:  Philip A. Nelson
  23                 Computer Science Department, 9062
  24                 Western Washington University
  25                 Bellingham, WA 98226-9062
  26
  27 *************************************************************************/
  28
  29 #include "bcdefs.h"
  30 #include "proto.h"
  31
  32 /* Storage used for special numbers. */
  33 bc_num _zero_;
  34 bc_num _one_;
  35 bc_num _two_;
  36
  37
  38 /* "Frees" a bc_num NUM.  Actually decreases reference count and only
  39    frees the storage if reference count is zero. */
  40
  41 void
  42 free_num (num)
  43     bc_num *num;
  44 {
  45   if (*num == NULL) return;
  46   (*num)->n_refs--;
  47   if ((*num)->n_refs == 0) free(*num);
  48   *num = NULL;
  49 }
  50
  51
  52 /* new_num allocates a number and sets fields to known values. */
  53
  54 bc_num
  55 new_num (length, scale)
  56      int length, scale;
  57 {
  58   bc_num temp;
  59
  60   temp = (bc_num) malloc (sizeof(bc_struct)+length+scale);
  61   if (temp == NULL) out_of_memory ();
  62   temp->n_sign = PLUS;
  63   temp->n_len = length;
  64   temp->n_scale = scale;
  65   temp->n_refs = 1;
  66   temp->n_value[0] = 0;
  67   return temp;
  68 }
  69
  70
  71 /* Intitialize the number package! */
  72
  73 void
  74 init_numbers ()
  75 {
  76   _zero_ = new_num (1,0);
  77   _one_  = new_num (1,0);
  78   _one_->n_value[0] = 1;
  79   _two_  = new_num (1,0);
  80   _two_->n_value[0] = 2;
  81 }
  82
  83
  84 /* Make a copy of a number!  Just increments the reference count! */
  85
  86 bc_num
  87 copy_num (num)
  88      bc_num num;
  89 {
  90   num->n_refs++;
  91   return num;
  92 }
  93
  94
  95 /* Initialize a number NUM by making it a copy of zero. */
  96
  97 void
  98 init_num (num)
  99      bc_num *num;
 100 {
 101   *num = copy_num (_zero_);
 102 }
 103
 104
 105 /* Convert an integer VAL to a bc number NUM. */
 106
 107 void
 108 int2num (num, val)
 109      bc_num *num;
 110      int val;
 111 {
 112   char buffer[30];
 113   char *bptr, *vptr;
 114   int  ix = 1;
 115   char neg = 0;
 116
 117   /* Sign. */
 118   if (val < 0)
 119     {
 120       neg = 1;
 121       val = -val;
 122     }
 123
 124   /* Get things going. */
 125   bptr = buffer;
 126   *bptr++ = val % 10;
 127   val = val / 10;
 128
 129   /* Extract remaining digits. */
 130   while (val != 0)
 131     {
 132       *bptr++ = val % 10;
 133       val = val / 10;
 134       ix++;             /* Count the digits. */
 135     }
 136
 137   /* Make the number. */
 138   free_num (num);
 139   *num = new_num (ix, 0);
 140   if (neg) (*num)->n_sign = MINUS;
 141
 142   /* Assign the digits. */
 143   vptr = (*num)->n_value;
 144   while (ix-- > 0)
 145     *vptr++ = *--bptr;
 146 }
 147
 148
 149 /* Convert a number NUM to a long.  The function returns only the integer
 150    part of the number.  For numbers that are too large to represent as
 151    a long, this function returns a zero.  This can be detected by checking
 152    the NUM for zero after having a zero returned. */
 153
 154 long
 155 num2long (num)
 156      bc_num num;
 157 {
 158   long val;
 159   char *nptr;
 160   int  index;
 161
 162   /* Extract the int value, ignore the fraction. */
 163   val = 0;
 164   nptr = num->n_value;
 165   for (index=num->n_len; (index>0) && (val<=(LONG_MAX/10)); index--)
 166     val = val*10 + *nptr++;
 167
 168   /* Check for overflow.  If overflow, return zero. */
 169   if (index>0) val = 0;
 170   if (val < 0) val = 0;
 171
 172   /* Return the value. */
 173   if (num->n_sign == PLUS)
 174     return (val);
 175   else
 176     return (-val);
 177 }
 178
 179
 180 /* The following are some math routines for numbers. */
 181 _PROTOTYPE(static int _do_compare, (bc_num n1, bc_num n2, int use_sign,
 182                                     int ignore_last));
 183 _PROTOTYPE(static void _rm_leading_zeros, (bc_num num));
 184 _PROTOTYPE(static bc_num _do_add, (bc_num n1, bc_num n2));
 185 _PROTOTYPE(static bc_num _do_sub, (bc_num n1, bc_num n2));
 186 _PROTOTYPE(static void _one_mult, (unsigned char *num, int size, int digit,
 187                                    unsigned char *result));
 188
 189
 190
 191 /* Compare two bc numbers.  Return value is 0 if equal, -1 if N1 is less
 192    than N2 and +1 if N1 is greater than N2.  If USE_SIGN is false, just
 193    compare the magnitudes. */
 194
 195 static int
 196 _do_compare (n1, n2, use_sign, ignore_last)
 197      bc_num n1, n2;
 198      int use_sign;
 199      int ignore_last;
 200 {
 201   char *n1ptr, *n2ptr;
 202   int  count;
 203
 204   /* First, compare signs. */
 205   if (use_sign && n1->n_sign != n2->n_sign)
 206     {
 207       if (n1->n_sign == PLUS)
 208         return (1);     /* Positive N1 > Negative N2 */
 209       else
 210         return (-1);    /* Negative N1 < Positive N1 */
 211     }
 212
 213   /* Now compare the magnitude. */
 214   if (n1->n_len != n2->n_len)
 215     {
 216       if (n1->n_len > n2->n_len)
 217         {
 218           /* Magnitude of n1 > n2. */
 219           if (!use_sign || n1->n_sign == PLUS)
 220             return (1);
 221           else
 222             return (-1);
 223         }
 224       else
 225         {
 226           /* Magnitude of n1 < n2. */
 227           if (!use_sign || n1->n_sign == PLUS)
 228             return (-1);
 229           else
 230             return (1);
 231         }
 232     }
 233
 234   /* If we get here, they have the same number of integer digits.
 235      check the integer part and the equal length part of the fraction. */
 236   count = n1->n_len + MIN (n1->n_scale, n2->n_scale);
 237   n1ptr = n1->n_value;
 238   n2ptr = n2->n_value;
 239
 240   while ((count > 0) && (*n1ptr == *n2ptr))
 241     {
 242       n1ptr++;
 243       n2ptr++;
 244       count--;
 245     }
 246   if (ignore_last && count == 1 && n1->n_scale == n2->n_scale)
 247     return (0);
 248   if (count != 0)
 249     {
 250       if (*n1ptr > *n2ptr)
 251         {
 252           /* Magnitude of n1 > n2. */
 253           if (!use_sign || n1->n_sign == PLUS)
 254             return (1);
 255           else
 256             return (-1);
 257         }
 258       else
 259         {
 260           /* Magnitude of n1 < n2. */
 261           if (!use_sign || n1->n_sign == PLUS)
 262             return (-1);
 263           else
 264             return (1);
 265         }
 266     }
 267
 268   /* They are equal up to the last part of the equal part of the fraction. */
 269   if (n1->n_scale != n2->n_scale)
 270     if (n1->n_scale > n2->n_scale)
 271       {
 272         for (count = n1->n_scale-n2->n_scale; count>0; count--)
 273           if (*n1ptr++ != 0)
 274             {
 275               /* Magnitude of n1 > n2. */
 276               if (!use_sign || n1->n_sign == PLUS)
 277                 return (1);
 278               else
 279                 return (-1);
 280             }
 281       }
 282     else
 283       {
 284         for (count = n2->n_scale-n1->n_scale; count>0; count--)
 285           if (*n2ptr++ != 0)
 286             {
 287               /* Magnitude of n1 < n2. */
 288               if (!use_sign || n1->n_sign == PLUS)
 289                 return (-1);
 290               else
 291                 return (1);
 292             }
 293       }
 294
 295   /* They must be equal! */
 296   return (0);
 297 }
 298
 299
 300 /* This is the "user callable" routine to compare numbers N1 and N2. */
 301
 302 int
 303 bc_compare (n1, n2)
 304      bc_num n1, n2;
 305 {
 306   return _do_compare (n1, n2, TRUE, FALSE);
 307 }
 308
 309
 310 /* In some places we need to check if the number NUM is zero. */
 311
 312 char
 313 is_zero (num)
 314      bc_num num;
 315 {
 316   int  count;
 317   char *nptr;
 318
 319   /* Quick check. */
 320   if (num == _zero_) return TRUE;
 321
 322   /* Initialize */
 323   count = num->n_len + num->n_scale;
 324   nptr = num->n_value;
 325
 326   /* The check */
 327   while ((count > 0) && (*nptr++ == 0)) count--;
 328
 329   if (count != 0)
 330     return FALSE;
 331   else
 332     return TRUE;
 333 }
 334
 335
 336 /* In some places we need to check if the number is negative. */
 337
 338 char
 339 is_neg (num)
 340      bc_num num;
 341 {
 342   return num->n_sign == MINUS;
 343 }
 344
 345
 346 /* For many things, we may have leading zeros in a number NUM.
 347    _rm_leading_zeros just moves the data to the correct
 348    place and adjusts the length. */
 349
 350 static void
 351 _rm_leading_zeros (num)
 352      bc_num num;
 353 {
 354   int bytes;
 355   char *dst, *src;
 356
 357   /* Do a quick check to see if we need to do it. */
 358   if (*num->n_value != 0) return;
 359
 360   /* The first digit is 0, find the first non-zero digit in the 10's or
 361      greater place. */
 362   bytes = num->n_len;
 363   src = num->n_value;
 364   while (bytes > 1 && *src == 0) src++, bytes--;
 365   num->n_len = bytes;
 366   bytes += num->n_scale;
 367   dst = num->n_value;
 368   while (bytes-- > 0) *dst++ = *src++;
 369
 370 }
 371
 372
 373 /* Perform addition: N1 is added to N2 and the value is
 374    returned.  The signs of N1 and N2 are ignored. */
 375
 376 static bc_num
 377 _do_add (n1, n2)
 378      bc_num n1, n2;
 379 {
 380   bc_num sum;
 381   int sum_scale, sum_digits;
 382   char *n1ptr, *n2ptr, *sumptr;
 383   int carry, n1bytes, n2bytes;
 384
 385   /* Prepare sum. */
 386   sum_scale = MAX (n1->n_scale, n2->n_scale);
 387   sum_digits = MAX (n1->n_len, n2->n_len) + 1;
 388   sum = new_num (sum_digits,sum_scale);
 389
 390   /* Start with the fraction part.  Initialize the pointers. */
 391   n1bytes = n1->n_scale;
 392   n2bytes = n2->n_scale;
 393   n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1);
 394   n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1);
 395   sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1);
 396
 397   /* Add the fraction part.  First copy the longer fraction.*/
 398   if (n1bytes != n2bytes)
 399     {
 400       if (n1bytes > n2bytes)
 401         while (n1bytes>n2bytes)
 402           { *sumptr-- = *n1ptr--; n1bytes--;}
 403       else
 404         while (n2bytes>n1bytes)
 405           { *sumptr-- = *n2ptr--; n2bytes--;}
 406     }
 407
 408   /* Now add the remaining fraction part and equal size integer parts. */
 409   n1bytes += n1->n_len;
 410   n2bytes += n2->n_len;
 411   carry = 0;
 412   while ((n1bytes > 0) && (n2bytes > 0))
 413     {
 414       *sumptr = *n1ptr-- + *n2ptr-- + carry;
 415       if (*sumptr > 9)
 416         {
 417            carry = 1;
 418            *sumptr -= 10;
 419         }
 420       else
 421         carry = 0;
 422       sumptr--;
 423       n1bytes--;
 424       n2bytes--;
 425     }
 426
 427   /* Now add carry the longer integer part. */
 428   if (n1bytes == 0)
 429     { n1bytes = n2bytes; n1ptr = n2ptr; }
 430   while (n1bytes-- > 0)
 431     {
 432       *sumptr = *n1ptr-- + carry;
 433       if (*sumptr > 9)
 434         {
 435            carry = 1;
 436            *sumptr -= 10;
 437          }
 438       else
 439         carry = 0;
 440       sumptr--;
 441     }
 442
 443   /* Set final carry. */
 444   if (carry == 1)
 445     *sumptr += 1;
 446
 447   /* Adjust sum and return. */
 448   _rm_leading_zeros (sum);
 449   return sum;
 450 }
 451
 452
 453 /* Perform subtraction: N2 is subtracted from N1 and the value is
 454    returned.  The signs of N1 and N2 are ignored.  Also, N1 is
 455    assumed to be larger than N2.  */
 456
 457 static bc_num
 458 _do_sub (n1, n2)
 459      bc_num n1, n2;
 460 {
 461   bc_num diff;
 462   int diff_scale, diff_len;
 463   int min_scale, min_len;
 464   char *n1ptr, *n2ptr, *diffptr;
 465   int borrow, count, val;
 466
 467   /* Allocate temporary storage. */
 468   diff_len = MAX (n1->n_len, n2->n_len);
 469   diff_scale = MAX (n1->n_scale, n2->n_scale);
 470   min_len = MIN  (n1->n_len, n2->n_len);
 471   min_scale = MIN (n1->n_scale, n2->n_scale);
 472   diff = new_num (diff_len, diff_scale);
 473
 474   /* Initialize the subtract. */
 475   n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1);
 476   n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1);
 477   diffptr = (char *) (diff->n_value + diff_len + diff_scale -1);
 478
 479   /* Subtract the numbers. */
 480   borrow = 0;
 481
 482   /* Take care of the longer scaled number. */
 483   if (n1->n_scale != min_scale)
 484     {
 485       /* n1 has the longer scale */
 486       for (count = n1->n_scale - min_scale; count > 0; count--)
 487         *diffptr-- = *n1ptr--;
 488     }
 489   else
 490     {
 491       /* n2 has the longer scale */
 492       for (count = n2->n_scale - min_scale; count > 0; count--)
 493         {
 494           val = - *n2ptr-- - borrow;
 495           if (val < 0)
 496             {
 497               val += 10;
 498               borrow = 1;
 499             }
 500           else
 501             borrow = 0;
 502           *diffptr-- = val;
 503         }
 504     }
 505
 506   /* Now do the equal length scale and integer parts. */
 507
 508   for (count = 0; count < min_len + min_scale; count++)
 509     {
 510       val = *n1ptr-- - *n2ptr-- - borrow;
 511       if (val < 0)
 512         {
 513           val += 10;
 514           borrow = 1;
 515         }
 516       else
 517         borrow = 0;
 518       *diffptr-- = val;
 519     }
 520
 521   /* If n1 has more digits then n2, we now do that subtract. */
 522   if (diff_len != min_len)
 523     {
 524       for (count = diff_len - min_len; count > 0; count--)
 525         {
 526           val = *n1ptr-- - borrow;
 527           if (val < 0)
 528             {
 529               val += 10;
 530               borrow = 1;
 531             }
 532           else
 533             borrow = 0;
 534           *diffptr-- = val;
 535         }
 536     }
 537
 538   /* Clean up and return. */
 539   _rm_leading_zeros (diff);
 540   return diff;
 541 }
 542
 543
 544 /* Here is the full add routine that takes care of negative numbers.
 545    N1 is added to N2 and the result placed into RESULT. */
 546
 547 void
 548 bc_add ( n1, n2, result)
 549      bc_num n1, n2, *result;
 550 {
 551   bc_num sum;
 552   int cmp_res;
 553
 554   if (n1->n_sign == n2->n_sign)
 555     {
 556       sum = _do_add (n1, n2);
 557       sum->n_sign = n1->n_sign;
 558     }
 559   else
 560     {
 561       /* subtraction must be done. */
 562       cmp_res = _do_compare (n1, n2, FALSE, FALSE);  /* Compare magnitudes. */
 563       switch (cmp_res)
 564         {
 565         case -1:
 566           /* n1 is less than n2, subtract n1 from n2. */
 567           sum = _do_sub (n2, n1);
 568           sum->n_sign = n2->n_sign;
 569           break;
 570         case  0:
 571           /* They are equal! return zero! */
 572           sum = copy_num (_zero_);
 573           break;
 574         case  1:
 575           /* n2 is less than n1, subtract n2 from n1. */
 576           sum = _do_sub (n1, n2);
 577           sum->n_sign = n1->n_sign;
 578         }
 579     }
 580
 581   /* Clean up and return. */
 582   free_num (result);
 583   *result = sum;
 584 }
 585
 586
 587 /* Here is the full subtract routine that takes care of negative numbers.
 588    N2 is subtracted from N1 and the result placed in RESULT. */
 589
 590 void
 591 bc_sub ( n1, n2, result)
 592      bc_num n1, n2, *result;
 593 {
 594   bc_num diff;
 595   int cmp_res;
 596
 597   if (n1->n_sign != n2->n_sign)
 598     {
 599       diff = _do_add (n1, n2);
 600       diff->n_sign = n1->n_sign;
 601     }
 602   else
 603     {
 604       /* subtraction must be done. */
 605       cmp_res = _do_compare (n1, n2, FALSE, FALSE);  /* Compare magnitudes. */
 606       switch (cmp_res)
 607         {
 608         case -1:
 609           /* n1 is less than n2, subtract n1 from n2. */
 610           diff = _do_sub (n2, n1);
 611           diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS);
 612           break;
 613         case  0:
 614           /* They are equal! return zero! */
 615           diff = copy_num (_zero_);
 616           break;
 617         case  1:
 618           /* n2 is less than n1, subtract n2 from n1. */
 619           diff = _do_sub (n1, n2);
 620           diff->n_sign = n1->n_sign;
 621           break;
 622         }
 623     }
 624
 625   /* Clean up and return. */
 626   free_num (result);
 627   *result = diff;
 628 }
 629
 630
 631 /* The multiply routine.  N2 time N1 is put int PROD with the scale of
 632    the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
 633    */
 634
 635 void
 636 bc_multiply (n1, n2, prod, scale)
 637      bc_num n1, n2, *prod;
 638      int scale;
 639 {
 640   bc_num pval;                  /* For the working storage. */
 641   char *n1ptr, *n2ptr, *pvptr;  /* Work pointers. */
 642   char *n1end, *n2end;          /* To the end of n1 and n2. */
 643
 644   int indx;
 645   int len1, len2, total_digits;
 646   long sum;
 647   int full_scale, prod_scale;
 648   int toss;
 649
 650   /* Initialize things. */
 651   len1 = n1->n_len + n1->n_scale;
 652   len2 = n2->n_len + n2->n_scale;
 653   total_digits = len1 + len2;
 654   full_scale = n1->n_scale + n2->n_scale;
 655   prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
 656   toss = full_scale - prod_scale;
 657   pval =  new_num (total_digits-full_scale, prod_scale);
 658   pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
 659   n1end = (char *) (n1->n_value + len1 - 1);
 660   n2end = (char *) (n2->n_value + len2 - 1);
 661   pvptr = (char *) (pval->n_value + total_digits - toss - 1);
 662   sum = 0;
 663
 664   /* Here are the loops... */
 665   for (indx = 0; indx < toss; indx++)
 666     {
 667       n1ptr = (char *) (n1end - MAX(0, indx-len2+1));
 668       n2ptr = (char *) (n2end - MIN(indx, len2-1));
 669       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
 670         sum += *n1ptr-- * *n2ptr++;
 671       sum = sum / 10;
 672     }
 673   for ( ; indx < total_digits-1; indx++)
 674     {
 675       n1ptr = (char *) (n1end - MAX(0, indx-len2+1));
 676       n2ptr = (char *) (n2end - MIN(indx, len2-1));
 677       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
 678         sum += *n1ptr-- * *n2ptr++;
 679       *pvptr-- = sum % 10;
 680       sum = sum / 10;
 681     }
 682   *pvptr-- = sum;
 683
 684   /* Assign to prod and clean up the number. */
 685   free_num (prod);
 686   *prod = pval;
 687   _rm_leading_zeros (*prod);
 688   if (is_zero (*prod))
 689     (*prod)->n_sign = PLUS;
 690 }
 691
 692
 693 /* Some utility routines for the divide:  First a one digit multiply.
 694    NUM (with SIZE digits) is multiplied by DIGIT and the result is
 695    placed into RESULT.  It is written so that NUM and RESULT can be
 696    the same pointers.  */
 697
 698 static void
 699 _one_mult (num, size, digit, result)
 700      unsigned char *num;
 701      int size, digit;
 702      unsigned char *result;
 703 {
 704   int carry, value;
 705   unsigned char *nptr, *rptr;
 706
 707   if (digit == 0)
 708     memset (result, 0, size);
 709   else
 710     {
 711       if (digit == 1)
 712         memcpy (result, num, size);
 713       else
 714         {
 715           /* Initialize */
 716           nptr = (unsigned char *) (num+size-1);
 717           rptr = (unsigned char *) (result+size-1);
 718           carry = 0;
 719
 720           while (size-- > 0)
 721             {
 722               value = *nptr-- * digit + carry;
 723               *rptr-- = value % 10;
 724               carry = value / 10;
 725             }
 726
 727           if (carry != 0) *rptr = carry;
 728         }
 729     }
 730 }
 731
 732
 733 /* The full division routine. This computes N1 / N2.  It returns
 734    0 if the division is ok and the result is in QUOT.  The number of
 735    digits after the decimal point is SCALE. It returns -1 if division
 736    by zero is tried.  The algorithm is found in Knuth Vol 2. p237. */
 737
 738 int
 739 bc_divide (n1, n2, quot, scale)
 740      bc_num n1, n2, *quot;
 741      int scale;
 742 {
 743   bc_num qval;
 744   unsigned char *num1, *num2;
 745   unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
 746   int  scale1, val;
 747   unsigned int  len1, len2, scale2, qdigits, extra, count;
 748   unsigned int  qdig, qguess, borrow, carry;
 749   unsigned char *mval;
 750   char zero;
 751   unsigned int  norm;
 752
 753   /* Test for divide by zero. */
 754   if (is_zero (n2)) return -1;
 755
 756   /* Test for divide by 1.  If it is we must truncate. */
 757   if (n2->n_scale == 0)
 758     {
 759       if (n2->n_len == 1 && *n2->n_value == 1)
 760         {
 761           qval = new_num (n1->n_len, scale);
 762           qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
 763           memset (&qval->n_value[n1->n_len],0,scale);
 764           memcpy (qval->n_value, n1->n_value,
 765                   n1->n_len + MIN(n1->n_scale,scale));
 766           free_num (quot);
 767           *quot = qval;
 768         }
 769     }
 770
 771   /* Set up the divide.  Move the decimal point on n1 by n2's scale.
 772      Remember, zeros on the end of num2 are wasted effort for dividing. */
 773   scale2 = n2->n_scale;
 774   n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
 775   while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
 776
 777   len1 = n1->n_len + scale2;
 778   scale1 = n1->n_scale - scale2;
 779   if (scale1 < scale)
 780     extra = scale - scale1;
 781   else
 782     extra = 0;
 783   num1 = (unsigned char *) malloc (n1->n_len+n1->n_scale+extra+2);
 784   if (num1 == NULL) out_of_memory();
 785   memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
 786   memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
 787
 788   len2 = n2->n_len + scale2;
 789   num2 = (unsigned char *) malloc (len2+1);
 790   if (num2 == NULL) out_of_memory();
 791   memcpy (num2, n2->n_value, len2);
 792   *(num2+len2) = 0;
 793   n2ptr = num2;
 794   while (*n2ptr == 0)
 795     {
 796       n2ptr++;
 797       len2--;
 798     }
 799
 800   /* Calculate the number of quotient digits. */
 801   if (len2 > len1+scale)
 802     {
 803       qdigits = scale+1;
 804       zero = TRUE;
 805     }
 806   else
 807     {
 808       zero = FALSE;
 809       if (len2>len1)
 810         qdigits = scale+1;      /* One for the zero integer part. */
 811       else
 812         qdigits = len1-len2+scale+1;
 813     }
 814
 815   /* Allocate and zero the storage for the quotient. */
 816   qval = new_num (qdigits-scale,scale);
 817   memset (qval->n_value, 0, qdigits);
 818
 819   /* Allocate storage for the temporary storage mval. */
 820   mval = (unsigned char *) malloc (len2+1);
 821   if (mval == NULL) out_of_memory ();
 822
 823   /* Now for the full divide algorithm. */
 824   if (!zero)
 825     {
 826       /* Normalize */
 827       norm =  10 / ((int)*n2ptr + 1);
 828       if (norm != 1)
 829         {
 830           _one_mult (num1, len1+scale1+extra+1, norm, num1);
 831           _one_mult (n2ptr, len2, norm, n2ptr);
 832         }
 833
 834       /* Initialize divide loop. */
 835       qdig = 0;
 836       if (len2 > len1)
 837         qptr = (unsigned char *) qval->n_value+len2-len1;
 838       else
 839         qptr = (unsigned char *) qval->n_value;
 840
 841       /* Loop */
 842       while (qdig <= len1+scale-len2)
 843         {
 844           /* Calculate the quotient digit guess. */
 845           if (*n2ptr == num1[qdig])
 846             qguess = 9;
 847           else
 848             qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
 849
 850           /* Test qguess. */
 851           if (n2ptr[1]*qguess >
 852               (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
 853                + num1[qdig+2])
 854             {
 855               qguess--;
 856               /* And again. */
 857               if (n2ptr[1]*qguess >
 858                   (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
 859                   + num1[qdig+2])
 860                 qguess--;
 861             }
 862
 863           /* Multiply and subtract. */
 864           borrow = 0;
 865           if (qguess != 0)
 866             {
 867               *mval = 0;
 868               _one_mult (n2ptr, len2, qguess, mval+1);
 869               ptr1 = (unsigned char *) num1+qdig+len2;
 870               ptr2 = (unsigned char *) mval+len2;
 871               for (count = 0; count < len2+1; count++)
 872                 {
 873                   val = (int) *ptr1 - (int) *ptr2-- - borrow;
 874                   if (val < 0)
 875                     {
 876                       val += 10;
 877                       borrow = 1;
 878                     }
 879                   else
 880                     borrow = 0;
 881                   *ptr1-- = val;
 882                 }
 883             }
 884
 885           /* Test for negative result. */
 886           if (borrow == 1)
 887             {
 888               qguess--;
 889               ptr1 = (unsigned char *) num1+qdig+len2;
 890               ptr2 = (unsigned char *) n2ptr+len2-1;
 891               carry = 0;
 892               for (count = 0; count < len2; count++)
 893                 {
 894                   val = (int) *ptr1 + (int) *ptr2-- + carry;
 895                   if (val > 9)
 896                     {
 897                       val -= 10;
 898                       carry = 1;
 899                     }
 900                   else
 901                     carry = 0;
 902                   *ptr1-- = val;
 903                 }
 904               if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
 905             }
 906
 907           /* We now know the quotient digit. */
 908           *qptr++ =  qguess;
 909           qdig++;
 910         }
 911     }
 912
 913   /* Clean up and return the number. */
 914   qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
 915   if (is_zero (qval)) qval->n_sign = PLUS;
 916   _rm_leading_zeros (qval);
 917   free_num (quot);
 918   *quot = qval;
 919
 920   /* Clean up temporary storage. */
 921   free (mval);
 922   free (num1);
 923   free (num2);
 924
 925   return 0;     /* Everything is OK. */
 926 }
 927
 928
 929 /* Modulo for numbers.  This computes NUM1 % NUM2  and puts the
 930    result in RESULT.   */
 931
 932 int
 933 bc_modulo (num1, num2, result, scale)
 934      bc_num num1, num2, *result;
 935      int scale;
 936 {
 937   bc_num temp;
 938   int rscale;
 939
 940   /* Check for correct numbers. */
 941   if (is_zero (num2)) return -1;
 942
 943   /* Calculate final scale. */
 944   rscale = MAX (num1->n_scale, num2->n_scale+scale);
 945   init_num (&temp);
 946
 947   /* Calculate it. */
 948   bc_divide (num1, num2, &temp, scale);
 949   bc_multiply (temp, num2, &temp, rscale);
 950   bc_sub (num1, temp, result);
 951   free_num (&temp);
 952
 953   return 0;     /* Everything is OK. */
 954 }
 955
 956
 957 /* Raise NUM1 to the NUM2 power.  The result is placed in RESULT.
 958    Maximum exponent is LONG_MAX.  If a NUM2 is not an integer,
 959    only the integer part is used.  */
 960
 961 void
 962 bc_raise (num1, num2, result, scale)
 963      bc_num num1, num2, *result;
 964      int scale;
 965 {
 966    bc_num temp, power;
 967    long exponent;
 968    int rscale;
 969    char neg;
 970
 971    /* Check the exponent for scale digits and convert to a long. */
 972    if (num2->n_scale != 0)
 973      rt_warn ("non-zero scale in exponent");
 974    exponent = num2long (num2);
 975    if (exponent == 0 && (num2->n_len > 1 || num2->n_value[0] != 0))
 976        rt_error ("exponent too large in raise");
 977
 978    /* Special case if exponent is a zero. */
 979    if (exponent == 0)
 980      {
 981        free_num (result);
 982        *result = copy_num (_one_);
 983        return;
 984      }
 985
 986    /* Other initializations. */
 987    if (exponent < 0)
 988      {
 989        neg = TRUE;
 990        exponent = -exponent;
 991        rscale = scale;
 992      }
 993    else
 994      {
 995        neg = FALSE;
 996        rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale));
 997      }
 998    temp = copy_num (_one_);
 999    power = copy_num (num1);
1000
1001    /* Do the calculation. */
1002    while (exponent != 0)
1003      {
1004        if (exponent & 1 != 0)
1005          bc_multiply (temp, power, &temp, rscale);
1006        bc_multiply (power, power, &power, rscale);
1007        exponent = exponent >> 1;
1008      }
1009
1010    /* Assign the value. */
1011    if (neg)
1012      {
1013        bc_divide (_one_, temp, result, rscale);
1014        free_num (&temp);
1015      }
1016    else
1017      {
1018        free_num (result);
1019        *result = temp;
1020      }
1021    free_num (&power);
1022 }
1023
1024
1025 /* Take the square root NUM and return it in NUM with SCALE digits
1026    after the decimal place. */
1027
1028 int
1029 bc_sqrt (num, scale)
1030      bc_num *num;
1031      int scale;
1032 {
1033   int rscale, cmp_res, done;
1034   int cscale;
1035   bc_num guess, guess1, point5;
1036
1037   /* Initial checks. */
1038   cmp_res = bc_compare (*num, _zero_);
1039   if (cmp_res < 0)
1040     return 0;           /* error */
1041   else
1042     {
1043       if (cmp_res == 0)
1044         {
1045           free_num (num);
1046           *num = copy_num (_zero_);
1047           return 1;
1048         }
1049     }
1050   cmp_res = bc_compare (*num, _one_);
1051   if (cmp_res == 0)
1052     {
1053       free_num (num);
1054       *num = copy_num (_one_);
1055       return 1;
1056     }
1057
1058   /* Initialize the variables. */
1059   rscale = MAX (scale, (*num)->n_scale);
1060   cscale = rscale + 2;
1061   init_num (&guess);
1062   init_num (&guess1);
1063   point5 = new_num (1,1);
1064   point5->n_value[1] = 5;
1065
1066
1067   /* Calculate the initial guess. */
1068   if (cmp_res < 0)
1069     /* The number is between 0 and 1.  Guess should start at 1. */
1070     guess = copy_num (_one_);
1071   else
1072     {
1073       /* The number is greater than 1.  Guess should start at 10^(exp/2). */
1074       int2num (&guess,10);
1075       int2num (&guess1,(*num)->n_len);
1076       bc_multiply (guess1, point5, &guess1, rscale);
1077       guess1->n_scale = 0;
1078       bc_raise (guess, guess1, &guess, rscale);
1079       free_num (&guess1);
1080     }
1081
1082   /* Find the square root using Newton's algorithm. */
1083   done = FALSE;
1084   while (!done)
1085     {
1086       free_num (&guess1);
1087       guess1 = copy_num (guess);
1088       bc_divide (*num,guess,&guess,cscale);
1089       bc_add (guess,guess1,&guess);
1090       bc_multiply (guess,point5,&guess,cscale);
1091       cmp_res = _do_compare (guess,guess1,FALSE,TRUE);
1092       if (cmp_res == 0) done = TRUE;
1093     }
1094
1095   /* Assign the number and clean up. */
1096   free_num (num);
1097   bc_divide (guess,_one_,num,rscale);
1098   free_num (&guess);
1099   free_num (&guess1);
1100   free_num (&point5);
1101   return 1;
1102 }
1103
1104
1105 /* The following routines provide output for bcd numbers package
1106    using the rules of POSIX bc for output. */
1107
1108 /* This structure is used for saving digits in the conversion process. */
1109 typedef struct stk_rec {
1110         long  digit;
1111         struct stk_rec *next;
1112 } stk_rec;
1113
1114 /* The reference string for digits. */
1115 char ref_str[] = "0123456789ABCDEF";
1116
1117
1118 /* A special output routine for "multi-character digits."  Exactly
1119    SIZE characters must be output for the value VAL.  If SPACE is
1120    non-zero, we must output one space before the number.  OUT_CHAR
1121    is the actual routine for writing the characters. */
1122
1123 void
1124 out_long (val, size, space, out_char)
1125      long val;
1126      int size, space;
1127 #ifdef __STDC__
1128      void (*out_char)(int);
1129 #else
1130      void (*out_char)();
1131 #endif
1132 {
1133   char digits[40];
1134   int len, ix;
1135
1136   if (space) (*out_char) (' ');
1137   sprintf (digits, "%ld", val);
1138   len = strlen (digits);
1139   while (size > len)
1140     {
1141       (*out_char) ('0');
1142       size--;
1143     }
1144   for (ix=0; ix < len; ix++)
1145     (*out_char) (digits[ix]);
1146 }
1147
1148 /* Output of a bcd number.  NUM is written in base O_BASE using OUT_CHAR
1149    as the routine to do the actual output of the characters. */
1150
1151 void
1152 out_num (num, o_base, out_char)
1153      bc_num num;
1154      int o_base;
1155 #ifdef __STDC__
1156      void (*out_char)(int);
1157 #else
1158      void (*out_char)();
1159 #endif
1160 {
1161   char *nptr;
1162   int  index, fdigit, pre_space;
1163   stk_rec *digits, *temp;
1164   bc_num int_part, frac_part, base, cur_dig, t_num, max_o_digit;
1165
1166   /* The negative sign if needed. */
1167   if (num->n_sign == MINUS) (*out_char) ('-');
1168
1169   /* Output the number. */
1170   if (is_zero (num))
1171     (*out_char) ('0');
1172   else
1173     if (o_base == 10)
1174       {
1175         /* The number is in base 10, do it the fast way. */
1176         nptr = num->n_value;
1177         if (num->n_len > 1 || *nptr != 0)
1178           for (index=num->n_len; index>0; index--)
1179             (*out_char) (BCD_CHAR(*nptr++));
1180         else
1181           nptr++;
1182
1183         /* Now the fraction. */
1184         if (num->n_scale > 0)
1185           {
1186             (*out_char) ('.');
1187             for (index=0; index<num->n_scale; index++)
1188               (*out_char) (BCD_CHAR(*nptr++));
1189           }
1190       }
1191     else
1192       {
1193         /* The number is some other base. */
1194         digits = NULL;
1195         init_num (&int_part);
1196         bc_divide (num, _one_, &int_part, 0);
1197         init_num (&frac_part);
1198         init_num (&cur_dig);
1199         init_num (&base);
1200         bc_sub (num, int_part, &frac_part);
1201         int2num (&base, o_base);
1202         init_num (&max_o_digit);
1203         int2num (&max_o_digit, o_base-1);
1204
1205
1206         /* Get the digits of the integer part and push them on a stack. */
1207         while (!is_zero (int_part))
1208           {
1209             bc_modulo (int_part, base, &cur_dig, 0);
1210             temp = (stk_rec *) malloc (sizeof(stk_rec));
1211             if (temp == NULL) out_of_memory();
1212             temp->digit = num2long (cur_dig);
1213             temp->next = digits;
1214             digits = temp;
1215             bc_divide (int_part, base, &int_part, 0);
1216           }
1217
1218         /* Print the digits on the stack. */
1219         if (digits != NULL)
1220           {
1221             /* Output the digits. */
1222             while (digits != NULL)
1223               {
1224                 temp = digits;
1225                 digits = digits->next;
1226                 if (o_base <= 16)
1227                   (*out_char) (ref_str[ (int) temp->digit]);
1228                 else
1229                   out_long (temp->digit, max_o_digit->n_len, 1, out_char);
1230                 free (temp);
1231               }
1232           }
1233
1234         /* Get and print the digits of the fraction part. */
1235         if (num->n_scale > 0)
1236           {
1237             (*out_char) ('.');
1238             pre_space = 0;
1239             t_num = copy_num (_one_);
1240             while (t_num->n_len <= num->n_scale) {
1241               bc_multiply (frac_part, base, &frac_part, num->n_scale);
1242               fdigit = num2long (frac_part);
1243               int2num (&int_part, fdigit);
1244               bc_sub (frac_part, int_part, &frac_part);
1245               if (o_base <= 16)
1246                 (*out_char) (ref_str[fdigit]);
1247               else {
1248                 out_long (fdigit, max_o_digit->n_len, pre_space, out_char);
1249                 pre_space = 1;
1250               }
1251               bc_multiply (t_num, base, &t_num, 0);
1252             }
1253           }
1254
1255         /* Clean up. */
1256         free_num (&int_part);
1257         free_num (&frac_part);
1258         free_num (&base);
1259         free_num (&cur_dig);
1260       }
1261 }
1262
1263
1264 #if DEBUG > 0
1265
1266 /* Debugging procedures.  Some are just so one can call them from the
1267    debugger.  */
1268
1269 /* p_n prints the number NUM in base 10. */
1270
1271 void
1272 p_n (num)
1273      bc_num num;
1274 {
1275   out_num (num, 10, out_char);
1276   return 0;
1277 }
1278
1279
1280 /* p_b prints a character array as if it was a string of bcd digits. */
1281 void
1282 p_v (name, num, len)
1283      char *name;
1284      unsigned char *num;
1285      int len;
1286 {
1287   int i;
1288   printf ("%s=", name);
1289   for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i]));
1290   printf ("\n");
1291 }
1292
1293
1294 /* Convert strings to bc numbers.  Base 10 only.*/
1295
1296 void
1297 str2num (num, str, scale)
1298      bc_num *num;
1299      char *str;
1300      int scale;
1301 {
1302   int digits, strscale;
1303   char *ptr, *nptr;
1304   char zero_int;
1305
1306   /* Prepare num. */
1307   free_num (num);
1308
1309   /* Check for valid number and count digits. */
1310   ptr = str;
1311   digits = 0;
1312   strscale = 0;
1313   zero_int = FALSE;
1314   if ( (*ptr == '+') || (*ptr == '-'))  ptr++;  /* Sign */
1315   while (*ptr == '0') ptr++;                    /* Skip leading zeros. */
1316   while (isdigit(*ptr)) ptr++, digits++;        /* digits */
1317   if (*ptr == '.') ptr++;                       /* decimal point */
1318   while (isdigit(*ptr)) ptr++, strscale++;      /* digits */
1319   if ((*ptr != '\0') || (digits+strscale == 0))
1320     {
1321       *num = copy_num (_zero_);
1322       return;
1323     }
1324
1325   /* Adjust numbers and allocate storage and initialize fields. */
1326   strscale = MIN(strscale, scale);
1327   if (digits == 0)
1328     {
1329       zero_int = TRUE;
1330       digits = 1;
1331     }
1332   *num = new_num (digits, strscale);
1333
1334   /* Build the whole number. */
1335   ptr = str;
1336   if (*ptr == '-')
1337     {
1338       (*num)->n_sign = MINUS;
1339       ptr++;
1340     }
1341   else
1342     {
1343       (*num)->n_sign = PLUS;
1344       if (*ptr == '+') ptr++;
1345     }
1346   while (*ptr == '0') ptr++;                    /* Skip leading zeros. */
1347   nptr = (*num)->n_value;
1348   if (zero_int)
1349     {
1350       *nptr++ = 0;
1351       digits = 0;
1352     }
1353   for (;digits > 0; digits--)
1354     *nptr++ = CH_VAL(*ptr++);
1355
1356
1357   /* Build the fractional part. */
1358   if (strscale > 0)
1359     {
1360       ptr++;  /* skip the decimal point! */
1361       for (;strscale > 0; strscale--)
1362         *nptr++ = CH_VAL(*ptr++);
1363     }
1364 }
1365
1366 /* Convert a numbers to a string.  Base 10 only.*/
1367
1368 char
1369 *num2str (num)
1370       bc_num num;
1371 {
1372   char *str, *sptr;
1373   char *nptr;
1374   int  index, signch;
1375
1376   /* Allocate the string memory. */
1377   signch = ( num->n_sign == PLUS ? 0 : 1 );  /* Number of sign chars. */
1378   if (num->n_scale > 0)
1379     str = (char *) malloc (num->n_len + num->n_scale + 2 + signch);
1380   else
1381     str = (char *) malloc (num->n_len + 1 + signch);
1382   if (str == NULL) out_of_memory();
1383
1384   /* The negative sign if needed. */
1385   sptr = str;
1386   if (signch) *sptr++ = '-';
1387
1388   /* Load the whole number. */
1389   nptr = num->n_value;
1390   for (index=num->n_len; index>0; index--)
1391     *sptr++ = BCD_CHAR(*nptr++);
1392
1393   /* Now the fraction. */
1394   if (num->n_scale > 0)
1395     {
1396       *sptr++ = '.';
1397       for (index=0; index<num->n_scale; index++)
1398         *sptr++ = BCD_CHAR(*nptr++);
1399     }
1400
1401   /* Terminate the string and return it! */
1402   *sptr = '\0';
1403   return (str);
1404 }
1405 #endif