/*
 *  M_APM  -  mapm_div.c
 *
 *  Copyright (C) 1999 - 2007   Michael C. Ring
 *
 *  Permission to use, copy, and distribute this software and its
 *  documentation for any purpose with or without fee is hereby granted,
 *  provided that the above copyright notice appear in all copies and
 *  that both that copyright notice and this permission notice appear
 *  in supporting documentation.
 *
 *  Permission to modify the software is granted. Permission to distribute
 *  the modified code is granted. Modifications are to be distributed by
 *  using the file 'license.txt' as a template to modify the file header.
 *  'license.txt' is available in the official MAPM distribution.
 *
 *  This software is provided "as is" without express or implied warranty.
 */

/*
 *
 *      This file contains the basic division functions
 *
 */

#include "pgAdmin3.h"
#include "pgscript/utilities/mapm-lib/m_apm_lc.h"

static	M_APM	M_div_worka;
static	M_APM	M_div_workb;
static	M_APM	M_div_tmp7;
static	M_APM	M_div_tmp8;
static	M_APM	M_div_tmp9;

static	int	M_div_firsttime = TRUE;

/****************************************************************************/
void	M_free_all_div()
{
	if (M_div_firsttime == FALSE)
	{
		m_apm_free(M_div_worka);
		m_apm_free(M_div_workb);
		m_apm_free(M_div_tmp7);
		m_apm_free(M_div_tmp8);
		m_apm_free(M_div_tmp9);

		M_div_firsttime = TRUE;
	}
}
/****************************************************************************/
void	m_apm_integer_div_rem(M_APM qq, M_APM rr, M_APM aa, M_APM bb)
{
	m_apm_integer_divide(qq, aa, bb);
	m_apm_multiply(M_div_tmp7, qq, bb);
	m_apm_subtract(rr, aa, M_div_tmp7);
}
/****************************************************************************/
void	m_apm_integer_divide(M_APM rr, M_APM aa, M_APM bb)
{
	/*
	 *    we must use this divide function since the
	 *    faster divide function using the reciprocal
	 *    will round the result (possibly changing
	 *    nnm.999999...  -->  nn(m+1).0000 which would
	 *    invalidate the 'integer_divide' goal).
	 */

	M_apm_sdivide(rr, 4, aa, bb);

	if (rr->m_apm_exponent <= 0)        /* result is 0 */
	{
		M_set_to_zero(rr);
	}
	else
	{
		if (rr->m_apm_datalength > rr->m_apm_exponent)
		{
			rr->m_apm_datalength = rr->m_apm_exponent;
			M_apm_normalize(rr);
		}
	}
}
/****************************************************************************/
void	M_apm_sdivide(M_APM r, int places, M_APM a, M_APM b)
{
	int	j, k, m, b0, sign, nexp, indexr, icompare, iterations;
	long    trial_numer;
	void	*vp;

	if (M_div_firsttime)
	{
		M_div_firsttime = FALSE;

		M_div_worka = m_apm_init();
		M_div_workb = m_apm_init();
		M_div_tmp7  = m_apm_init();
		M_div_tmp8  = m_apm_init();
		M_div_tmp9  = m_apm_init();
	}

	sign = a->m_apm_sign * b->m_apm_sign;

	if (sign == 0)      /* one number is zero, result is zero */
	{
		if (b->m_apm_sign == 0)
		{
			M_apm_log_error_msg(M_APM_RETURN, "\'M_apm_sdivide\', Divide by 0");
		}

		M_set_to_zero(r);
		return;
	}

	/*
	 *  Knuth step D1. Since base = 100, base / 2 = 50.
	 *  (also make the working copies positive)
	 */

	if (b->m_apm_data[0] >= 50)
	{
		m_apm_absolute_value(M_div_worka, a);
		m_apm_absolute_value(M_div_workb, b);
	}
	else       /* 'normal' step D1 */
	{
		k = 100 / (b->m_apm_data[0] + 1);
		m_apm_set_long(M_div_tmp9, (long)k);

		m_apm_multiply(M_div_worka, M_div_tmp9, a);
		m_apm_multiply(M_div_workb, M_div_tmp9, b);

		M_div_worka->m_apm_sign = 1;
		M_div_workb->m_apm_sign = 1;
	}

	/* setup trial denominator for step D3 */

	b0 = 100 * (int)M_div_workb->m_apm_data[0];

	if (M_div_workb->m_apm_datalength >= 3)
		b0 += M_div_workb->m_apm_data[1];

	nexp = M_div_worka->m_apm_exponent - M_div_workb->m_apm_exponent;

	if (nexp > 0)
		iterations = nexp + places + 1;
	else
		iterations = places + 1;

	k = (iterations + 1) >> 1;     /* required size of result, in bytes */

	if (k > r->m_apm_malloclength)
	{
		if ((vp = MAPM_REALLOC(r->m_apm_data, (k + 32))) == NULL)
		{
			/* fatal, this does not return */

			M_apm_log_error_msg(M_APM_FATAL, "\'M_apm_sdivide\', Out of memory");
		}

		r->m_apm_malloclength = k + 28;
		r->m_apm_data = (UCHAR *)vp;
	}

	/* clear the exponent in the working copies */

	M_div_worka->m_apm_exponent = 0;
	M_div_workb->m_apm_exponent = 0;

	/* if numbers are equal, ratio == 1.00000... */

	if ((icompare = m_apm_compare(M_div_worka, M_div_workb)) == 0)
	{
		iterations = 1;
		r->m_apm_data[0] = 10;
		nexp++;
	}
	else			           /* ratio not 1, do the real division */
	{
		if (icompare == 1)                        /* numerator > denominator */
		{
			nexp++;                           /* to adjust the final exponent */
			M_div_worka->m_apm_exponent += 1;     /* multiply numerator by 10 */
		}
		else                                      /* numerator < denominator */
		{
			M_div_worka->m_apm_exponent += 2;    /* multiply numerator by 100 */
		}

		indexr = 0;
		m      = 0;

		while (TRUE)
		{
			/*
			 *  Knuth step D3. Only use the 3rd -> 6th digits if the number
			 *  actually has that many digits.
			 */

			trial_numer = 10000L * (long)M_div_worka->m_apm_data[0];

			if (M_div_worka->m_apm_datalength >= 5)
			{
				trial_numer += 100 * M_div_worka->m_apm_data[1]
				               + M_div_worka->m_apm_data[2];
			}
			else
			{
				if (M_div_worka->m_apm_datalength >= 3)
					trial_numer += 100 * M_div_worka->m_apm_data[1];
			}

			j = (int)(trial_numer / b0);

			/*
			 *    Since the library 'normalizes' all the results, we need
			 *    to look at the exponent of the number to decide if we
			 *    have a lead in 0n or 00.
			 */

			if ((k = 2 - M_div_worka->m_apm_exponent) > 0)
			{
				while (TRUE)
				{
					j /= 10;
					if (--k == 0)
						break;
				}
			}

			if (j == 100)     /* qhat == base ??      */
				j = 99;         /* if so, decrease by 1 */

			m_apm_set_long(M_div_tmp8, (long)j);
			m_apm_multiply(M_div_tmp7, M_div_tmp8, M_div_workb);

			/*
			 *    Compare our q-hat (j) against the desired number.
			 *    j is either correct, 1 too large, or 2 too large
			 *    per Theorem B on pg 272 of Art of Compter Programming,
			 *    Volume 2, 3rd Edition.
			 *
			 *    The above statement is only true if using the 2 leading
			 *    digits of the numerator and the leading digit of the
			 *    denominator. Since we are using the (3) leading digits
			 *    of the numerator and the (2) leading digits of the
			 *    denominator, we eliminate the case where our q-hat is
			 *    2 too large, (and q-hat being 1 too large is quite remote).
			 */

			if (m_apm_compare(M_div_tmp7, M_div_worka) == 1)
			{
				j--;
				m_apm_subtract(M_div_tmp8, M_div_tmp7, M_div_workb);
				m_apm_copy(M_div_tmp7, M_div_tmp8);
			}

			/*
			 *  Since we know q-hat is correct, step D6 is unnecessary.
			 *
			 *  Store q-hat, step D5. Since D6 is unnecessary, we can
			 *  do D5 before D4 and decide if we are done.
			 */

			r->m_apm_data[indexr++] = (UCHAR)j;    /* j == 'qhat' */
			m += 2;

			if (m >= iterations)
				break;

			/* step D4 */

			m_apm_subtract(M_div_tmp9, M_div_worka, M_div_tmp7);

			/*
			 *  if the subtraction yields zero, the division is exact
			 *  and we are done early.
			 */

			if (M_div_tmp9->m_apm_sign == 0)
			{
				iterations = m;
				break;
			}

			/* multiply by 100 and re-save */
			M_div_tmp9->m_apm_exponent += 2;
			m_apm_copy(M_div_worka, M_div_tmp9);
		}
	}

	r->m_apm_sign       = sign;
	r->m_apm_exponent   = nexp;
	r->m_apm_datalength = iterations;

	M_apm_normalize(r);
}
/****************************************************************************/