/* number.c: Implements arbitrary precision numbers. */

/*

    Copyright (C) 1991, 1992, 1993, 1994, 1997, 2000, 2012-2017 Free Software Foundation, Inc.

    This program is free software; you can redistribute it and/or modify

    it under the terms of the GNU General Public License as published by

    the Free Software Foundation; either version 3 of the License , or

    (at your option) any later version.

    This program is distributed in the hope that it will be useful,

    but WITHOUT ANY WARRANTY; without even the implied warranty of

    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License

    along with this program; see the file COPYING.  If not, write to:

      The Free Software Foundation, Inc.

      51 Franklin Street, Fifth Floor

      Boston, MA 02110-1301  USA

    You may contact the author by:

       e-mail:  philnelson@acm.org

      us-mail:  Philip A. Nelson

                Computer Science Department, 9062

                Western Washington University

                Bellingham, WA 98226-9062

*************************************************************************/

#include <stdio.h>

#include <config.h>

#include <number.h>

#include <assert.h>

#ifdef HAVE_STDLIB_H

#include <stdlib.h>

#endif

#ifdef HAVE_STRING_H

#include <string.h>

#endif

#include <ctype.h>

/* Prototypes needed for external utility routines. */

#define bc_rt_warn rt_warn

#define bc_rt_error rt_error

#define bc_out_of_memory out_of_memory

void rt_warn (const char *mesg ,...);

void rt_error (const char *mesg ,...);

void out_of_memory (void);

/* Storage used for special numbers. */

bc_num _zero_;

bc_num _one_;

bc_num _two_;

static bc_num _bc_Free_list = NULL;

/* new_num allocates a number and sets fields to known values. */

bc_num

bc_new_num (int length, int scale)

{

  bc_num temp;

  if (_bc_Free_list != NULL) {

    temp = _bc_Free_list;

    _bc_Free_list = temp->n_next;

  } else {

    temp = (bc_num) malloc (sizeof(bc_struct));

    if (temp == NULL) bc_out_of_memory ();

  }

  temp->n_sign = PLUS;

  temp->n_len = length;

  temp->n_scale = scale;

  temp->n_refs = 1;

  temp->n_ptr = (char *) malloc (length+scale);

  if (temp->n_ptr == NULL) bc_out_of_memory();

  temp->n_value = temp->n_ptr;

  memset (temp->n_ptr, 0, length+scale);

  return temp;

}

/* "Frees" a bc_num NUM.  Actually decreases reference count and only

   frees the storage if reference count is zero. */

void

bc_free_num (bc_num *num)

{

  if (*num == NULL) return;

  (*num)->n_refs--;

  if ((*num)->n_refs == 0) {

    if ((*num)->n_ptr)

      free ((*num)->n_ptr);

    (*num)->n_next = _bc_Free_list;

    _bc_Free_list = *num;

  }

  *num = NULL;

}

/* Intitialize the number package! */

void

bc_init_numbers (void)

{

  _zero_ = bc_new_num (1,0);

  _one_  = bc_new_num (1,0);

  _one_->n_value[0] = 1;

  _two_  = bc_new_num (1,0);

  _two_->n_value[0] = 2;

}

/* Make a copy of a number!  Just increments the reference count! */

bc_num

bc_copy_num (bc_num num)

{

  num->n_refs++;

  return num;

}

/* Initialize a number NUM by making it a copy of zero. */

void

bc_init_num (bc_num *num)

{

  *num = bc_copy_num (_zero_);

}

/* For many things, we may have leading zeros in a number NUM.

   _bc_rm_leading_zeros just moves the data "value" pointer to the

   correct place and adjusts the length. */

static void

_bc_rm_leading_zeros (bc_num num)

{

  /* We can move n_value to point to the first non zero digit! */

  while (*num->n_value == 0 && num->n_len > 1) {

    num->n_value++;

    num->n_len--;

  }

}

/* Compare two bc numbers.  Return value is 0 if equal, -1 if N1 is less

   than N2 and +1 if N1 is greater than N2.  If USE_SIGN is false, just

   compare the magnitudes. */

