Logo Search packages:      
Sourcecode: mathomatic version File versions

gcd.c

/*
 * Algebraic manipulator numeric GCD routines.
 *
 * Copyright (c) 1996 George Gesslein II.
 */

#include "am.h"
#include "externs.h"

/*
 * Return the Greatest Common Divisor of doubles "d1" and "d2".
 * Will work with non-integers, but there may be some floating point error.
 * Always works correctly with integers.
 * Returns 0 if it failed.
 */
double
gcd(d1, d2)
double      d1, d2;
{
      int   count;
      double      larger;
      double      divisor;
      double      r1;
      double      d;

      d1 = fabs(d1);
      d2 = fabs(d2);
      if (d1 > d2) {
            larger = d1;
            divisor = d2;
      } else {
            larger = d2;
            divisor = d1;
      }
      d = larger * epsilon;
      if ((2 * d) >= divisor)
            return 0.0;
      for (count = 1; count < 50; count++) {
            r1 = fmod(larger, divisor);
            if (r1 <= d)
                  return divisor;
            if ((divisor - r1) <= d)
                  return divisor;
            larger = divisor;
            divisor = r1;
      }
      return 0.0;
}

/*
 * Round a double to the nearest integer.
 */
double
my_round(d1)
double      d1;
{
      if (d1 >= 0.0) {
            modf(d1 + 0.5, &d1);
      } else {
            modf(d1 - 0.5, &d1);
      }
      return d1;
}

/*
 * Convert the passed double "d" to a fully reduced fraction.
 * "d1p" points to the numerator and "d2p" points to the denominator.
 *
 * Return true with integer numerator and denominator
 * if conversion was successful.
 * True return indicates "d" is rational.
 */
int
f_to_fraction(d, d1p, d2p)
double      d;
double      *d1p;
double      *d2p;
{
      double      divisor;
      double      d1, d2;
      double      k3, k4;

      *d1p = d;
      *d2p = 1.0;
      if (fmod(d, 1.0) == 0.0)
            return true;
      if ((divisor = gcd(1.0, d)) > epsilon) {
            d1 = d / divisor;
            d2 = 1.0 / divisor;
            d1 = my_round(d1);
            d2 = my_round(d2);
            if (fabs(d1) >= 1.0e10)
                  return false;
            if (d2 >= 1.0e10)
                  return false;
            if (d2 < 1.5)
                  return false;
            divisor = gcd(d1, d2);
            if (fmod(divisor, 1.0) != 0.0) {
                  printf(_("Error in gcd() function!\n"));
                  return false;
            }
            if (divisor > 1.0) {
                  d1 = d1 / divisor;
                  d2 = d2 / divisor;
            }
            k3 = (d1 / d2);
            k4 = d;
#if   EPSILON
            if (fabs(k3 - k4) > (small_epsilon * fabs(k3))) {
                  return false;
            }
#else
            if (k3 != k4)
                  return false;
#endif
            *d1p = d1;
            *d2p = d2;
            return true;
      }
      return false;
}

/*
 * Convert non-integer constants in an expression
 * to fractions when appropriate.
 */
make_fractions(equation, np)
token_type  *equation;
int         *np;
{
      int   i, j, k;
      int   level;
      double      d1, d2;
      int   inc_level;

      for (i = 0; i < *np; i += 2) {
            if (equation[i].kind == CONSTANT) {
                  level = equation[i].level;
                  if (i > 0 && equation[i-1].level == level
                      && equation[i-1].token.operatr == DIVIDE)
                        continue;
                  if (!f_to_fraction((double) equation[i].token.constant, &d1, &d2)
                      || d2 == 1.0)
                        continue;
                  if ((*np + 2) > n_tokens) {
                        error_huge();
                  }
                  inc_level = (*np > 1);
                  if ((i + 1) < *np && equation[i+1].level == level) {
                        switch (equation[i+1].token.operatr) {
                        case TIMES:
                              if (d1 == 1.0) {
                                    blt(&equation[i], &equation[i+2], (*np - (i + 2)) * sizeof(*equation));
                                    *np -= 2;
                              } else {
                                    equation[i].token.constant = d1;
                              }
                              for (j = i + 1; j < *np && equation[j].level >= level; j += 2) {
                                    if (equation[j].level == level && equation[j].token.operatr == DIVIDE) {
                                          break;
                                    }
                              }
                              blt(&equation[j+2], &equation[j], (*np - j) * sizeof(*equation));
                              *np += 2;
                              equation[j].level = level;
                              equation[j].kind = OPERATOR;
                              equation[j].token.operatr = DIVIDE;
                              j++;
                              equation[j].level = level;
                              equation[j].kind = CONSTANT;
                              equation[j].token.constant = d2;
                              if (d1 == 1.0) {
                                    i -= 2;
                              }
                              continue;
                        case DIVIDE:
                              inc_level = false;
                        default:
                              break;
                        }
                  }
                  j = i;
                  blt(&equation[i+3], &equation[i+1], (*np - (i + 1)) * sizeof(*equation));
                  *np += 2;
                  equation[j].token.constant = d1;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = DIVIDE;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = d2;
                  if (inc_level) {
                        for (k = i; k <= j; k++)
                              equation[k].level++;
                  }
            }
      }
}

/*
 * Perfect output of GCD and LCM of doubles d1 and d2.
 *
 * Return true if the GCD was found and displayed.
 */
int
ngcd(d1, d2)
double      d1, d2;
{
      double      d3, d4;
      double      d5, d6;

      if (fabs(d1) >= MAX_K_INTEGER || fabs(d2) >= MAX_K_INTEGER) {
            return false;
      }
      d3 = gcd(d1, d2);
      if (d3 == 0.0) {
            printf(_("No GCD found.\n"));
            return false;
      }
      if (fabs(d1 / d3) < (1.0 - epsilon)
          || fabs(d2 / d3) < (1.0 - epsilon)) {
            printf(_("The computed GCD is too large!\n"));
            return false;
      }
      modf(fabs(d1 / d3) + 0.5, &d4);
      d4 = fabs(d1 / d4);
      if (fabs(modf(fabs(d1 / d4), &d5) - 0.5) < (0.5 - epsilon)
          || fabs(modf(fabs(d2 / d4), &d6) - 0.5) < (0.5 - epsilon)) {
            printf(_("The computed GCD is too inaccurate!\n"));
            return false;
      }
      modf(fabs(d1 / d4) + 0.5, &d5);
      modf(fabs(d2 / d4) + 0.5, &d6);
      if (gcd(d5, d6) != 1.0) {
            printf(_("The computed GCD is too small!\n"));
            return false;
      }
      printf(_("Greatest Common Divisor (GCD) = %.14g\n"), d4);
      printf(_("Least Common Multiple (LCM) = %.14g\n"), fabs((d1 * d2) / d4));
      return true;
}

Generated by  Doxygen 1.6.0   Back to index