/* LibTomFloat, multiple-precision floating-point library
 *
 * LibTomFloat is a library that provides multiple-precision
 * floating-point artihmetic as well as trigonometric functionality.
 *
 * This library requires the public domain LibTomMath to be installed.
 * 
 * This library is free for all purposes without any express
 * gurantee it works
 *
 * Tom St Denis, tomstdenis@iahu.ca, http://float.libtomcrypt.org
 */
#include <tomfloat.h>

int  mpf_sub(mp_float *a, mp_float *b, mp_float *c)
{
   int err;
   mp_float tmp;
   long diff;

   if (mpf_iszero(a)) {
      diff = c->radix;
      if ((err = mpf_neg(b, c)) != MP_OKAY) {
         return err;
      }
      return mpf_normalize_to(c, diff);
   } else if (mpf_iszero(b)) {
      diff = c->radix;
      if ((err = mpf_copy(a, c)) != MP_OKAY) {
         return err;
      }
      return mpf_normalize_to(c, diff);
   }

   if (a->exp < b->exp) {
      /* tmp == a normalize to b's exp */
      if ((err = mpf_init_copy(a, &tmp)) != MP_OKAY) {
         return err;
      }

      /* now make tmp.exp == b.exp by dividing tmp by 2^(b.exp - tmp.exp) */
      diff = b->exp - tmp.exp;
      tmp.exp = b->exp;
      if ((err = mp_div_2d(&(tmp.mantissa), diff, (&tmp.mantissa), NULL)) != MP_OKAY)  { goto __TMP; }
      if ((err = mp_sub(&(tmp.mantissa), &(b->mantissa), &(c->mantissa))) != MP_OKAY)  { goto __TMP; }
      c->exp = b->exp;
   } else {
      /* tmp == b normalize to a's radix */
      if ((err = mpf_init_copy(b, &tmp)) != MP_OKAY) {
         return err;
      }

      /* now make tmp.exp == a.exp by dividing tmp by 2^(a.exp - tmp.exp) */
      diff = a->exp - tmp.exp;
      tmp.exp = a->exp;
      if ((err = mp_div_2d(&(tmp.mantissa), diff, (&tmp.mantissa), NULL)) != MP_OKAY)   { goto __TMP; }
		if ((err = mp_sub(&(a->mantissa), &(tmp.mantissa), &(c->mantissa))) != MP_OKAY)  { goto __TMP; }
      c->exp = a->exp;
   }

   err    = mpf_normalize(c);

__TMP:  mpf_clear(&tmp);
   return err;
}