static int

_bc_do_compare ( bc_num n1, bc_num n2, int use_sign, int ignore_last )

{

  char *n1ptr, *n2ptr;

  int  count;

  /* First, compare signs. */

  if (use_sign && n1->n_sign != n2->n_sign)

    {

      if (n1->n_sign == PLUS)

	return (1);	/* Positive N1 > Negative N2 */

      else

	return (-1);	/* Negative N1 < Positive N1 */

    }

  /* Now compare the magnitude. */

  if (n1->n_len != n2->n_len)

    {

      if (n1->n_len > n2->n_len)

	{

	  /* Magnitude of n1 > n2. */

	  if (!use_sign || n1->n_sign == PLUS)

	    return (1);

	  else

	    return (-1);

	}

      else

	{

	  /* Magnitude of n1 < n2. */

	  if (!use_sign || n1->n_sign == PLUS)

	    return (-1);

	  else

	    return (1);

	}

    }

  /* If we get here, they have the same number of integer digits.

     check the integer part and the equal length part of the fraction. */

  count = n1->n_len + MIN (n1->n_scale, n2->n_scale);

  n1ptr = n1->n_value;

  n2ptr = n2->n_value;

  while ((count > 0) && (*n1ptr == *n2ptr))

    {

      n1ptr++;

      n2ptr++;

      count--;

    }

  if (ignore_last && count == 1 && n1->n_scale == n2->n_scale)

    return (0);

  if (count != 0)

    {

      if (*n1ptr > *n2ptr)

	{

	  /* Magnitude of n1 > n2. */

	  if (!use_sign || n1->n_sign == PLUS)

	    return (1);

	  else

	    return (-1);

	}

      else

	{

	  /* Magnitude of n1 < n2. */

	  if (!use_sign || n1->n_sign == PLUS)

	    return (-1);

	  else

	    return (1);

	}

    }

  /* They are equal up to the last part of the equal part of the fraction. */

  if (n1->n_scale != n2->n_scale)

    {

      if (n1->n_scale > n2->n_scale)

	{

	  for (count = n1->n_scale-n2->n_scale; count>0; count--)

	    if (*n1ptr++ != 0)

	      {

		/* Magnitude of n1 > n2. */

		if (!use_sign || n1->n_sign == PLUS)

		  return (1);

		else

		  return (-1);

	      }

	}

      else

	{

	  for (count = n2->n_scale-n1->n_scale; count>0; count--)

	    if (*n2ptr++ != 0)

	      {

		/* Magnitude of n1 < n2. */

		if (!use_sign || n1->n_sign == PLUS)

		  return (-1);

		else

		  return (1);

	      }

	}

    }

  /* They must be equal! */

  return (0);

}

/* This is the "user callable" routine to compare numbers N1 and N2. */

int

bc_compare ( bc_num n1, bc_num n2 )

{

  return _bc_do_compare (n1, n2, TRUE, FALSE);

}

/* In some places we need to check if the number is negative. */

char

bc_is_neg (bc_num num)

{

  return num->n_sign == MINUS;

}

/* In some places we need to check if the number NUM is zero. */

char

bc_is_zero (bc_num num)

{

  int  count;

  char *nptr;

  /* Quick check. */

  if (num == _zero_) return TRUE;

  /* Initialize */

  count = num->n_len + num->n_scale;

  nptr = num->n_value;

  /* The check */

  while ((count > 0) && (*nptr++ == 0)) count--;

  if (count != 0)

    return FALSE;

  else

    return TRUE;

}

/* In some places we need to check if the number NUM is almost zero.

   Specifically, all but the last digit is 0 and the last digit is 1.

   Last digit is defined by scale. */

char

bc_is_near_zero (bc_num num, int scale)

{

  int  count;

  char *nptr;

  /* Error checking */

  if (scale > num->n_scale)

    scale = num->n_scale;

  /* Initialize */

  count = num->n_len + scale;

  nptr = num->n_value;

  /* The check */

  while ((count > 0) && (*nptr++ == 0)) count--;

  if (count != 0 && (count != 1 || *--nptr != 1))

    return FALSE;

  else

    return TRUE;

}

