/[dtapublic]/projs/trunk/shared_source/c_datd/gmp_rats.c

Diff of /projs/trunk/shared_source/c_datd/gmp_rats.c

Parent Directory | Revision Log | View Patch Patch

-revision 56 by dashley,
Sat Oct 29 01:53:01 2016 UTC
+revision 71 by dashley,
Sat Nov  5 11:07:06 2016 UTC
 Line 1
  //$Header$
  //-------------------------------------------------------------------------------------------------
  //This file is part of "David T. Ashley's Shared Source Code", a set of shared components
  //integrated into many of David T. Ashley's projects.
  //-------------------------------------------------------------------------------------------------
  //This source code and any program in which it is compiled/used is provided under the MIT License,
  //reproduced below.
  //-------------------------------------------------------------------------------------------------
  //Permission is hereby granted, free of charge, to any person obtaining a copy of
  //this software and associated documentation files(the "Software"), to deal in the
  //Software without restriction, including without limitation the rights to use,
  //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
  //Software, and to permit persons to whom the Software is furnished to do so,
  //subject to the following conditions :
  //
  //The above copyright notice and this permission notice shall be included in all
  //copies or substantial portions of the Software.
  //
  //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
  //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  //SOFTWARE.
  //-------------------------------------------------------------------------------------------------
  #define MODULE_GMP_RATS
  #include <assert.h>
  #include <stdio.h>
  #include <string.h>
  #include "bstrfunc.h"
  #include "charfunc.h"
  #include "gmp_ints.h"
  #include "gmp_rats.h"
  #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
     #include "ccmalloc.h"
  #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
     #include "tclalloc.h"
  #else
     #include <malloc.h>
  #endif
  /******************************************************************/
  /***  STATUS FUNCTIONS  *******************************************/
  /******************************************************************/
  //Functions in this category provide information about rational
  //numbers.
  //08/08/01:  Visual inspection OK.
  int GMP_RATS_mpq_is_nan(const GMP_RATS_mpq_struct *rn)
     {
     assert(rn != NULL);
     //A rational number is NAN in one of two
     //circumstances.  If either of the integer components
     //is NAN, or else if there is a zero denominator.
     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
        {
        return(1);
        }
     if (GMP_INTS_mpz_is_zero(&(rn->den)))
        {
        return(1);
        }
     //We're clean ...
     return(0);
     }
  /******************************************************************/
  /***  INITIALIZATION, CLEARING, AND SETTING FUNCTIONS  ************/
  /******************************************************************/
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_init(GMP_RATS_mpq_struct *arg)
     {
     //Eyeball the input parameter.
     assert(arg != NULL);
     //Initialize the numerator and denominator.
     GMP_INTS_mpz_init(&(arg->num));
     GMP_INTS_mpz_init(&(arg->den));
     //Canonically, we must start off as 0/1--canonical zero.
     GMP_INTS_mpz_set_ui(&(arg->num), 0);
     GMP_INTS_mpz_set_ui(&(arg->den), 1);
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_clear(GMP_RATS_mpq_struct *arg)
     {
     //Eyeball the input parameter.
     assert(arg != NULL);
     //Clear the numerator and denominator.  The called functions
     //will check for NULL pointers and so forth.
     GMP_INTS_mpz_clear(&(arg->num));
     GMP_INTS_mpz_clear(&(arg->den));
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_set_si(GMP_RATS_mpq_struct *arg,
                           int num,
                           int den)
     {
     //Eyeball the input parameters.
     assert(arg != NULL);
     //Set the numerator and denominator.
     GMP_INTS_mpz_set_si(&(arg->num), num);
     GMP_INTS_mpz_set_si(&(arg->den), den);
     }
  //08/08/01:  Visual inspection OK.
  void GMP_RATS_mpq_copy(      GMP_RATS_mpq_struct *dst,
                         const GMP_RATS_mpq_struct *src)
     {
     assert(dst != NULL);
     assert(src != NULL);
     GMP_INTS_mpz_copy(&(dst->num), &(src->num));
     GMP_INTS_mpz_copy(&(dst->den), &(src->den));
     }
  //08/13/01:  Visual inspection OK.
  void GMP_RATS_mpq_swap( GMP_RATS_mpq_struct *a,
                          GMP_RATS_mpq_struct *b)
     {
     assert(a != NULL);
     assert(b != NULL);
     //Handle the swap by swapping integer components.
     GMP_INTS_mpz_swap(&(a->num), &(b->num));
     GMP_INTS_mpz_swap(&(a->den), &(b->den));
     }
  //08/13/01:  Visual inspection OK.
  void GMP_RATS_mpq_swap_components(GMP_RATS_mpq_struct *arg)
     {
     assert(arg != NULL);
     GMP_INTS_mpz_swap(&(arg->num), &(arg->den));
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_set_complex_slash_sepd_rat_num(const char *s,
                                                   int *failure,
                                                   GMP_RATS_mpq_struct *rn)
     {
     char *slash_posn, *numerator, *denominator;
     int s_len, numerator_len, denominator_len;
     int i;
     //Eyeball the input parameters.
     assert(s != NULL);
     assert(failure != NULL);
     assert(rn != NULL);
     //Start off believing there is no failure.
     *failure = 0;
     //Figure out if there is one and only one slash in the
     //string.  If this condition isn't met, we cannot
     //go further.
     slash_posn = strchr(s, '/');
     if (!slash_posn)
        {
        *failure = 1;
        return;
        }
     if (strchr(slash_posn + 1, '/')) //There is a second occurence.
        {
        *failure = 1;
        return;
        }
     //At this point we have one and only one slash.
     //Crack the string in two.  We must do this because the
     //input is a constant string.  We are not allowed to touch it
     //in the logical domain because of the "const" keyword.  We can't
     //do this in the physical domain because the debugger will nail
     //us for it.
     s_len = strlen(s);
     numerator_len = slash_posn - s;
     denominator_len = strlen(slash_posn + 1);
     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
        numerator   = CCMALLOC_malloc(sizeof(char) * (numerator_len   + 1));
        denominator = CCMALLOC_malloc(sizeof(char) * (denominator_len + 1));
     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
        numerator   = TclpAlloc(sizeof(char) * (numerator_len   + 1));
        denominator = TclpAlloc(sizeof(char) * (denominator_len + 1));
     #else
        numerator   = malloc(sizeof(char) * (numerator_len   + 1));
        denominator = malloc(sizeof(char) * (denominator_len + 1));
     #endif
     assert(numerator   != NULL);
     assert(denominator != NULL);
     for (i=0; i<numerator_len; i++)
        {
        numerator[i] = s[i];
        }
     numerator[numerator_len] = 0;
     for (i=0; i<denominator_len; i++)
        {
        denominator[i] = s[slash_posn - s + 1 + i];
        }
     denominator[denominator_len] = 0;
     //Try to parse out the numerator as an arbitrary integer.
     //If this can't be done, it is an immediate failure.
     GMP_INTS_mpz_set_general_int(&(rn->num),
                                  failure,
                                  numerator);
     if (*failure)
        {
        *failure = 1;  //Clamp to 1, don't know what non-zero value
                       //was there.
        goto ret_pt;
        }
     //Try to parse out the denominator.
     GMP_INTS_mpz_set_general_int(&(rn->den),
                                  failure,
                                  denominator);
     if (*failure)
        {
        *failure = 1;  //Clamp to 1, don't know what non-zero value
                       //was there.
        goto ret_pt;
        }
     //At this point, we have both a numerator and denominator.
     //Clean up and return.
     ret_pt:
     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
        CCMALLOC_free(numerator);
        CCMALLOC_free(denominator);
     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
        TclpFree(numerator);
        TclpFree(denominator);
     #else
        free(numerator);
        free(denominator);
     #endif
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_set_sci_not_rat_num(const char *s,
                                        int *failure,
                                        GMP_RATS_mpq_struct *rn)
     {
     int parse_failure;
        //Return code from the floating point parsing
        //function.
     char mant_sign;
        //Sign character, if any, from the mantissa,
        //or N otherwise.
     size_t mant_bdp;
        //The index to the start of the mantissa before
        //the decimal point.
     size_t mant_bdp_len;
        //The length of the mantissa before the decimal
        //point.  Zero means not defined, i.e. that
        //no characters were parsed and interpreted as
        //that part of a floating point number.
     size_t mant_adp;
     size_t mant_adp_len;
        //Similar fields for after the decimal point.
     char exp_sign;
        //Sign of the exponent, if any, or N otherwise.
     size_t exp;
     size_t exp_len;
        //Similar fields as to the mantissa, but for the
        //exponent.
     size_t si;
        //Iteration variable.
     int exponent_val;
        //The value of the exponent.  We can't accept
        //an exponent outside the range of a 24-bit
        //signed integer.  The 24-bit limit is arbitrary.
        //For one thing, it gives room to detect overflow
        //as are adding and multiplying by 10.
     //Eyeball the input parameters.
     assert(s != NULL);
     assert(failure != NULL);
     assert(rn != NULL);
        //Subcomponents of the rational number will be checked as
        //we make integer calls, if we're in debug mode.
     //Start off believing no failure.
     *failure = 0;
     //Set the output to 0/1.  This is the default case for some
     //steps below.
     GMP_RATS_mpq_set_si(rn, 0, 1);
     //Attempt to parse the number as a general number
     //in scientific notation.
     BSTRFUNC_parse_gen_sci_not_num(s,
                                    &parse_failure,
                                    &mant_sign,
                                    &mant_bdp,
                                    &mant_bdp_len,
                                    &mant_adp,
                                    &mant_adp_len,
                                    &exp_sign,
                                    &exp,
                                    &exp_len);
     //If it wouldn't parse as a general number, can't go further.
     if (parse_failure)
        {
        *failure = 1;
        return;
        }
     else
        {
        //The number parsed out.  The general strategy is to form a rational number
        //consisting of the mantissa, with the decimal point shifted fully right, over
        //a denominator of 1.  From there, we process the exponent and combine it with
        //the number of characters after the decimal point to form a virtual exponent.
        //If the exponent is positive, we multiply the numerator by the power of 10.
        //If the exponent is negative, we multiply the denominator by that power of 10.
        //We want to trim the trailing zeros off of the portion of the mantissa after the
        //decimal point.  We only need to back up indices, no need to make copies, etc.
        //Note that it is possible that there are only zeros, in which case we'll end
        //up with a length of zero.
        while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))
           mant_adp_len--;
        //Trim the leading zeros off of the portion of the mantissa before the
        //decimal point.  Note that it is possible that there is only a zero,
        //so we may trim it down to nothing.
        while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))
           {
           mant_bdp++;
           mant_bdp_len--;
           }
        //If we have only zeros in the mantissa, both before the
        //decimal point and after, then we return 0.
        if ((mant_bdp_len + mant_adp_len) == 0)
           {
           *failure = 0;
           return;
           }
        //Convert the numerator to an integer which represents the
        //part before the mantissa and the part after the mantissa
        //concatenated as an integer.  We could call a function to do
        //this, but the function is not really any better in algorithm.
        //We can do it ourselves.
        GMP_INTS_mpz_set_ui(&(rn->num), 0);
        for (si = 0; si < mant_bdp_len; si++)
           {
           int val;
           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);
           val = CHARFUNC_digit_to_val(s[mant_bdp + si]);
           if (val >= 0)
              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);
           }
        for (si = 0; si < mant_adp_len; si++)
           {
           int val;
           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);
           val = CHARFUNC_digit_to_val(s[mant_adp + si]);
           if (val >= 0)
              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);
           }
        //The numerator should now have an integer which is
        //The absolute value of the mantissa.  Process the possible
        //sign.
        if (mant_sign == '-')
           GMP_INTS_mpz_negate(&(rn->num));
        //We now need to form a value from the exponent, if any.
        //First, tackle the exponent.  Process the
        //exponent into a signed integer.  We have to
        //balk at anything outside of 24 bits.  The
        //procedure used automatically handles
        //leading zeros correctly.
        exponent_val = 0;
        for (si=exp; si<(exp+exp_len); si++)
           {
           int val;
           val = CHARFUNC_digit_to_val(s[si]);
           assert(val >= 0 && val <= 9);
           exponent_val *= 10;
           exponent_val += val;
           if (((exp_sign=='-') && (exponent_val>8388608))
              ||
              ((exp_sign != '-') && (exponent_val>8388607)))
              {
              *failure = 1;
              return;
              }
           }
        //If we're here, the exponent has been computed and
        //is within 24 bits.  However, we need to adjust for
        //the sign.
        if (exp_sign == '-')
           exponent_val = -exponent_val;
        //We need to adjust the exponent for the number of digits
        //after the decimal point.
        exponent_val -= mant_adp_len;
        //Again, clip for size.
        if ((exponent_val < -8388608) || (exponent_val > 8388607))
           {
           *failure = 1;
           return;
           }
        //There are two cases to consider.  If the exponent
        //is positive, we need to multiply the numerator
        //by 10 exponentiated to the power of the exponent.
        //If the exponent is negative, we need to do the
        //same thing to the denominator.  If the exponent
        //is negative, we don't need to do anything.
        if (exponent_val > 0)
           {
           GMP_INTS_mpz_struct k10, k10_exponentiated;
           GMP_INTS_mpz_init(&k10);
           GMP_INTS_mpz_init(&k10_exponentiated);
           GMP_INTS_mpz_set_ui(&k10, 10);
           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, exponent_val);
           GMP_INTS_mpz_mul(&(rn->num), &(rn->num), &k10_exponentiated);
           GMP_INTS_mpz_clear(&k10);
           GMP_INTS_mpz_clear(&k10_exponentiated);
           *failure = 0;
           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
              *failure = 1;
           return;
           }
        else if (exponent_val < 0)
           {
           GMP_INTS_mpz_struct k10, k10_exponentiated;
           GMP_INTS_mpz_init(&k10);
           GMP_INTS_mpz_init(&k10_exponentiated);
           GMP_INTS_mpz_set_ui(&k10, 10);
           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, -exponent_val);
           GMP_INTS_mpz_mul(&(rn->den), &(rn->den), &k10_exponentiated);
           GMP_INTS_mpz_clear(&k10);
           GMP_INTS_mpz_clear(&k10_exponentiated);
           *failure = 0;
           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
              *failure = 1;
           return;
           }
        }
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_set_all_format_rat_num(const char *s,
                                           int *failure,
                                           GMP_RATS_mpq_struct *rn)
     {
     //Eyeball the input parameters.
     assert(s != NULL);
     assert(failure != NULL);
     assert(rn != NULL);
     //Assume no failure.
     *failure = 0;
     //Try in order to parse as integers with slash then
     //as number in scientific notation.
     GMP_RATS_mpq_set_complex_slash_sepd_rat_num(s,
                                                 failure,
                                                 rn);
     if (!*failure)
        return;
     GMP_RATS_mpq_set_sci_not_rat_num(s,
                                      failure,
                                      rn);
     if (*failure)
        *failure = 1;  //Clamp output.
     }
  /******************************************************************/
  /***  NORMALIZATION FUNCTIONS  ************************************/
  /******************************************************************/
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_normalize_sign(GMP_RATS_mpq_struct *rn)
     {
     //Eyeball the input.
     assert(rn != NULL);
     if (GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))
        {
        //Both negative, can negate both, this leaves both positive,
        //which is the normalized form for a positive rational
        //number.
        GMP_INTS_mpz_negate(&rn->num);
        GMP_INTS_mpz_negate(&rn->den);
        }
     else if (!GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))
        {
        //Denominator neg, numerator non-neg, can negate both.  This
        //will leave numerator neg, denominator pos, which is
        //normalized form for negative rational number.
        GMP_INTS_mpz_negate(&rn->num);
        GMP_INTS_mpz_negate(&rn->den);
        }
     }
  //08/07/01:  Visual inspection OK.
  void GMP_RATS_mpq_normalize(GMP_RATS_mpq_struct *rn)
     {
     //Eyeball the input.
     assert(rn != NULL);
     //Cover some special cases.  If either component has flags
     //set, don't even touch it.
     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
        {
        return;
        }
     //If the denominator is zero, normalize it to 1/0, the canonical
     //for for an illegal rational number.
     else if (GMP_INTS_mpz_is_zero(&(rn->den)))
        {
        GMP_RATS_mpq_set_si(rn, 1, 0);
        return;
        }
     //If the numerator is zero, convert the number to the canonical
     //form for zero of 0/1.
     else if (GMP_INTS_mpz_is_zero(&(rn->num)))
        {
        GMP_RATS_mpq_set_si(rn, 0, 1);
        return;
        }
     else
        {
        int num_is_neg;
        int den_is_neg;
        GMP_INTS_mpz_struct gcd, quotient, remainder;
        //Allocate space for the integers used.
        GMP_INTS_mpz_init(&gcd);
        GMP_INTS_mpz_init(&quotient);
        GMP_INTS_mpz_init(&remainder);
        //This is the most normal case, where we need to
        //look at reducing the numerator and denominator.
        //One way to do it would be to obtain the g.c.d.
        //and divide this out, and this is the route
        //we'll take.  However, must grab out the sign.
        if (GMP_INTS_mpz_is_neg(&(rn->num)))
           {
           num_is_neg = 1;
           GMP_INTS_mpz_negate(&(rn->num));
           }
        else
           {
           num_is_neg = 0;
           }
        if (GMP_INTS_mpz_is_neg(&(rn->den)))
           {
           den_is_neg = 1;
           GMP_INTS_mpz_negate(&(rn->den));
           }
        else
           {
           den_is_neg = 0;
           }
        //Calculate the GCD.
        GMP_INTS_mpz_gcd(&gcd, &(rn->num), &(rn->den));
        //Divide the numerator by the GCD and store it
        //back.
        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,
                             &(rn->num), &gcd);
        GMP_INTS_mpz_copy(&(rn->num), &quotient);
        //Divide the denominator by the GCD and store it
        //back.
        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,
                             &(rn->den), &gcd);
        GMP_INTS_mpz_copy(&(rn->den), &quotient);
        //We now need to adjust the sign.  Both the
        //numerator and denominator are definitely
        //positive.  Need to make the numerator
        //negative if either but not both of the
        //original signs were negative.
        if ((num_is_neg && !den_is_neg) || (!num_is_neg && den_is_neg))
           {
           GMP_INTS_mpz_negate(&(rn->num));
           }
        //Deallocate space for the integers used.
        GMP_INTS_mpz_clear(&gcd);
        GMP_INTS_mpz_clear(&quotient);
        GMP_INTS_mpz_clear(&remainder);
        return;
        }
     }
  /******************************************************************/
  /***  ARITHMETIC FUNCTIONS  ***************************************/
  /******************************************************************/
  //08/08/01:  Visual inspection OK.
  void GMP_RATS_mpq_add(      GMP_RATS_mpq_struct *result,
                        const GMP_RATS_mpq_struct *arg1,
                        const GMP_RATS_mpq_struct *arg2)
     {
     GMP_RATS_mpq_struct rv;
     GMP_INTS_mpz_struct temp;
     //Eyeball the input parameters.
     assert(result != NULL);
     assert(arg1   != NULL);
     assert(arg2   != NULL);
     //Generally speaking, we do not want to require that
     //the arguments and the result be distinct, as this is
     //too much of a restriction on the caller.  The approach
     //taken, somewhat wasteful, is to allocate a place for
     //the return value.
     //
     //For addition, if we are adding a/b and c/d, the
     //result is necessarily algebraically
     //(ad + cb)/bd.
     //
     //If either rational number in the input is invalid,
     //flag the result as invalid.
     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
        {
        GMP_RATS_mpq_set_si(result, 1, 0);
        }
     else
        {
        //Both rational numbers are OK.  Can simply stage the
        //result by the algebraic identity and then
        //normalize it.  Only need one temporary variable.
        //
        //Initialize the rational number that we will use to
        //hold return value in case it is the same as one
        //or both of the arguments.
        GMP_RATS_mpq_init(&rv);
        //Initialize the temporary integer.
        GMP_INTS_mpz_init(&temp);
        //numerator = a * d
        GMP_INTS_mpz_mul(&(rv.num), &(arg1->num), &(arg2->den));
        //temp = c * b
        GMP_INTS_mpz_mul(&temp, &(arg2->num), &(arg1->den));
        //numerator = a * d + c * b
        GMP_INTS_mpz_add(&(rv.num), &(rv.num), &temp);
        //denominator = b * d
        GMP_INTS_mpz_mul(&(rv.den), &(arg1->den), &(arg2->den));
        //Copy the temporary result to the actual return value.
        //Had to wait until now in case result was the same
        //as either or both args.
        GMP_RATS_mpq_copy(result, &rv);
        //Normalize the result.
        GMP_RATS_mpq_normalize(result);
        //Free dynamic memory.
        GMP_RATS_mpq_clear(&rv);
        GMP_INTS_mpz_clear(&temp);
        }
     }
  //08/08/01:  Visual inspection OK.
  void GMP_RATS_mpq_sub(      GMP_RATS_mpq_struct *result,
                        const GMP_RATS_mpq_struct *arg1,
                        const GMP_RATS_mpq_struct *arg2)
     {
     GMP_RATS_mpq_struct negated_arg_2;
     //Eyeball the input parameters.
     assert(result != NULL);
     assert(arg1   != NULL);
     assert(arg2   != NULL);
     //For the subtract function, we could do it directly,
     //but might as well just define it recursively
     //in terms of add.  We can't modify the inputs,
     //so copy the second off and negate it.  All error
     //flags and so forth will propagate automatically.
     //
     //Allocate space for the negated arg 2.
     GMP_RATS_mpq_init(&negated_arg_2);
     //Copy from the original.
     GMP_RATS_mpq_copy(&negated_arg_2, arg2);
     //Negate the copy.  Negating the numerator will
     //do it.
     GMP_INTS_mpz_negate(&(negated_arg_2.num));
     //Make the add, which now is really a subtract.
     GMP_RATS_mpq_add(result, arg1, &negated_arg_2);
     //Destroy the temporary variable.
     GMP_RATS_mpq_clear(&negated_arg_2);
     }
  //08/16/01: Visual inspection OK.
  void GMP_RATS_mpq_mul(      GMP_RATS_mpq_struct *result,
                        const GMP_RATS_mpq_struct *arg1,
                        const GMP_RATS_mpq_struct *arg2)
     {
     //Eyeball the input parameters.
     assert(result != NULL);
     assert(arg1   != NULL);
     assert(arg2   != NULL);
     //If either rational number in the input is invalid,
     //flag the result as invalid.
     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
        {
        GMP_RATS_mpq_set_si(result, 1, 0);
        }
     else
        {
        //Rational number multiplication is a simple matter.
        //Just multiply components.  Don't need to worry
        //about rational numbers overlapping, as numerator
        //operations and denominator operations are separate.
        GMP_INTS_mpz_mul(&(result->num),
                         &(arg1->num),
                         &(arg2->num));
        GMP_INTS_mpz_mul(&(result->den),
                         &(arg1->den),
                         &(arg2->den));
        //Normalize it.
        GMP_RATS_mpq_normalize(result);
        }
     }
  //08/16/01: Visual inspection OK.
  void GMP_RATS_mpq_div(      GMP_RATS_mpq_struct *result,
                        const GMP_RATS_mpq_struct *arg1,
                        const GMP_RATS_mpq_struct *arg2)
     {
     GMP_RATS_mpq_struct rv;
     //Eyeball the input parameters.
     assert(result != NULL);
     assert(arg1   != NULL);
     assert(arg2   != NULL);
     //If either rational number in the input is invalid,
     //flag the result as invalid.
     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
        {
        GMP_RATS_mpq_set_si(result, 1, 0);
        }
     else
        {
        //Rational number division is a simple matter.
        //Just multiply components.  We do need to worry
        //about rational numbers overlapping, so must
        //make a copy of the return value.  If denominator
        //of return value is zero, it is NAN, but caller
        //should detect this.
        //
        //Allocate return value.
        GMP_RATS_mpq_init(&rv);
        //Calculate quotient.
        GMP_INTS_mpz_mul(&(rv.num),
                         &(arg1->num),
                         &(arg2->den));
        GMP_INTS_mpz_mul(&(rv.den),
                         &(arg1->den),
                         &(arg2->num));
        //Normalize quotient.
        GMP_RATS_mpq_normalize(&rv);
        //Copy to its destination.
        GMP_RATS_mpq_copy(result, &rv);
        //Deallocate temporary return value.
        GMP_RATS_mpq_clear(&rv);
        }
     }
  /******************************************************************/
  /***  COMPARISON FUNCTIONS  ***************************************/
  /******************************************************************/
  //08/16/01: Visual inspection OK.
  int GMP_RATS_mpq_cmp(const  GMP_RATS_mpq_struct *arg1,
                       const  GMP_RATS_mpq_struct *arg2,
                       int   *failure)
     {
     int arg1_sgn;
     int arg2_sgn;
     int rv, failure_rv;
     GMP_INTS_mpz_struct prod1, prod2;
     //Eyeball the input parameters.  Note that the third
     //parameter may be NULL.
     assert(arg1 != NULL);
     assert(arg2 != NULL);
     //If either of the input arguments are NAN, we
     //cannot compare arguments.  We return 0, and it
     //depends on the caller whether it is important
     //that the comparison is bogus.
     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
        {
        if (failure != NULL)
           *failure = 1;
        return(0);
        }
     //Calculate the sign of the left argument.  The encoding
     //we'll use is -1 means negative, 0 means zero, and
     //1 means positive.
     if (GMP_INTS_mpz_is_zero(&(arg1->num)))
        {
        arg1_sgn = 0;
        }
     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))
        {
        arg1_sgn = 1;
        }
     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))
        {
        arg1_sgn = -1;
        }
     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))
        {
        arg1_sgn = -1;
        }
     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))
        {
        arg1_sgn = 1;
        }
     //Calculate the sign of the right argument.  The encoding
     //we'll use is -1 means negative, 0 means zero, and
     //1 means positive.
     if (GMP_INTS_mpz_is_zero(&(arg2->num)))
        {
        arg2_sgn = 0;
        }
     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))
        {
        arg2_sgn = 1;
        }
     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))
        {
        arg2_sgn = -1;
        }
     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))
        {
        arg2_sgn = -1;
        }
     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))
        {
        arg2_sgn = 1;
        }
     //OK, can handle some simple cases where the signs of the
     //operands are different or both are zero.
     if ((arg1_sgn == 0) && (arg2_sgn == 0))
        {
        if (failure != NULL)
           *failure = 0;
        return(0);
        }
     else if ((arg1_sgn == -1) && (arg2_sgn > -1))
        {
        if (failure != NULL)
           *failure = 0;
        return(-1);
        }
     else if ((arg1_sgn == 0) && (arg2_sgn < 0))
        {
        if (failure != NULL)
           *failure = 0;
        return(1);
        }
     else if ((arg1_sgn == 0) && (arg2_sgn > 0))
        {
        if (failure != NULL)
           *failure = 0;
        return(-1);
        }
     else if ((arg1_sgn == 1) && (arg2_sgn < 1))
        {
        if (failure != NULL)
           *failure = 0;
        return(1);
        }
     //OK at this point, we cannot make a simple determination
     //as to the relative ordering.  The signs of arg1 and
     //arg2 are both the same, either both positive or both
     //negative.  We have to do a multiplication to sort
     //it out.
     //
     //Allocate the two integers to hold multiplication
     //results.
     GMP_INTS_mpz_init(&prod1);
     GMP_INTS_mpz_init(&prod2);
     //Cross-multiply to get relative magnitudes.
     GMP_INTS_mpz_mul(&prod1, &(arg1->num), &(arg2->den));
     GMP_INTS_mpz_mul(&prod2, &(arg1->den), &(arg2->num));
     //Take absolute values.
     GMP_INTS_mpz_abs(&prod1);
     GMP_INTS_mpz_abs(&prod2);
     //If we overflowed either multiplication and generated
     //a NAN, we cannot complete the compare.
     if (GMP_INTS_mpz_get_flags(&prod1) || GMP_INTS_mpz_get_flags(&prod2))
        {
        failure_rv = 1;
        rv = 0;
        }
     //If both rational numbers were effectively positive, we can
     //use the relative ordering of the products as the relative
     //ordering of the rational numbers.
     else if (arg1_sgn == 1)
        {
        //Compare the integers.
        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);
        //Clamp the return value.
        if (rv < 0)
           rv = -1;
        else if (rv == 0)
           rv = 0;
        else if (rv > 0)
           rv = 1;
        //There was no error.
        failure_rv = 0;
        }
     else
        {
        //The only case that *should* allow us to be here is
        //if the sign of both numbers is neg.
        assert(arg1_sgn == -1);
        //Compare the integers.
        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);
        //Invert and clamp the return value.
        if (rv < 0)
           rv = 1;
        else if (rv == 0)
           rv = 0;
        else if (rv > 0)
           rv = -1;
        //There was no error.
        failure_rv = 0;
        }
     //Deallocate the two integers.
     GMP_INTS_mpz_clear(&prod1);
     GMP_INTS_mpz_clear(&prod2);
     //Return the return values.
     if (failure != NULL)
        *failure = failure_rv;
     return(rv);
     }
  /******************************************************************/
  /***  VERSION CONTROL REPORTING FUNCTIONS  ************************/
  /******************************************************************/
  //08/07/01:  Visual inspection OK.
  const char *GMP_RATS_cvcinfo(void)
     {
     return("$Header$");
     }
  //08/07/01:  Visual inspection OK.
  const char *GMP_RATS_hvcinfo(void)
     {
     return(GMP_RATS_H_VERSION);
     }
  //End of gmp_rats.c.

 Legend:



Removed from v.56
 


changed lines


 
Added in v.71
 Legend:



Removed from v.56
 


changed lines


 
Added in v.71
-Removed from v.56
+Added in v.71

dashley@gmail.com	ViewVC Help
Powered by ViewVC 1.1.25