## File: gcd.c

package info (click to toggle)
mathomatic 14.0.6-2
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239` ``````/* * General floating point GCD routine and associated code for Mathomatic. * * Copyright (C) 1987-2008 George Gesslein II. */ #include "includes.h" /* * Return the Greatest Common Divisor (GCD) of doubles "d1" and "d2", * by using the Euclidean GCD algorithm. * * The GCD is defined as the largest positive number that evenly divides both "d1" and "d2". * Will usually work with non-integers, but there may be some floating point error. * Always works perfectly with integers up to MAX_K_INTEGER. * * Return 0 on failure or if either parameter is 0, otherwise return the positive GCD. */ double gcd(d1, d2) double d1, d2; { int count; double larger; double divisor; double r1; double d; if (!isfinite(d1) || !isfinite(d2)) { return 0.0; /* operands must be finite */ } d1 = fabs(d1); d2 = fabs(d2); if (d1 > d2) { larger = d1; divisor = d2; } else { larger = d2; divisor = d1; } if (divisor <= 0.0 || larger >= MAX_K_INTEGER) { return 0.0; /* out of range */ } d = larger * epsilon; if (d >= divisor) { return 0.0; /* result would be too inaccurate */ } for (count = 1; count < 50; count++) { r1 = fmod(larger, divisor); if (r1 <= d || (divisor - r1) <= d) { if (r1 != 0.0 && divisor <= (100.0 * d)) return 0.0; return divisor; } larger = divisor; divisor = r1; } return 0.0; } /* * Return the verified exact Greatest Common Divisor (GCD) of doubles "d1" and "d2". * * Return 0 on failure or if either parameter is 0, otherwise return the positive GCD. */ double gcd_verified(d1, d2) double d1, d2; { double divisor; divisor = gcd(d1, d2); if (divisor != 0.0) { if (gcd(d1 / divisor, d2 / divisor) != 1.0) return 0.0; } return divisor; } /* * Return a floating point double rounded to the nearest integer. */ double my_round(d1) double d1; /* value to round */ { 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. * This done by the following simple algorithm: * * divisor = gcd(d, 1.0) * numerator = d / divisor * denominator = 1.0 / divisor * * Return true with integers in numerator and denominator if conversion was successful. * Otherwise return false with numerator = "d" and denominator = "1.0". * * True return indicates "d" is rational and finite, otherwise "d" is probably irrational. */ int f_to_fraction(d, numeratorp, denominatorp) double d; /* floating point number to convert */ double *numeratorp; /* returned numerator */ double *denominatorp; /* returned denominator */ { double divisor; double numerator, denominator; double k3, k4; *numeratorp = d; *denominatorp = 1.0; if (!isfinite(d)) { return false; } /* see if "d" is an integer, or very close to an integer: */ if (fmod(d, 1.0) == 0.0) { return true; } k3 = fabs(d) * small_epsilon; k4 = my_round(d); if (k4 != 0.0 && fabs(k4 - d) <= k3) { *numeratorp = k4; return true; } /* try to convert non-integer floating point value in "d" to a fraction: */ if ((divisor = gcd(1.0, d)) > epsilon) { numerator = my_round(d / divisor); denominator = my_round(1.0 / divisor); /* don't allow more than 11 digits in the numerator or denominator: */ if (fabs(numerator) >= 1.0e12) return false; if (denominator >= 1.0e12 || denominator < 2.0) return false; /* make sure the result is a fully reduced fraction: */ divisor = gcd(numerator, denominator); if (divisor > 1.0) { /* shouldn't happen, but just in case */ numerator /= divisor; denominator /= divisor; } k3 = (numerator / denominator); k4 = d; if (fabs(k3 - k4) > (small_epsilon * fabs(k3))) { return false; /* result is too inaccurate */ } /* numerator and denominator are integral */ *numeratorp = numerator; *denominatorp = denominator; return true; } return false; } /* * Convert non-integer constants in an equation side to fractions, when appropriate. */ void make_fractions(equation, np) token_type *equation; /* equation side pointer */ int *np; /* pointer to length of equation side */ { int i, j, k; int level; double numerator, denominator; 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(equation[i].token.constant, &numerator, &denominator)) continue; if (denominator == 1.0) { equation[i].token.constant = numerator; 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: for (j = i + 3; j < *np && equation[j].level >= level; j += 2) { if (equation[j].level == level && equation[j].token.operatr == DIVIDE) { break; } } if (numerator == 1.0) { blt(&equation[i], &equation[i+2], (j - (i + 2)) * sizeof(token_type)); j -= 2; } else { equation[i].token.constant = numerator; blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type)); *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 = denominator; if (numerator == 1.0) { i -= 2; } continue; case DIVIDE: inc_level = false; break; } } j = i; blt(&equation[i+3], &equation[i+1], (*np - (i + 1)) * sizeof(token_type)); *np += 2; equation[j].token.constant = numerator; 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 = denominator; if (inc_level) { for (k = i; k <= j; k++) equation[k].level++; } } } } ``````