/* Perform addition: N1 is added to N2 and the value is

   returned.  The signs of N1 and N2 are ignored.

   SCALE_MIN is to set the minimum scale of the result. */

static bc_num

_bc_do_add (bc_num n1, bc_num n2, int scale_min)

{

  bc_num sum;

  int sum_scale, sum_digits;

  char *n1ptr, *n2ptr, *sumptr;

  int carry, n1bytes, n2bytes;

  int count;

  /* Prepare sum. */

  sum_scale = MAX (n1->n_scale, n2->n_scale);

  sum_digits = MAX (n1->n_len, n2->n_len) + 1;

  sum = bc_new_num (sum_digits, MAX(sum_scale, scale_min));

  /* Zero extra digits made by scale_min. */

  if (scale_min > sum_scale)

    {

      sumptr = (char *) (sum->n_value + sum_scale + sum_digits);

      for (count = scale_min - sum_scale; count > 0; count--)

	*sumptr++ = 0;

    }

  /* Start with the fraction part.  Initialize the pointers. */

  n1bytes = n1->n_scale;

  n2bytes = n2->n_scale;

  n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1);

  n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1);

  sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1);

  /* Add the fraction part.  First copy the longer fraction.*/

  if (n1bytes != n2bytes)

    {

      if (n1bytes > n2bytes)

	while (n1bytes>n2bytes)

	  { *sumptr-- = *n1ptr--; n1bytes--;}

      else

	while (n2bytes>n1bytes)

	  { *sumptr-- = *n2ptr--; n2bytes--;}

    }

  /* Now add the remaining fraction part and equal size integer parts. */

  n1bytes += n1->n_len;

  n2bytes += n2->n_len;

  carry = 0;

  while ((n1bytes > 0) && (n2bytes > 0))

    {

      *sumptr = *n1ptr-- + *n2ptr-- + carry;

      if (*sumptr > (BASE-1))

	{

	   carry = 1;

	   *sumptr -= BASE;

	}

      else

	carry = 0;

      sumptr--;

      n1bytes--;

      n2bytes--;

    }

  /* Now add carry the longer integer part. */

  if (n1bytes == 0)

    { n1bytes = n2bytes; n1ptr = n2ptr; }

  while (n1bytes-- > 0)

    {

      *sumptr = *n1ptr-- + carry;

      if (*sumptr > (BASE-1))

	{

	   carry = 1;

	   *sumptr -= BASE;

	 }

      else

	carry = 0;

      sumptr--;

    }

  /* Set final carry. */

  if (carry == 1)

    *sumptr += 1;

  /* Adjust sum and return. */

  _bc_rm_leading_zeros (sum);

  return sum;

}

/* Perform subtraction: N2 is subtracted from N1 and the value is

   returned.  The signs of N1 and N2 are ignored.  Also, N1 is

   assumed to be larger than N2.  SCALE_MIN is the minimum scale

   of the result. */

static bc_num

_bc_do_sub (bc_num n1, bc_num n2, int scale_min)

