mpf_sub.c

   1 /* LibTomFloat, multiple-precision floating-point library
   2  *
   3  * LibTomFloat is a library that provides multiple-precision
   4  * floating-point artihmetic as well as trigonometric functionality.
   5  *
   6  * This library requires the public domain LibTomMath to be installed.
   7  *
   8  * This library is free for all purposes without any express
   9  * gurantee it works
  10  *
  11  * Tom St Denis, tomstdenis@iahu.ca, http://float.libtomcrypt.org
  12  */
  13 #include <tomfloat.h>
  14
  15 int  mpf_sub(mp_float *a, mp_float *b, mp_float *c)
  16 {
  17    int err;
  18    mp_float tmp;
  19    long diff;
  20
  21    if (mpf_iszero(a)) {
  22       diff = c->radix;
  23       if ((err = mpf_neg(b, c)) != MP_OKAY) {
  24          return err;
  25       }
  26       return mpf_normalize_to(c, diff);
  27    } else if (mpf_iszero(b)) {
  28       diff = c->radix;
  29       if ((err = mpf_copy(a, c)) != MP_OKAY) {
  30          return err;
  31       }
  32       return mpf_normalize_to(c, diff);
  33    }
  34
  35    if (a->exp < b->exp) {
  36       /* tmp == a normalize to b's exp */
  37       if ((err = mpf_init_copy(a, &tmp)) != MP_OKAY) {
  38          return err;
  39       }
  40
  41       /* now make tmp.exp == b.exp by dividing tmp by 2^(b.exp - tmp.exp) */
  42       diff = b->exp - tmp.exp;
  43       tmp.exp = b->exp;
  44       if ((err = mp_div_2d(&(tmp.mantissa), diff, (&tmp.mantissa), NULL)) != MP_OKAY)  { goto __TMP; }
  45       if ((err = mp_sub(&(tmp.mantissa), &(b->mantissa), &(c->mantissa))) != MP_OKAY)  { goto __TMP; }
  46       c->exp = b->exp;
  47    } else {
  48       /* tmp == b normalize to a's radix */
  49       if ((err = mpf_init_copy(b, &tmp)) != MP_OKAY) {
  50          return err;
  51       }
  52
  53       /* now make tmp.exp == a.exp by dividing tmp by 2^(a.exp - tmp.exp) */
  54       diff = a->exp - tmp.exp;
  55       tmp.exp = a->exp;
  56       if ((err = mp_div_2d(&(tmp.mantissa), diff, (&tmp.mantissa), NULL)) != MP_OKAY)   { goto __TMP; }
  57                 if ((err = mp_sub(&(a->mantissa), &(tmp.mantissa), &(c->mantissa))) != MP_OKAY)  { goto __TMP; }
  58       c->exp = a->exp;
  59    }
  60
  61    err    = mpf_normalize(c);
  62
  63 __TMP:  mpf_clear(&tmp);
  64    return err;
  65 }