{

  bc_num diff;

  int diff_scale, diff_len;

  int min_scale, min_len;

  char *n1ptr, *n2ptr, *diffptr;

  int borrow, count, val;

  /* Allocate temporary storage. */

  diff_len = MAX (n1->n_len, n2->n_len);

  diff_scale = MAX (n1->n_scale, n2->n_scale);

  min_len = MIN  (n1->n_len, n2->n_len);

  min_scale = MIN (n1->n_scale, n2->n_scale);

  diff = bc_new_num (diff_len, MAX(diff_scale, scale_min));

  /* Zero extra digits made by scale_min. */

  if (scale_min > diff_scale)

    {

      diffptr = (char *) (diff->n_value + diff_len + diff_scale);

      for (count = scale_min - diff_scale; count > 0; count--)

	*diffptr++ = 0;

    }

  /* Initialize the subtract. */

  n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1);

  n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1);

  diffptr = (char *) (diff->n_value + diff_len + diff_scale -1);

  /* Subtract the numbers. */

  borrow = 0;

  /* Take care of the longer scaled number. */

  if (n1->n_scale != min_scale)

    {

      /* n1 has the longer scale */

      for (count = n1->n_scale - min_scale; count > 0; count--)

	*diffptr-- = *n1ptr--;

    }

  else

    {

      /* n2 has the longer scale */

      for (count = n2->n_scale - min_scale; count > 0; count--)

	{

	  val = - *n2ptr-- - borrow;

	  if (val < 0)

	    {

	      val += BASE;

	      borrow = 1;

	    }

	  else

	    borrow = 0;

	  *diffptr-- = val;

	}

    }

  /* Now do the equal length scale and integer parts. */

  for (count = 0; count < min_len + min_scale; count++)

    {

      val = *n1ptr-- - *n2ptr-- - borrow;

      if (val < 0)

	{

	  val += BASE;

	  borrow = 1;

	}

      else

	borrow = 0;

      *diffptr-- = val;

    }

  /* If n1 has more digits then n2, we now do that subtract. */

  if (diff_len != min_len)

    {

      for (count = diff_len - min_len; count > 0; count--)

	{

	  val = *n1ptr-- - borrow;

	  if (val < 0)

	    {

	      val += BASE;

	      borrow = 1;

	    }

	  else

	    borrow = 0;

	  *diffptr-- = val;

	}

    }

  /* Clean up and return. */

  _bc_rm_leading_zeros (diff);

  return diff;

}

/* Here is the full subtract routine that takes care of negative numbers.

   N2 is subtracted from N1 and the result placed in RESULT.  SCALE_MIN

   is the minimum scale for the result. */

void

bc_sub (bc_num n1, bc_num n2, bc_num *result, int scale_min)

{

  bc_num diff = NULL;

  int cmp_res;

  int res_scale;

  if (n1->n_sign != n2->n_sign)

    {

      diff = _bc_do_add (n1, n2, scale_min);

      diff->n_sign = n1->n_sign;

    }

  else

    {

      /* subtraction must be done. */

      /* Compare magnitudes. */

      cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);

      switch (cmp_res)

	{

	case -1:

	  /* n1 is less than n2, subtract n1 from n2. */

	  diff = _bc_do_sub (n2, n1, scale_min);

	  diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS);

	  break;

	case  0:

	  /* They are equal! return zero! */

	  res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));

	  diff = bc_new_num (1, res_scale);

	  memset (diff->n_value, 0, res_scale+1);

	  break;

	case  1:

	  /* n2 is less than n1, subtract n2 from n1. */

	  diff = _bc_do_sub (n1, n2, scale_min);

	  diff->n_sign = n1->n_sign;

	  break;

	}

    }

  /* Clean up and return. */

  bc_free_num (result);

  *result = diff;

}

/* Here is the full add routine that takes care of negative numbers.

   N1 is added to N2 and the result placed into RESULT.  SCALE_MIN

   is the minimum scale for the result. */

void

bc_add (bc_num n1, bc_num n2, bc_num *result, int scale_min)

{

  bc_num sum = NULL;

  int cmp_res;

  int res_scale;

  if (n1->n_sign == n2->n_sign)

    {

      sum = _bc_do_add (n1, n2, scale_min);

      sum->n_sign = n1->n_sign;

    }

  else

    {

      /* subtraction must be done. */

      cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);  /* Compare magnitudes. */

      switch (cmp_res)

	{

	case -1:

	  /* n1 is less than n2, subtract n1 from n2. */

	  sum = _bc_do_sub (n2, n1, scale_min);

	  sum->n_sign = n2->n_sign;

	  break;

	case  0:

	  /* They are equal! return zero with the correct scale! */

	  res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));

	  sum = bc_new_num (1, res_scale);

	  memset (sum->n_value, 0, res_scale+1);

	  break;

	case  1:

	  /* n2 is less than n1, subtract n2 from n1. */

	  sum = _bc_do_sub (n1, n2, scale_min);

	  sum->n_sign = n1->n_sign;

	}

    }

  /* Clean up and return. */

  bc_free_num (result);

  *result = sum;

}

/* Recursive vs non-recursive multiply crossover ranges. */

#if defined(MULDIGITS)

#include "muldigits.h"

#else

#define MUL_BASE_DIGITS 80

#endif

int mul_base_digits = MUL_BASE_DIGITS;

#define MUL_SMALL_DIGITS mul_base_digits/4

/* Multiply utility routines */

static bc_num

new_sub_num (int length, int scale, char *value)

{

  bc_num temp;

  if (_bc_Free_list != NULL) {

    temp = _bc_Free_list;

    _bc_Free_list = temp->n_next;

  } else {

    temp = (bc_num) malloc (sizeof(bc_struct));

    if (temp == NULL) bc_out_of_memory ();

  }

  temp->n_sign = PLUS;

  temp->n_len = length;

  temp->n_scale = scale;

  temp->n_refs = 1;

  temp->n_ptr = NULL;

  temp->n_value = value;

  return temp;

}

static void

_bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod)

{

  char *n1ptr, *n2ptr, *pvptr;

  char *n1end, *n2end;		/* To the end of n1 and n2. */

  int indx, sum, prodlen;

  prodlen = n1len+n2len+1;

  *prod = bc_new_num (prodlen, 0);

  n1end = (char *) (n1->n_value + n1len - 1);

  n2end = (char *) (n2->n_value + n2len - 1);

  pvptr = (char *) ((*prod)->n_value + prodlen - 1);

  sum = 0;

  /* Here is the loop... */

  for (indx = 0; indx < prodlen-1; indx++)

    {

      n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));

      n2ptr = (char *) (n2end - MIN(indx, n2len-1));

      while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))

	sum += *n1ptr-- * *n2ptr++;

      *pvptr-- = sum % BASE;

      sum = sum / BASE;

    }

  *pvptr = sum;

}

/* A special adder/subtractor for the recursive divide and conquer

   multiply algorithm.  Note: if sub is called, accum must

   be larger that what is being subtracted.  Also, accum and val

   must have n_scale = 0.  (e.g. they must look like integers. *) */

static void

_bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)

{

  signed char *accp, *valp;

  int  count, carry;

  count = val->n_len;

  if (val->n_value[0] == 0)

    count--;

  assert (accum->n_len+accum->n_scale >= shift+count);

  /* Set up pointers and others */

  accp = (signed char *)(accum->n_value +

			 accum->n_len + accum->n_scale - shift - 1);

  valp = (signed char *)(val->n_value + val->n_len - 1);

  carry = 0;

  if (sub) {

    /* Subtraction, carry is really borrow. */

    while (count--) {

      *accp -= *valp-- + carry;

      if (*accp < 0) {

	carry = 1;

        *accp-- += BASE;

      } else {

	carry = 0;

	accp--;

      }

    }

    while (carry) {

      *accp -= carry;

      if (*accp < 0)

	*accp-- += BASE;

      else

	carry = 0;

    }

  } else {

    /* Addition */

    while (count--) {

      *accp += *valp-- + carry;

      if (*accp > (BASE-1)) {

	carry = 1;

        *accp-- -= BASE;

      } else {

	carry = 0;

	accp--;

      }

    }

    while (carry) {

      *accp += carry;

      if (*accp > (BASE-1))

	*accp-- -= BASE;

      else

	carry = 0;

    }

  }

}

/* Recursive divide and conquer multiply algorithm.  

   Based on 

   Let u = u0 + u1*(b^n)

   Let v = v0 + v1*(b^n)

   Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0

   B is the base of storage, number of digits in u1,u0 close to equal.

*/

static void

_bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod)

{ 

  bc_num u0, u1, v0, v1;

  bc_num m1, m2, m3, d1, d2;

  int n, prodlen, m1zero;

  int d1len, d2len;

  /* Base case? */

  if ((ulen+vlen) < mul_base_digits

      || ulen < MUL_SMALL_DIGITS

      || vlen < MUL_SMALL_DIGITS ) {

    _bc_simp_mul (u, ulen, v, vlen, prod);

    return;

  }

  /* Calculate n -- the u and v split point in digits. */

  n = (MAX(ulen, vlen)+1) / 2;

  /* Split u and v. */

  if (ulen < n) {

    u1 = bc_copy_num (_zero_);

    u0 = new_sub_num (ulen,0, u->n_value);

  } else {

    u1 = new_sub_num (ulen-n, 0, u->n_value);

    u0 = new_sub_num (n, 0, u->n_value+ulen-n);

  }

  if (vlen < n) {

    v1 = bc_copy_num (_zero_);

    v0 = new_sub_num (vlen,0, v->n_value);

  } else {

    v1 = new_sub_num (vlen-n, 0, v->n_value);

    v0 = new_sub_num (n, 0, v->n_value+vlen-n);

    }

  _bc_rm_leading_zeros (u1);

  _bc_rm_leading_zeros (u0);

  _bc_rm_leading_zeros (v1);

  _bc_rm_leading_zeros (v0);

  m1zero = bc_is_zero(u1) || bc_is_zero(v1);

  /* Calculate sub results ... */

  bc_init_num(&d1;;

  bc_init_num(&d2;;

  bc_sub (u1, u0, &d1, 0);

  d1len = d1->n_len;

  bc_sub (v0, v1, &d2, 0);

  d2len = d2->n_len;

  /* Do recursive multiplies and shifted adds. */

  if (m1zero)

    m1 = bc_copy_num (_zero_);

  else

    _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1;;

  if (bc_is_zero(d1) || bc_is_zero(d2))

    m2 = bc_copy_num (_zero_);

  else

    _bc_rec_mul (d1, d1len, d2, d2len, &m2;;

  if (bc_is_zero(u0) || bc_is_zero(v0))

    m3 = bc_copy_num (_zero_);

  else

    _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3;;

  /* Initialize product */

  prodlen = ulen+vlen+1;

  *prod = bc_new_num(prodlen, 0);

  if (!m1zero) {

    _bc_shift_addsub (*prod, m1, 2*n, 0);

    _bc_shift_addsub (*prod, m1, n, 0);

  }

  _bc_shift_addsub (*prod, m3, n, 0);

  _bc_shift_addsub (*prod, m3, 0, 0);

  _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);

  /* Now clean up! */

  bc_free_num (&u1;;

  bc_free_num (&u0;;

  bc_free_num (&v1;;

  bc_free_num (&m1;;

  bc_free_num (&v0;;

  bc_free_num (&m2;;

  bc_free_num (&m3;;

  bc_free_num (&d1;;

  bc_free_num (&d2;;

}

/* The multiply routine.  N2 times N1 is put int PROD with the scale of

   the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).

   */

void

bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale)

{

  bc_num pval; 

  int len1, len2;

  int full_scale, prod_scale;

  /* Initialize things. */

  len1 = n1->n_len + n1->n_scale;

  len2 = n2->n_len + n2->n_scale;

  full_scale = n1->n_scale + n2->n_scale;

  prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));

  /* Do the multiply */

  _bc_rec_mul (n1, len1, n2, len2, &pval);

  /* Assign to prod and clean up the number. */

  pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );

  pval->n_value = pval->n_ptr;

  pval->n_len = len2 + len1 + 1 - full_scale;

  pval->n_scale = prod_scale;

  _bc_rm_leading_zeros (pval);

  if (bc_is_zero (pval))

    pval->n_sign = PLUS;

  bc_free_num (prod);

  *prod = pval;

}

/* Some utility routines for the divide:  First a one digit multiply.

   NUM (with SIZE digits) is multiplied by DIGIT and the result is

   placed into RESULT.  It is written so that NUM and RESULT can be

   the same pointers.  */

static void

_one_mult (unsigned char *num, int size, int digit, unsigned char *result)

{

  int carry, value;

  unsigned char *nptr, *rptr;

  if (digit == 0)

    memset (result, 0, size);

  else

    {

      if (digit == 1)

	memcpy (result, num, size);

      else

	{