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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 70 by dashley, Sat Oct 29 01:53:01 2016 UTC revision 71 by dashley, Sat Nov 5 11:07:06 2016 UTC
# Line 1  Line 1 
1  //$Header$  //$Header$
2  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
3  //This file is part of "David T. Ashley's Shared Source Code", a set of shared components  //This file is part of "David T. Ashley's Shared Source Code", a set of shared components
4  //integrated into many of David T. Ashley's projects.  //integrated into many of David T. Ashley's projects.
5  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
6  //This source code and any program in which it is compiled/used is provided under the MIT License,  //This source code and any program in which it is compiled/used is provided under the MIT License,
7  //reproduced below.  //reproduced below.
8  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
9  //Permission is hereby granted, free of charge, to any person obtaining a copy of  //Permission is hereby granted, free of charge, to any person obtaining a copy of
10  //this software and associated documentation files(the "Software"), to deal in the  //this software and associated documentation files(the "Software"), to deal in the
11  //Software without restriction, including without limitation the rights to use,  //Software without restriction, including without limitation the rights to use,
12  //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the  //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
13  //Software, and to permit persons to whom the Software is furnished to do so,  //Software, and to permit persons to whom the Software is furnished to do so,
14  //subject to the following conditions :  //subject to the following conditions :
15  //  //
16  //The above copyright notice and this permission notice shall be included in all  //The above copyright notice and this permission notice shall be included in all
17  //copies or substantial portions of the Software.  //copies or substantial portions of the Software.
18  //  //
19  //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,  //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE  //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
22  //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER  //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  //SOFTWARE.  //SOFTWARE.
26  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
27  #define MODULE_GMP_RATS  #define MODULE_GMP_RATS
28    
29  #include <assert.h>  #include <assert.h>
30  #include <stdio.h>  #include <stdio.h>
31  #include <string.h>  #include <string.h>
32    
33  #include "bstrfunc.h"  #include "bstrfunc.h"
34  #include "charfunc.h"  #include "charfunc.h"
35  #include "gmp_ints.h"  #include "gmp_ints.h"
36  #include "gmp_rats.h"  #include "gmp_rats.h"
37    
38  #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)  #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
39     #include "ccmalloc.h"     #include "ccmalloc.h"
40  #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)  #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
41     #include "tclalloc.h"     #include "tclalloc.h"
42  #else  #else
43     #include <malloc.h>     #include <malloc.h>
44  #endif  #endif
45    
46    
47  /******************************************************************/  /******************************************************************/
48  /***  STATUS FUNCTIONS  *******************************************/  /***  STATUS FUNCTIONS  *******************************************/
49  /******************************************************************/  /******************************************************************/
50  //Functions in this category provide information about rational  //Functions in this category provide information about rational
51  //numbers.  //numbers.
52  //08/08/01:  Visual inspection OK.  //08/08/01:  Visual inspection OK.
53  int GMP_RATS_mpq_is_nan(const GMP_RATS_mpq_struct *rn)  int GMP_RATS_mpq_is_nan(const GMP_RATS_mpq_struct *rn)
54     {     {
55     assert(rn != NULL);     assert(rn != NULL);
56    
57     //A rational number is NAN in one of two     //A rational number is NAN in one of two
58     //circumstances.  If either of the integer components     //circumstances.  If either of the integer components
59     //is NAN, or else if there is a zero denominator.     //is NAN, or else if there is a zero denominator.
60     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
61        {        {
62        return(1);        return(1);
63        }        }
64     if (GMP_INTS_mpz_is_zero(&(rn->den)))     if (GMP_INTS_mpz_is_zero(&(rn->den)))
65        {        {
66        return(1);        return(1);
67        }        }
68    
69     //We're clean ...     //We're clean ...
70     return(0);     return(0);
71     }     }
72    
73    
74  /******************************************************************/  /******************************************************************/
75  /***  INITIALIZATION, CLEARING, AND SETTING FUNCTIONS  ************/  /***  INITIALIZATION, CLEARING, AND SETTING FUNCTIONS  ************/
76  /******************************************************************/  /******************************************************************/
77  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
78  void GMP_RATS_mpq_init(GMP_RATS_mpq_struct *arg)  void GMP_RATS_mpq_init(GMP_RATS_mpq_struct *arg)
79     {     {
80     //Eyeball the input parameter.     //Eyeball the input parameter.
81     assert(arg != NULL);     assert(arg != NULL);
82    
83     //Initialize the numerator and denominator.     //Initialize the numerator and denominator.
84     GMP_INTS_mpz_init(&(arg->num));     GMP_INTS_mpz_init(&(arg->num));
85     GMP_INTS_mpz_init(&(arg->den));     GMP_INTS_mpz_init(&(arg->den));
86    
87     //Canonically, we must start off as 0/1--canonical zero.     //Canonically, we must start off as 0/1--canonical zero.
88     GMP_INTS_mpz_set_ui(&(arg->num), 0);     GMP_INTS_mpz_set_ui(&(arg->num), 0);
89     GMP_INTS_mpz_set_ui(&(arg->den), 1);     GMP_INTS_mpz_set_ui(&(arg->den), 1);
90     }     }
91    
92    
93  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
94  void GMP_RATS_mpq_clear(GMP_RATS_mpq_struct *arg)  void GMP_RATS_mpq_clear(GMP_RATS_mpq_struct *arg)
95     {     {
96     //Eyeball the input parameter.     //Eyeball the input parameter.
97     assert(arg != NULL);     assert(arg != NULL);
98    
99     //Clear the numerator and denominator.  The called functions     //Clear the numerator and denominator.  The called functions
100     //will check for NULL pointers and so forth.     //will check for NULL pointers and so forth.
101     GMP_INTS_mpz_clear(&(arg->num));     GMP_INTS_mpz_clear(&(arg->num));
102     GMP_INTS_mpz_clear(&(arg->den));     GMP_INTS_mpz_clear(&(arg->den));
103     }     }
104    
105    
106  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
107  void GMP_RATS_mpq_set_si(GMP_RATS_mpq_struct *arg,  void GMP_RATS_mpq_set_si(GMP_RATS_mpq_struct *arg,
108                           int num,                           int num,
109                           int den)                           int den)
110     {     {
111     //Eyeball the input parameters.     //Eyeball the input parameters.
112     assert(arg != NULL);     assert(arg != NULL);
113    
114     //Set the numerator and denominator.     //Set the numerator and denominator.
115     GMP_INTS_mpz_set_si(&(arg->num), num);     GMP_INTS_mpz_set_si(&(arg->num), num);
116     GMP_INTS_mpz_set_si(&(arg->den), den);     GMP_INTS_mpz_set_si(&(arg->den), den);
117     }     }
118    
119    
120  //08/08/01:  Visual inspection OK.  //08/08/01:  Visual inspection OK.
121  void GMP_RATS_mpq_copy(      GMP_RATS_mpq_struct *dst,  void GMP_RATS_mpq_copy(      GMP_RATS_mpq_struct *dst,
122                         const GMP_RATS_mpq_struct *src)                         const GMP_RATS_mpq_struct *src)
123     {     {
124     assert(dst != NULL);     assert(dst != NULL);
125     assert(src != NULL);     assert(src != NULL);
126    
127     GMP_INTS_mpz_copy(&(dst->num), &(src->num));     GMP_INTS_mpz_copy(&(dst->num), &(src->num));
128     GMP_INTS_mpz_copy(&(dst->den), &(src->den));     GMP_INTS_mpz_copy(&(dst->den), &(src->den));
129     }     }
130    
131    
132  //08/13/01:  Visual inspection OK.  //08/13/01:  Visual inspection OK.
133  void GMP_RATS_mpq_swap( GMP_RATS_mpq_struct *a,  void GMP_RATS_mpq_swap( GMP_RATS_mpq_struct *a,
134                          GMP_RATS_mpq_struct *b)                          GMP_RATS_mpq_struct *b)
135     {     {
136     assert(a != NULL);     assert(a != NULL);
137     assert(b != NULL);     assert(b != NULL);
138    
139     //Handle the swap by swapping integer components.     //Handle the swap by swapping integer components.
140     GMP_INTS_mpz_swap(&(a->num), &(b->num));     GMP_INTS_mpz_swap(&(a->num), &(b->num));
141     GMP_INTS_mpz_swap(&(a->den), &(b->den));     GMP_INTS_mpz_swap(&(a->den), &(b->den));
142     }     }
143    
144    
145  //08/13/01:  Visual inspection OK.  //08/13/01:  Visual inspection OK.
146  void GMP_RATS_mpq_swap_components(GMP_RATS_mpq_struct *arg)  void GMP_RATS_mpq_swap_components(GMP_RATS_mpq_struct *arg)
147     {     {
148     assert(arg != NULL);     assert(arg != NULL);
149    
150     GMP_INTS_mpz_swap(&(arg->num), &(arg->den));     GMP_INTS_mpz_swap(&(arg->num), &(arg->den));
151     }     }
152    
153    
154  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
155  void GMP_RATS_mpq_set_complex_slash_sepd_rat_num(const char *s,  void GMP_RATS_mpq_set_complex_slash_sepd_rat_num(const char *s,
156                                                   int *failure,                                                   int *failure,
157                                                   GMP_RATS_mpq_struct *rn)                                                   GMP_RATS_mpq_struct *rn)
158     {     {
159     char *slash_posn, *numerator, *denominator;     char *slash_posn, *numerator, *denominator;
160     int s_len, numerator_len, denominator_len;     int s_len, numerator_len, denominator_len;
161     int i;     int i;
162    
163     //Eyeball the input parameters.     //Eyeball the input parameters.
164     assert(s != NULL);     assert(s != NULL);
165     assert(failure != NULL);     assert(failure != NULL);
166     assert(rn != NULL);     assert(rn != NULL);
167    
168     //Start off believing there is no failure.     //Start off believing there is no failure.
169     *failure = 0;     *failure = 0;
170    
171     //Figure out if there is one and only one slash in the     //Figure out if there is one and only one slash in the
172     //string.  If this condition isn't met, we cannot     //string.  If this condition isn't met, we cannot
173     //go further.     //go further.
174     slash_posn = strchr(s, '/');     slash_posn = strchr(s, '/');
175     if (!slash_posn)     if (!slash_posn)
176        {        {
177        *failure = 1;        *failure = 1;
178        return;        return;
179        }        }
180     if (strchr(slash_posn + 1, '/')) //There is a second occurence.     if (strchr(slash_posn + 1, '/')) //There is a second occurence.
181        {        {
182        *failure = 1;        *failure = 1;
183        return;        return;
184        }        }
185    
186     //At this point we have one and only one slash.     //At this point we have one and only one slash.
187     //Crack the string in two.  We must do this because the     //Crack the string in two.  We must do this because the
188     //input is a constant string.  We are not allowed to touch it     //input is a constant string.  We are not allowed to touch it
189     //in the logical domain because of the "const" keyword.  We can't     //in the logical domain because of the "const" keyword.  We can't
190     //do this in the physical domain because the debugger will nail     //do this in the physical domain because the debugger will nail
191     //us for it.     //us for it.
192     s_len = strlen(s);     s_len = strlen(s);
193     numerator_len = slash_posn - s;     numerator_len = slash_posn - s;
194     denominator_len = strlen(slash_posn + 1);     denominator_len = strlen(slash_posn + 1);
195     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
196        numerator   = CCMALLOC_malloc(sizeof(char) * (numerator_len   + 1));        numerator   = CCMALLOC_malloc(sizeof(char) * (numerator_len   + 1));
197        denominator = CCMALLOC_malloc(sizeof(char) * (denominator_len + 1));        denominator = CCMALLOC_malloc(sizeof(char) * (denominator_len + 1));
198     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
199        numerator   = TclpAlloc(sizeof(char) * (numerator_len   + 1));        numerator   = TclpAlloc(sizeof(char) * (numerator_len   + 1));
200        denominator = TclpAlloc(sizeof(char) * (denominator_len + 1));        denominator = TclpAlloc(sizeof(char) * (denominator_len + 1));
201     #else     #else
202        numerator   = malloc(sizeof(char) * (numerator_len   + 1));        numerator   = malloc(sizeof(char) * (numerator_len   + 1));
203        denominator = malloc(sizeof(char) * (denominator_len + 1));        denominator = malloc(sizeof(char) * (denominator_len + 1));
204     #endif     #endif
205    
206     assert(numerator   != NULL);     assert(numerator   != NULL);
207     assert(denominator != NULL);     assert(denominator != NULL);
208    
209     for (i=0; i<numerator_len; i++)     for (i=0; i<numerator_len; i++)
210        {        {
211        numerator[i] = s[i];        numerator[i] = s[i];
212        }          }  
213     numerator[numerator_len] = 0;     numerator[numerator_len] = 0;
214        
215     for (i=0; i<denominator_len; i++)     for (i=0; i<denominator_len; i++)
216        {        {
217        denominator[i] = s[slash_posn - s + 1 + i];        denominator[i] = s[slash_posn - s + 1 + i];
218        }        }
219     denominator[denominator_len] = 0;     denominator[denominator_len] = 0;
220    
221     //Try to parse out the numerator as an arbitrary integer.     //Try to parse out the numerator as an arbitrary integer.
222     //If this can't be done, it is an immediate failure.     //If this can't be done, it is an immediate failure.
223     GMP_INTS_mpz_set_general_int(&(rn->num),     GMP_INTS_mpz_set_general_int(&(rn->num),
224                                  failure,                                  failure,
225                                  numerator);                                  numerator);
226     if (*failure)     if (*failure)
227        {        {
228        *failure = 1;  //Clamp to 1, don't know what non-zero value        *failure = 1;  //Clamp to 1, don't know what non-zero value
229                       //was there.                       //was there.
230        goto ret_pt;        goto ret_pt;
231        }        }
232    
233     //Try to parse out the denominator.     //Try to parse out the denominator.
234     GMP_INTS_mpz_set_general_int(&(rn->den),     GMP_INTS_mpz_set_general_int(&(rn->den),
235                                  failure,                                  failure,
236                                  denominator);                                  denominator);
237     if (*failure)     if (*failure)
238        {        {
239        *failure = 1;  //Clamp to 1, don't know what non-zero value        *failure = 1;  //Clamp to 1, don't know what non-zero value
240                       //was there.                       //was there.
241        goto ret_pt;        goto ret_pt;
242        }        }
243    
244     //At this point, we have both a numerator and denominator.     //At this point, we have both a numerator and denominator.
245     //Clean up and return.     //Clean up and return.
246     ret_pt:     ret_pt:
247     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
248        CCMALLOC_free(numerator);        CCMALLOC_free(numerator);
249        CCMALLOC_free(denominator);        CCMALLOC_free(denominator);
250     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
251        TclpFree(numerator);        TclpFree(numerator);
252        TclpFree(denominator);        TclpFree(denominator);
253     #else     #else
254        free(numerator);        free(numerator);
255        free(denominator);        free(denominator);
256     #endif     #endif
257     }     }
258    
259    
260  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
261  void GMP_RATS_mpq_set_sci_not_rat_num(const char *s,  void GMP_RATS_mpq_set_sci_not_rat_num(const char *s,
262                                        int *failure,                                        int *failure,
263                                        GMP_RATS_mpq_struct *rn)                                        GMP_RATS_mpq_struct *rn)
264     {     {
265     int parse_failure;     int parse_failure;
266        //Return code from the floating point parsing        //Return code from the floating point parsing
267        //function.        //function.
268     char mant_sign;     char mant_sign;
269        //Sign character, if any, from the mantissa,        //Sign character, if any, from the mantissa,
270        //or N otherwise.        //or N otherwise.
271     size_t mant_bdp;     size_t mant_bdp;
272        //The index to the start of the mantissa before        //The index to the start of the mantissa before
273        //the decimal point.        //the decimal point.
274     size_t mant_bdp_len;     size_t mant_bdp_len;
275        //The length of the mantissa before the decimal        //The length of the mantissa before the decimal
276        //point.  Zero means not defined, i.e. that        //point.  Zero means not defined, i.e. that
277        //no characters were parsed and interpreted as        //no characters were parsed and interpreted as
278        //that part of a floating point number.        //that part of a floating point number.
279     size_t mant_adp;     size_t mant_adp;
280     size_t mant_adp_len;     size_t mant_adp_len;
281        //Similar fields for after the decimal point.        //Similar fields for after the decimal point.
282     char exp_sign;     char exp_sign;
283        //Sign of the exponent, if any, or N otherwise.        //Sign of the exponent, if any, or N otherwise.
284     size_t exp;     size_t exp;
285     size_t exp_len;     size_t exp_len;
286        //Similar fields as to the mantissa, but for the        //Similar fields as to the mantissa, but for the
287        //exponent.        //exponent.
288     size_t si;     size_t si;
289        //Iteration variable.        //Iteration variable.
290     int exponent_val;     int exponent_val;
291        //The value of the exponent.  We can't accept        //The value of the exponent.  We can't accept
292        //an exponent outside the range of a 24-bit        //an exponent outside the range of a 24-bit
293        //signed integer.  The 24-bit limit is arbitrary.        //signed integer.  The 24-bit limit is arbitrary.
294        //For one thing, it gives room to detect overflow        //For one thing, it gives room to detect overflow
295        //as are adding and multiplying by 10.        //as are adding and multiplying by 10.
296    
297     //Eyeball the input parameters.     //Eyeball the input parameters.
298     assert(s != NULL);     assert(s != NULL);
299     assert(failure != NULL);     assert(failure != NULL);
300     assert(rn != NULL);     assert(rn != NULL);
301        //Subcomponents of the rational number will be checked as        //Subcomponents of the rational number will be checked as
302        //we make integer calls, if we're in debug mode.        //we make integer calls, if we're in debug mode.
303    
304     //Start off believing no failure.     //Start off believing no failure.
305     *failure = 0;     *failure = 0;
306    
307     //Set the output to 0/1.  This is the default case for some     //Set the output to 0/1.  This is the default case for some
308     //steps below.     //steps below.
309     GMP_RATS_mpq_set_si(rn, 0, 1);     GMP_RATS_mpq_set_si(rn, 0, 1);
310    
311     //Attempt to parse the number as a general number     //Attempt to parse the number as a general number
312     //in scientific notation.     //in scientific notation.
313     BSTRFUNC_parse_gen_sci_not_num(s,     BSTRFUNC_parse_gen_sci_not_num(s,
314                                    &parse_failure,                                    &parse_failure,
315                                    &mant_sign,                                    &mant_sign,
316                                    &mant_bdp,                                    &mant_bdp,
317                                    &mant_bdp_len,                                    &mant_bdp_len,
318                                    &mant_adp,                                    &mant_adp,
319                                    &mant_adp_len,                                    &mant_adp_len,
320                                    &exp_sign,                                    &exp_sign,
321                                    &exp,                                    &exp,
322                                    &exp_len);                                    &exp_len);
323    
324     //If it wouldn't parse as a general number, can't go further.     //If it wouldn't parse as a general number, can't go further.
325     if (parse_failure)     if (parse_failure)
326        {        {
327        *failure = 1;        *failure = 1;
328        return;        return;
329        }        }
330     else     else
331        {        {
332        //The number parsed out.  The general strategy is to form a rational number        //The number parsed out.  The general strategy is to form a rational number
333        //consisting of the mantissa, with the decimal point shifted fully right, over        //consisting of the mantissa, with the decimal point shifted fully right, over
334        //a denominator of 1.  From there, we process the exponent and combine it with        //a denominator of 1.  From there, we process the exponent and combine it with
335        //the number of characters after the decimal point to form a virtual exponent.        //the number of characters after the decimal point to form a virtual exponent.
336        //If the exponent is positive, we multiply the numerator by the power of 10.        //If the exponent is positive, we multiply the numerator by the power of 10.
337        //If the exponent is negative, we multiply the denominator by that power of 10.        //If the exponent is negative, we multiply the denominator by that power of 10.
338    
339        //We want to trim the trailing zeros off of the portion of the mantissa after the        //We want to trim the trailing zeros off of the portion of the mantissa after the
340        //decimal point.  We only need to back up indices, no need to make copies, etc.        //decimal point.  We only need to back up indices, no need to make copies, etc.
341        //Note that it is possible that there are only zeros, in which case we'll end        //Note that it is possible that there are only zeros, in which case we'll end
342        //up with a length of zero.        //up with a length of zero.
343        while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))        while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))
344           mant_adp_len--;           mant_adp_len--;
345    
346        //Trim the leading zeros off of the portion of the mantissa before the        //Trim the leading zeros off of the portion of the mantissa before the
347        //decimal point.  Note that it is possible that there is only a zero,        //decimal point.  Note that it is possible that there is only a zero,
348        //so we may trim it down to nothing.        //so we may trim it down to nothing.
349        while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))        while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))
350           {           {
351           mant_bdp++;           mant_bdp++;
352           mant_bdp_len--;           mant_bdp_len--;
353           }           }
354    
355        //If we have only zeros in the mantissa, both before the        //If we have only zeros in the mantissa, both before the
356        //decimal point and after, then we return 0.        //decimal point and after, then we return 0.
357        if ((mant_bdp_len + mant_adp_len) == 0)        if ((mant_bdp_len + mant_adp_len) == 0)
358           {           {
359           *failure = 0;           *failure = 0;
360           return;           return;
361           }           }
362                
363        //Convert the numerator to an integer which represents the        //Convert the numerator to an integer which represents the
364        //part before the mantissa and the part after the mantissa        //part before the mantissa and the part after the mantissa
365        //concatenated as an integer.  We could call a function to do        //concatenated as an integer.  We could call a function to do
366        //this, but the function is not really any better in algorithm.        //this, but the function is not really any better in algorithm.
367        //We can do it ourselves.        //We can do it ourselves.
368        GMP_INTS_mpz_set_ui(&(rn->num), 0);        GMP_INTS_mpz_set_ui(&(rn->num), 0);
369        for (si = 0; si < mant_bdp_len; si++)        for (si = 0; si < mant_bdp_len; si++)
370           {           {
371           int val;           int val;
372    
373           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);
374           val = CHARFUNC_digit_to_val(s[mant_bdp + si]);           val = CHARFUNC_digit_to_val(s[mant_bdp + si]);
375           if (val >= 0)           if (val >= 0)
376              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);
377           }           }
378        for (si = 0; si < mant_adp_len; si++)        for (si = 0; si < mant_adp_len; si++)
379           {           {
380           int val;           int val;
381    
382           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);           GMP_INTS_mpz_mul_si(&(rn->num), &(rn->num), 10);
383           val = CHARFUNC_digit_to_val(s[mant_adp + si]);           val = CHARFUNC_digit_to_val(s[mant_adp + si]);
384           if (val >= 0)           if (val >= 0)
385              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);              GMP_INTS_mpz_add_ui(&(rn->num), &(rn->num), val);
386           }           }
387    
388        //The numerator should now have an integer which is        //The numerator should now have an integer which is
389        //The absolute value of the mantissa.  Process the possible        //The absolute value of the mantissa.  Process the possible
390        //sign.        //sign.
391        if (mant_sign == '-')        if (mant_sign == '-')
392           GMP_INTS_mpz_negate(&(rn->num));           GMP_INTS_mpz_negate(&(rn->num));
393    
394        //We now need to form a value from the exponent, if any.        //We now need to form a value from the exponent, if any.
395        //First, tackle the exponent.  Process the        //First, tackle the exponent.  Process the
396        //exponent into a signed integer.  We have to        //exponent into a signed integer.  We have to
397        //balk at anything outside of 24 bits.  The        //balk at anything outside of 24 bits.  The
398        //procedure used automatically handles        //procedure used automatically handles
399        //leading zeros correctly.        //leading zeros correctly.
400        exponent_val = 0;        exponent_val = 0;
401        for (si=exp; si<(exp+exp_len); si++)        for (si=exp; si<(exp+exp_len); si++)
402           {           {
403           int val;           int val;
404    
405           val = CHARFUNC_digit_to_val(s[si]);           val = CHARFUNC_digit_to_val(s[si]);
406    
407           assert(val >= 0 && val <= 9);           assert(val >= 0 && val <= 9);
408    
409           exponent_val *= 10;           exponent_val *= 10;
410           exponent_val += val;           exponent_val += val;
411    
412           if (((exp_sign=='-') && (exponent_val>8388608))           if (((exp_sign=='-') && (exponent_val>8388608))
413              ||              ||
414              ((exp_sign != '-') && (exponent_val>8388607)))              ((exp_sign != '-') && (exponent_val>8388607)))
415              {              {
416              *failure = 1;              *failure = 1;
417              return;              return;
418              }              }
419           }           }
420    
421        //If we're here, the exponent has been computed and        //If we're here, the exponent has been computed and
422        //is within 24 bits.  However, we need to adjust for        //is within 24 bits.  However, we need to adjust for
423        //the sign.        //the sign.
424        if (exp_sign == '-')        if (exp_sign == '-')
425           exponent_val = -exponent_val;           exponent_val = -exponent_val;
426    
427        //We need to adjust the exponent for the number of digits        //We need to adjust the exponent for the number of digits
428        //after the decimal point.        //after the decimal point.
429        exponent_val -= mant_adp_len;          exponent_val -= mant_adp_len;  
430                
431        //Again, clip for size.        //Again, clip for size.
432        if ((exponent_val < -8388608) || (exponent_val > 8388607))        if ((exponent_val < -8388608) || (exponent_val > 8388607))
433           {           {
434           *failure = 1;           *failure = 1;
435           return;           return;
436           }           }
437                    
438        //There are two cases to consider.  If the exponent        //There are two cases to consider.  If the exponent
439        //is positive, we need to multiply the numerator        //is positive, we need to multiply the numerator
440        //by 10 exponentiated to the power of the exponent.        //by 10 exponentiated to the power of the exponent.
441        //If the exponent is negative, we need to do the        //If the exponent is negative, we need to do the
442        //same thing to the denominator.  If the exponent        //same thing to the denominator.  If the exponent
443        //is negative, we don't need to do anything.        //is negative, we don't need to do anything.
444        if (exponent_val > 0)        if (exponent_val > 0)
445           {           {
446           GMP_INTS_mpz_struct k10, k10_exponentiated;           GMP_INTS_mpz_struct k10, k10_exponentiated;
447    
448           GMP_INTS_mpz_init(&k10);           GMP_INTS_mpz_init(&k10);
449           GMP_INTS_mpz_init(&k10_exponentiated);           GMP_INTS_mpz_init(&k10_exponentiated);
450    
451           GMP_INTS_mpz_set_ui(&k10, 10);           GMP_INTS_mpz_set_ui(&k10, 10);
452    
453           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, exponent_val);           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, exponent_val);
454    
455           GMP_INTS_mpz_mul(&(rn->num), &(rn->num), &k10_exponentiated);           GMP_INTS_mpz_mul(&(rn->num), &(rn->num), &k10_exponentiated);
456    
457           GMP_INTS_mpz_clear(&k10);           GMP_INTS_mpz_clear(&k10);
458           GMP_INTS_mpz_clear(&k10_exponentiated);           GMP_INTS_mpz_clear(&k10_exponentiated);
459    
460           *failure = 0;           *failure = 0;
461    
462           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
463              *failure = 1;              *failure = 1;
464    
465           return;           return;
466           }           }
467        else if (exponent_val < 0)        else if (exponent_val < 0)
468           {           {
469           GMP_INTS_mpz_struct k10, k10_exponentiated;           GMP_INTS_mpz_struct k10, k10_exponentiated;
470    
471           GMP_INTS_mpz_init(&k10);           GMP_INTS_mpz_init(&k10);
472           GMP_INTS_mpz_init(&k10_exponentiated);           GMP_INTS_mpz_init(&k10_exponentiated);
473    
474           GMP_INTS_mpz_set_ui(&k10, 10);           GMP_INTS_mpz_set_ui(&k10, 10);
475    
476           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, -exponent_val);           GMP_INTS_mpz_pow_ui(&k10_exponentiated, &k10, -exponent_val);
477    
478           GMP_INTS_mpz_mul(&(rn->den), &(rn->den), &k10_exponentiated);           GMP_INTS_mpz_mul(&(rn->den), &(rn->den), &k10_exponentiated);
479    
480           GMP_INTS_mpz_clear(&k10);           GMP_INTS_mpz_clear(&k10);
481           GMP_INTS_mpz_clear(&k10_exponentiated);           GMP_INTS_mpz_clear(&k10_exponentiated);
482    
483           *failure = 0;           *failure = 0;
484    
485           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))           if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
486              *failure = 1;              *failure = 1;
487    
488           return;           return;
489           }           }
490        }        }
491     }     }
492    
493    
494  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
495  void GMP_RATS_mpq_set_all_format_rat_num(const char *s,  void GMP_RATS_mpq_set_all_format_rat_num(const char *s,
496                                           int *failure,                                           int *failure,
497                                           GMP_RATS_mpq_struct *rn)                                           GMP_RATS_mpq_struct *rn)
498     {     {
499     //Eyeball the input parameters.     //Eyeball the input parameters.
500     assert(s != NULL);     assert(s != NULL);
501     assert(failure != NULL);     assert(failure != NULL);
502     assert(rn != NULL);     assert(rn != NULL);
503    
504     //Assume no failure.     //Assume no failure.
505     *failure = 0;     *failure = 0;
506    
507     //Try in order to parse as integers with slash then     //Try in order to parse as integers with slash then
508     //as number in scientific notation.     //as number in scientific notation.
509     GMP_RATS_mpq_set_complex_slash_sepd_rat_num(s,     GMP_RATS_mpq_set_complex_slash_sepd_rat_num(s,
510                                                 failure,                                                 failure,
511                                                 rn);                                                 rn);
512     if (!*failure)     if (!*failure)
513        return;        return;
514    
515     GMP_RATS_mpq_set_sci_not_rat_num(s,     GMP_RATS_mpq_set_sci_not_rat_num(s,
516                                      failure,                                      failure,
517                                      rn);                                      rn);
518    
519     if (*failure)     if (*failure)
520        *failure = 1;  //Clamp output.        *failure = 1;  //Clamp output.
521     }     }
522    
523    
524  /******************************************************************/  /******************************************************************/
525  /***  NORMALIZATION FUNCTIONS  ************************************/  /***  NORMALIZATION FUNCTIONS  ************************************/
526  /******************************************************************/  /******************************************************************/
527  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
528  void GMP_RATS_mpq_normalize_sign(GMP_RATS_mpq_struct *rn)  void GMP_RATS_mpq_normalize_sign(GMP_RATS_mpq_struct *rn)
529     {     {
530     //Eyeball the input.     //Eyeball the input.
531     assert(rn != NULL);     assert(rn != NULL);
532    
533     if (GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))     if (GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))
534        {        {
535        //Both negative, can negate both, this leaves both positive,        //Both negative, can negate both, this leaves both positive,
536        //which is the normalized form for a positive rational        //which is the normalized form for a positive rational
537        //number.        //number.
538        GMP_INTS_mpz_negate(&rn->num);        GMP_INTS_mpz_negate(&rn->num);
539        GMP_INTS_mpz_negate(&rn->den);        GMP_INTS_mpz_negate(&rn->den);
540        }        }
541     else if (!GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))     else if (!GMP_INTS_mpz_is_neg(&rn->num) && GMP_INTS_mpz_is_neg(&rn->den))
542        {        {
543        //Denominator neg, numerator non-neg, can negate both.  This        //Denominator neg, numerator non-neg, can negate both.  This
544        //will leave numerator neg, denominator pos, which is        //will leave numerator neg, denominator pos, which is
545        //normalized form for negative rational number.        //normalized form for negative rational number.
546        GMP_INTS_mpz_negate(&rn->num);        GMP_INTS_mpz_negate(&rn->num);
547        GMP_INTS_mpz_negate(&rn->den);        GMP_INTS_mpz_negate(&rn->den);
548        }            }    
549     }     }
550    
551    
552  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
553  void GMP_RATS_mpq_normalize(GMP_RATS_mpq_struct *rn)  void GMP_RATS_mpq_normalize(GMP_RATS_mpq_struct *rn)
554     {     {
555     //Eyeball the input.     //Eyeball the input.
556     assert(rn != NULL);     assert(rn != NULL);
557    
558     //Cover some special cases.  If either component has flags     //Cover some special cases.  If either component has flags
559     //set, don't even touch it.     //set, don't even touch it.
560     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))     if (GMP_INTS_mpz_get_flags(&(rn->num)) || GMP_INTS_mpz_get_flags(&(rn->den)))
561        {        {
562        return;        return;
563        }        }
564     //If the denominator is zero, normalize it to 1/0, the canonical     //If the denominator is zero, normalize it to 1/0, the canonical
565     //for for an illegal rational number.     //for for an illegal rational number.
566     else if (GMP_INTS_mpz_is_zero(&(rn->den)))     else if (GMP_INTS_mpz_is_zero(&(rn->den)))
567        {        {
568        GMP_RATS_mpq_set_si(rn, 1, 0);        GMP_RATS_mpq_set_si(rn, 1, 0);
569        return;        return;
570        }        }
571     //If the numerator is zero, convert the number to the canonical     //If the numerator is zero, convert the number to the canonical
572     //form for zero of 0/1.     //form for zero of 0/1.
573     else if (GMP_INTS_mpz_is_zero(&(rn->num)))     else if (GMP_INTS_mpz_is_zero(&(rn->num)))
574        {        {
575        GMP_RATS_mpq_set_si(rn, 0, 1);        GMP_RATS_mpq_set_si(rn, 0, 1);
576        return;        return;
577        }        }
578     else     else
579        {        {
580        int num_is_neg;        int num_is_neg;
581        int den_is_neg;        int den_is_neg;
582        GMP_INTS_mpz_struct gcd, quotient, remainder;        GMP_INTS_mpz_struct gcd, quotient, remainder;
583    
584        //Allocate space for the integers used.        //Allocate space for the integers used.
585        GMP_INTS_mpz_init(&gcd);        GMP_INTS_mpz_init(&gcd);
586        GMP_INTS_mpz_init(&quotient);        GMP_INTS_mpz_init(&quotient);
587        GMP_INTS_mpz_init(&remainder);        GMP_INTS_mpz_init(&remainder);
588    
589        //This is the most normal case, where we need to        //This is the most normal case, where we need to
590        //look at reducing the numerator and denominator.        //look at reducing the numerator and denominator.
591        //One way to do it would be to obtain the g.c.d.        //One way to do it would be to obtain the g.c.d.
592        //and divide this out, and this is the route        //and divide this out, and this is the route
593        //we'll take.  However, must grab out the sign.        //we'll take.  However, must grab out the sign.
594        if (GMP_INTS_mpz_is_neg(&(rn->num)))        if (GMP_INTS_mpz_is_neg(&(rn->num)))
595           {           {
596           num_is_neg = 1;           num_is_neg = 1;
597           GMP_INTS_mpz_negate(&(rn->num));           GMP_INTS_mpz_negate(&(rn->num));
598           }           }
599        else        else
600           {           {
601           num_is_neg = 0;           num_is_neg = 0;
602           }           }
603    
604        if (GMP_INTS_mpz_is_neg(&(rn->den)))        if (GMP_INTS_mpz_is_neg(&(rn->den)))
605           {           {
606           den_is_neg = 1;           den_is_neg = 1;
607           GMP_INTS_mpz_negate(&(rn->den));           GMP_INTS_mpz_negate(&(rn->den));
608           }           }
609        else        else
610           {           {
611           den_is_neg = 0;           den_is_neg = 0;
612           }           }
613    
614        //Calculate the GCD.        //Calculate the GCD.
615        GMP_INTS_mpz_gcd(&gcd, &(rn->num), &(rn->den));        GMP_INTS_mpz_gcd(&gcd, &(rn->num), &(rn->den));
616    
617        //Divide the numerator by the GCD and store it        //Divide the numerator by the GCD and store it
618        //back.        //back.
619        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,
620                             &(rn->num), &gcd);                             &(rn->num), &gcd);
621        GMP_INTS_mpz_copy(&(rn->num), &quotient);        GMP_INTS_mpz_copy(&(rn->num), &quotient);
622    
623        //Divide the denominator by the GCD and store it        //Divide the denominator by the GCD and store it
624        //back.        //back.
625        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,        GMP_INTS_mpz_tdiv_qr(&quotient, &remainder,
626                             &(rn->den), &gcd);                             &(rn->den), &gcd);
627        GMP_INTS_mpz_copy(&(rn->den), &quotient);        GMP_INTS_mpz_copy(&(rn->den), &quotient);
628    
629        //We now need to adjust the sign.  Both the        //We now need to adjust the sign.  Both the
630        //numerator and denominator are definitely        //numerator and denominator are definitely
631        //positive.  Need to make the numerator        //positive.  Need to make the numerator
632        //negative if either but not both of the        //negative if either but not both of the
633        //original signs were negative.        //original signs were negative.
634        if ((num_is_neg && !den_is_neg) || (!num_is_neg && den_is_neg))        if ((num_is_neg && !den_is_neg) || (!num_is_neg && den_is_neg))
635           {           {
636           GMP_INTS_mpz_negate(&(rn->num));           GMP_INTS_mpz_negate(&(rn->num));
637           }           }
638    
639        //Deallocate space for the integers used.        //Deallocate space for the integers used.
640        GMP_INTS_mpz_clear(&gcd);        GMP_INTS_mpz_clear(&gcd);
641        GMP_INTS_mpz_clear(&quotient);        GMP_INTS_mpz_clear(&quotient);
642        GMP_INTS_mpz_clear(&remainder);        GMP_INTS_mpz_clear(&remainder);
643    
644        return;        return;
645        }            }    
646     }     }
647    
648    
649  /******************************************************************/  /******************************************************************/
650  /***  ARITHMETIC FUNCTIONS  ***************************************/  /***  ARITHMETIC FUNCTIONS  ***************************************/
651  /******************************************************************/  /******************************************************************/
652  //08/08/01:  Visual inspection OK.  //08/08/01:  Visual inspection OK.
653  void GMP_RATS_mpq_add(      GMP_RATS_mpq_struct *result,  void GMP_RATS_mpq_add(      GMP_RATS_mpq_struct *result,
654                        const GMP_RATS_mpq_struct *arg1,                        const GMP_RATS_mpq_struct *arg1,
655                        const GMP_RATS_mpq_struct *arg2)                        const GMP_RATS_mpq_struct *arg2)
656     {     {
657     GMP_RATS_mpq_struct rv;     GMP_RATS_mpq_struct rv;
658     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
659    
660     //Eyeball the input parameters.     //Eyeball the input parameters.
661     assert(result != NULL);     assert(result != NULL);
662     assert(arg1   != NULL);     assert(arg1   != NULL);
663     assert(arg2   != NULL);     assert(arg2   != NULL);
664    
665     //Generally speaking, we do not want to require that     //Generally speaking, we do not want to require that
666     //the arguments and the result be distinct, as this is     //the arguments and the result be distinct, as this is
667     //too much of a restriction on the caller.  The approach     //too much of a restriction on the caller.  The approach
668     //taken, somewhat wasteful, is to allocate a place for     //taken, somewhat wasteful, is to allocate a place for
669     //the return value.     //the return value.
670     //     //
671     //For addition, if we are adding a/b and c/d, the     //For addition, if we are adding a/b and c/d, the
672     //result is necessarily algebraically     //result is necessarily algebraically
673     //(ad + cb)/bd.     //(ad + cb)/bd.
674     //     //
675     //If either rational number in the input is invalid,     //If either rational number in the input is invalid,
676     //flag the result as invalid.     //flag the result as invalid.
677     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
678        {        {
679        GMP_RATS_mpq_set_si(result, 1, 0);        GMP_RATS_mpq_set_si(result, 1, 0);
680        }        }
681     else     else
682        {        {
683        //Both rational numbers are OK.  Can simply stage the        //Both rational numbers are OK.  Can simply stage the
684        //result by the algebraic identity and then        //result by the algebraic identity and then
685        //normalize it.  Only need one temporary variable.        //normalize it.  Only need one temporary variable.
686        //        //
687        //Initialize the rational number that we will use to        //Initialize the rational number that we will use to
688        //hold return value in case it is the same as one        //hold return value in case it is the same as one
689        //or both of the arguments.        //or both of the arguments.
690        GMP_RATS_mpq_init(&rv);        GMP_RATS_mpq_init(&rv);
691    
692        //Initialize the temporary integer.        //Initialize the temporary integer.
693        GMP_INTS_mpz_init(&temp);        GMP_INTS_mpz_init(&temp);
694    
695        //numerator = a * d        //numerator = a * d
696        GMP_INTS_mpz_mul(&(rv.num), &(arg1->num), &(arg2->den));        GMP_INTS_mpz_mul(&(rv.num), &(arg1->num), &(arg2->den));
697    
698        //temp = c * b        //temp = c * b
699        GMP_INTS_mpz_mul(&temp, &(arg2->num), &(arg1->den));        GMP_INTS_mpz_mul(&temp, &(arg2->num), &(arg1->den));
700    
701        //numerator = a * d + c * b        //numerator = a * d + c * b
702        GMP_INTS_mpz_add(&(rv.num), &(rv.num), &temp);        GMP_INTS_mpz_add(&(rv.num), &(rv.num), &temp);
703    
704        //denominator = b * d        //denominator = b * d
705        GMP_INTS_mpz_mul(&(rv.den), &(arg1->den), &(arg2->den));        GMP_INTS_mpz_mul(&(rv.den), &(arg1->den), &(arg2->den));
706    
707        //Copy the temporary result to the actual return value.        //Copy the temporary result to the actual return value.
708        //Had to wait until now in case result was the same        //Had to wait until now in case result was the same
709        //as either or both args.        //as either or both args.
710        GMP_RATS_mpq_copy(result, &rv);        GMP_RATS_mpq_copy(result, &rv);
711    
712        //Normalize the result.        //Normalize the result.
713        GMP_RATS_mpq_normalize(result);        GMP_RATS_mpq_normalize(result);
714    
715        //Free dynamic memory.        //Free dynamic memory.
716        GMP_RATS_mpq_clear(&rv);        GMP_RATS_mpq_clear(&rv);
717        GMP_INTS_mpz_clear(&temp);        GMP_INTS_mpz_clear(&temp);
718        }        }
719     }     }
720    
721    
722  //08/08/01:  Visual inspection OK.  //08/08/01:  Visual inspection OK.
723  void GMP_RATS_mpq_sub(      GMP_RATS_mpq_struct *result,  void GMP_RATS_mpq_sub(      GMP_RATS_mpq_struct *result,
724                        const GMP_RATS_mpq_struct *arg1,                        const GMP_RATS_mpq_struct *arg1,
725                        const GMP_RATS_mpq_struct *arg2)                        const GMP_RATS_mpq_struct *arg2)
726     {     {
727     GMP_RATS_mpq_struct negated_arg_2;     GMP_RATS_mpq_struct negated_arg_2;
728    
729     //Eyeball the input parameters.     //Eyeball the input parameters.
730     assert(result != NULL);     assert(result != NULL);
731     assert(arg1   != NULL);     assert(arg1   != NULL);
732     assert(arg2   != NULL);     assert(arg2   != NULL);
733    
734     //For the subtract function, we could do it directly,     //For the subtract function, we could do it directly,
735     //but might as well just define it recursively     //but might as well just define it recursively
736     //in terms of add.  We can't modify the inputs,     //in terms of add.  We can't modify the inputs,
737     //so copy the second off and negate it.  All error     //so copy the second off and negate it.  All error
738     //flags and so forth will propagate automatically.     //flags and so forth will propagate automatically.
739     //     //
740     //Allocate space for the negated arg 2.     //Allocate space for the negated arg 2.
741     GMP_RATS_mpq_init(&negated_arg_2);     GMP_RATS_mpq_init(&negated_arg_2);
742    
743     //Copy from the original.     //Copy from the original.
744     GMP_RATS_mpq_copy(&negated_arg_2, arg2);     GMP_RATS_mpq_copy(&negated_arg_2, arg2);
745    
746     //Negate the copy.  Negating the numerator will     //Negate the copy.  Negating the numerator will
747     //do it.     //do it.
748     GMP_INTS_mpz_negate(&(negated_arg_2.num));     GMP_INTS_mpz_negate(&(negated_arg_2.num));
749    
750     //Make the add, which now is really a subtract.     //Make the add, which now is really a subtract.
751     GMP_RATS_mpq_add(result, arg1, &negated_arg_2);     GMP_RATS_mpq_add(result, arg1, &negated_arg_2);
752    
753     //Destroy the temporary variable.     //Destroy the temporary variable.
754     GMP_RATS_mpq_clear(&negated_arg_2);     GMP_RATS_mpq_clear(&negated_arg_2);
755     }     }
756    
757    
758  //08/16/01: Visual inspection OK.  //08/16/01: Visual inspection OK.
759  void GMP_RATS_mpq_mul(      GMP_RATS_mpq_struct *result,  void GMP_RATS_mpq_mul(      GMP_RATS_mpq_struct *result,
760                        const GMP_RATS_mpq_struct *arg1,                        const GMP_RATS_mpq_struct *arg1,
761                        const GMP_RATS_mpq_struct *arg2)                        const GMP_RATS_mpq_struct *arg2)
762     {     {
763     //Eyeball the input parameters.     //Eyeball the input parameters.
764     assert(result != NULL);     assert(result != NULL);
765     assert(arg1   != NULL);     assert(arg1   != NULL);
766     assert(arg2   != NULL);     assert(arg2   != NULL);
767    
768     //If either rational number in the input is invalid,     //If either rational number in the input is invalid,
769     //flag the result as invalid.     //flag the result as invalid.
770     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
771        {        {
772        GMP_RATS_mpq_set_si(result, 1, 0);        GMP_RATS_mpq_set_si(result, 1, 0);
773        }        }
774     else     else
775        {        {
776        //Rational number multiplication is a simple matter.        //Rational number multiplication is a simple matter.
777        //Just multiply components.  Don't need to worry        //Just multiply components.  Don't need to worry
778        //about rational numbers overlapping, as numerator        //about rational numbers overlapping, as numerator
779        //operations and denominator operations are separate.        //operations and denominator operations are separate.
780        GMP_INTS_mpz_mul(&(result->num),        GMP_INTS_mpz_mul(&(result->num),
781                         &(arg1->num),                         &(arg1->num),
782                         &(arg2->num));                         &(arg2->num));
783        GMP_INTS_mpz_mul(&(result->den),        GMP_INTS_mpz_mul(&(result->den),
784                         &(arg1->den),                         &(arg1->den),
785                         &(arg2->den));                         &(arg2->den));
786    
787        //Normalize it.        //Normalize it.
788        GMP_RATS_mpq_normalize(result);        GMP_RATS_mpq_normalize(result);
789        }        }
790     }     }
791    
792    
793  //08/16/01: Visual inspection OK.  //08/16/01: Visual inspection OK.
794  void GMP_RATS_mpq_div(      GMP_RATS_mpq_struct *result,  void GMP_RATS_mpq_div(      GMP_RATS_mpq_struct *result,
795                        const GMP_RATS_mpq_struct *arg1,                        const GMP_RATS_mpq_struct *arg1,
796                        const GMP_RATS_mpq_struct *arg2)                        const GMP_RATS_mpq_struct *arg2)
797     {     {
798     GMP_RATS_mpq_struct rv;     GMP_RATS_mpq_struct rv;
799    
800     //Eyeball the input parameters.     //Eyeball the input parameters.
801     assert(result != NULL);     assert(result != NULL);
802     assert(arg1   != NULL);     assert(arg1   != NULL);
803     assert(arg2   != NULL);     assert(arg2   != NULL);
804    
805     //If either rational number in the input is invalid,     //If either rational number in the input is invalid,
806     //flag the result as invalid.     //flag the result as invalid.
807     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
808        {        {
809        GMP_RATS_mpq_set_si(result, 1, 0);        GMP_RATS_mpq_set_si(result, 1, 0);
810        }        }
811     else     else
812        {        {
813        //Rational number division is a simple matter.        //Rational number division is a simple matter.
814        //Just multiply components.  We do need to worry        //Just multiply components.  We do need to worry
815        //about rational numbers overlapping, so must        //about rational numbers overlapping, so must
816        //make a copy of the return value.  If denominator        //make a copy of the return value.  If denominator
817        //of return value is zero, it is NAN, but caller        //of return value is zero, it is NAN, but caller
818        //should detect this.        //should detect this.
819        //        //
820        //Allocate return value.        //Allocate return value.
821        GMP_RATS_mpq_init(&rv);        GMP_RATS_mpq_init(&rv);
822    
823        //Calculate quotient.        //Calculate quotient.
824        GMP_INTS_mpz_mul(&(rv.num),        GMP_INTS_mpz_mul(&(rv.num),
825                         &(arg1->num),                         &(arg1->num),
826                         &(arg2->den));                         &(arg2->den));
827        GMP_INTS_mpz_mul(&(rv.den),        GMP_INTS_mpz_mul(&(rv.den),
828                         &(arg1->den),                         &(arg1->den),
829                         &(arg2->num));                         &(arg2->num));
830    
831        //Normalize quotient.        //Normalize quotient.
832        GMP_RATS_mpq_normalize(&rv);        GMP_RATS_mpq_normalize(&rv);
833    
834        //Copy to its destination.        //Copy to its destination.
835        GMP_RATS_mpq_copy(result, &rv);        GMP_RATS_mpq_copy(result, &rv);
836    
837        //Deallocate temporary return value.        //Deallocate temporary return value.
838        GMP_RATS_mpq_clear(&rv);        GMP_RATS_mpq_clear(&rv);
839        }        }
840     }     }
841    
842    
843  /******************************************************************/  /******************************************************************/
844  /***  COMPARISON FUNCTIONS  ***************************************/  /***  COMPARISON FUNCTIONS  ***************************************/
845  /******************************************************************/  /******************************************************************/
846  //08/16/01: Visual inspection OK.  //08/16/01: Visual inspection OK.
847  int GMP_RATS_mpq_cmp(const  GMP_RATS_mpq_struct *arg1,  int GMP_RATS_mpq_cmp(const  GMP_RATS_mpq_struct *arg1,
848                       const  GMP_RATS_mpq_struct *arg2,                       const  GMP_RATS_mpq_struct *arg2,
849                       int   *failure)                       int   *failure)
850     {     {
851     int arg1_sgn;     int arg1_sgn;
852     int arg2_sgn;     int arg2_sgn;
853     int rv, failure_rv;     int rv, failure_rv;
854     GMP_INTS_mpz_struct prod1, prod2;     GMP_INTS_mpz_struct prod1, prod2;
855    
856     //Eyeball the input parameters.  Note that the third     //Eyeball the input parameters.  Note that the third
857     //parameter may be NULL.     //parameter may be NULL.
858     assert(arg1 != NULL);     assert(arg1 != NULL);
859     assert(arg2 != NULL);     assert(arg2 != NULL);
860    
861     //If either of the input arguments are NAN, we     //If either of the input arguments are NAN, we
862     //cannot compare arguments.  We return 0, and it     //cannot compare arguments.  We return 0, and it
863     //depends on the caller whether it is important     //depends on the caller whether it is important
864     //that the comparison is bogus.     //that the comparison is bogus.
865     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))     if (GMP_RATS_mpq_is_nan(arg1) || GMP_RATS_mpq_is_nan(arg2))
866        {        {
867        if (failure != NULL)        if (failure != NULL)
868           *failure = 1;           *failure = 1;
869        return(0);        return(0);
870        }        }
871    
872     //Calculate the sign of the left argument.  The encoding     //Calculate the sign of the left argument.  The encoding
873     //we'll use is -1 means negative, 0 means zero, and     //we'll use is -1 means negative, 0 means zero, and
874     //1 means positive.     //1 means positive.
875     if (GMP_INTS_mpz_is_zero(&(arg1->num)))     if (GMP_INTS_mpz_is_zero(&(arg1->num)))
876        {        {
877        arg1_sgn = 0;        arg1_sgn = 0;
878        }        }
879     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))
880        {        {
881        arg1_sgn = 1;        arg1_sgn = 1;
882        }        }
883     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))     else if (GMP_INTS_mpz_is_neg(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))
884        {        {
885        arg1_sgn = -1;        arg1_sgn = -1;
886        }        }
887     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_neg(&(arg1->den)))
888        {        {
889        arg1_sgn = -1;        arg1_sgn = -1;
890        }        }
891     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))     else if (GMP_INTS_mpz_is_pos(&(arg1->num)) && GMP_INTS_mpz_is_pos(&(arg1->den)))
892        {        {
893        arg1_sgn = 1;        arg1_sgn = 1;
894        }        }
895    
896     //Calculate the sign of the right argument.  The encoding     //Calculate the sign of the right argument.  The encoding
897     //we'll use is -1 means negative, 0 means zero, and     //we'll use is -1 means negative, 0 means zero, and
898     //1 means positive.     //1 means positive.
899     if (GMP_INTS_mpz_is_zero(&(arg2->num)))     if (GMP_INTS_mpz_is_zero(&(arg2->num)))
900        {        {
901        arg2_sgn = 0;        arg2_sgn = 0;
902        }        }
903     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))
904        {        {
905        arg2_sgn = 1;        arg2_sgn = 1;
906        }        }
907     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))     else if (GMP_INTS_mpz_is_neg(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))
908        {        {
909        arg2_sgn = -1;        arg2_sgn = -1;
910        }        }
911     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_neg(&(arg2->den)))
912        {        {
913        arg2_sgn = -1;        arg2_sgn = -1;
914        }        }
915     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))     else if (GMP_INTS_mpz_is_pos(&(arg2->num)) && GMP_INTS_mpz_is_pos(&(arg2->den)))
916        {        {
917        arg2_sgn = 1;        arg2_sgn = 1;
918        }        }
919    
920     //OK, can handle some simple cases where the signs of the     //OK, can handle some simple cases where the signs of the
921     //operands are different or both are zero.     //operands are different or both are zero.
922     if ((arg1_sgn == 0) && (arg2_sgn == 0))     if ((arg1_sgn == 0) && (arg2_sgn == 0))
923        {        {
924        if (failure != NULL)        if (failure != NULL)
925           *failure = 0;           *failure = 0;
926        return(0);        return(0);
927        }        }
928     else if ((arg1_sgn == -1) && (arg2_sgn > -1))     else if ((arg1_sgn == -1) && (arg2_sgn > -1))
929        {        {
930        if (failure != NULL)        if (failure != NULL)
931           *failure = 0;           *failure = 0;
932        return(-1);        return(-1);
933        }        }
934     else if ((arg1_sgn == 0) && (arg2_sgn < 0))     else if ((arg1_sgn == 0) && (arg2_sgn < 0))
935        {        {
936        if (failure != NULL)        if (failure != NULL)
937           *failure = 0;           *failure = 0;
938        return(1);        return(1);
939        }        }
940     else if ((arg1_sgn == 0) && (arg2_sgn > 0))     else if ((arg1_sgn == 0) && (arg2_sgn > 0))
941        {        {
942        if (failure != NULL)        if (failure != NULL)
943           *failure = 0;           *failure = 0;
944        return(-1);        return(-1);
945        }        }
946     else if ((arg1_sgn == 1) && (arg2_sgn < 1))     else if ((arg1_sgn == 1) && (arg2_sgn < 1))
947        {        {
948        if (failure != NULL)        if (failure != NULL)
949           *failure = 0;           *failure = 0;
950        return(1);        return(1);
951        }        }
952    
953     //OK at this point, we cannot make a simple determination     //OK at this point, we cannot make a simple determination
954     //as to the relative ordering.  The signs of arg1 and     //as to the relative ordering.  The signs of arg1 and
955     //arg2 are both the same, either both positive or both     //arg2 are both the same, either both positive or both
956     //negative.  We have to do a multiplication to sort     //negative.  We have to do a multiplication to sort
957     //it out.     //it out.
958     //     //
959     //Allocate the two integers to hold multiplication     //Allocate the two integers to hold multiplication
960     //results.     //results.
961     GMP_INTS_mpz_init(&prod1);     GMP_INTS_mpz_init(&prod1);
962     GMP_INTS_mpz_init(&prod2);     GMP_INTS_mpz_init(&prod2);
963    
964     //Cross-multiply to get relative magnitudes.     //Cross-multiply to get relative magnitudes.
965     GMP_INTS_mpz_mul(&prod1, &(arg1->num), &(arg2->den));     GMP_INTS_mpz_mul(&prod1, &(arg1->num), &(arg2->den));
966     GMP_INTS_mpz_mul(&prod2, &(arg1->den), &(arg2->num));     GMP_INTS_mpz_mul(&prod2, &(arg1->den), &(arg2->num));
967        
968     //Take absolute values.     //Take absolute values.
969     GMP_INTS_mpz_abs(&prod1);     GMP_INTS_mpz_abs(&prod1);
970     GMP_INTS_mpz_abs(&prod2);     GMP_INTS_mpz_abs(&prod2);
971    
972     //If we overflowed either multiplication and generated     //If we overflowed either multiplication and generated
973     //a NAN, we cannot complete the compare.     //a NAN, we cannot complete the compare.
974     if (GMP_INTS_mpz_get_flags(&prod1) || GMP_INTS_mpz_get_flags(&prod2))     if (GMP_INTS_mpz_get_flags(&prod1) || GMP_INTS_mpz_get_flags(&prod2))
975        {        {
976        failure_rv = 1;        failure_rv = 1;
977        rv = 0;        rv = 0;
978        }        }
979     //If both rational numbers were effectively positive, we can     //If both rational numbers were effectively positive, we can
980     //use the relative ordering of the products as the relative     //use the relative ordering of the products as the relative
981     //ordering of the rational numbers.     //ordering of the rational numbers.
982     else if (arg1_sgn == 1)     else if (arg1_sgn == 1)
983        {        {
984        //Compare the integers.        //Compare the integers.
985        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);
986    
987        //Clamp the return value.        //Clamp the return value.
988        if (rv < 0)        if (rv < 0)
989           rv = -1;           rv = -1;
990        else if (rv == 0)        else if (rv == 0)
991           rv = 0;           rv = 0;
992        else if (rv > 0)        else if (rv > 0)
993           rv = 1;           rv = 1;
994    
995        //There was no error.        //There was no error.
996        failure_rv = 0;        failure_rv = 0;
997        }        }
998     else     else
999        {        {
1000        //The only case that *should* allow us to be here is        //The only case that *should* allow us to be here is
1001        //if the sign of both numbers is neg.        //if the sign of both numbers is neg.
1002        assert(arg1_sgn == -1);        assert(arg1_sgn == -1);
1003    
1004        //Compare the integers.        //Compare the integers.
1005        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);        rv = GMP_INTS_mpz_cmp(&prod1, &prod2);
1006    
1007        //Invert and clamp the return value.        //Invert and clamp the return value.
1008        if (rv < 0)        if (rv < 0)
1009           rv = 1;           rv = 1;
1010        else if (rv == 0)        else if (rv == 0)
1011           rv = 0;           rv = 0;
1012        else if (rv > 0)        else if (rv > 0)
1013           rv = -1;           rv = -1;
1014    
1015        //There was no error.        //There was no error.
1016        failure_rv = 0;        failure_rv = 0;
1017        }        }
1018    
1019     //Deallocate the two integers.     //Deallocate the two integers.
1020     GMP_INTS_mpz_clear(&prod1);     GMP_INTS_mpz_clear(&prod1);
1021     GMP_INTS_mpz_clear(&prod2);     GMP_INTS_mpz_clear(&prod2);
1022    
1023     //Return the return values.     //Return the return values.
1024     if (failure != NULL)     if (failure != NULL)
1025        *failure = failure_rv;        *failure = failure_rv;
1026     return(rv);     return(rv);
1027     }     }
1028    
1029    
1030  /******************************************************************/  /******************************************************************/
1031  /***  VERSION CONTROL REPORTING FUNCTIONS  ************************/  /***  VERSION CONTROL REPORTING FUNCTIONS  ************************/
1032  /******************************************************************/  /******************************************************************/
1033  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
1034  const char *GMP_RATS_cvcinfo(void)  const char *GMP_RATS_cvcinfo(void)
1035     {     {
1036     return("$Header$");     return("$Header$");
1037     }     }
1038    
1039    
1040  //08/07/01:  Visual inspection OK.  //08/07/01:  Visual inspection OK.
1041  const char *GMP_RATS_hvcinfo(void)  const char *GMP_RATS_hvcinfo(void)
1042     {     {
1043     return(GMP_RATS_H_VERSION);     return(GMP_RATS_H_VERSION);
1044     }     }
1045    
1046  //End of gmp_rats.c.  //End of gmp_rats.c.

Legend:
Removed from v.70  
changed lines
  Added in v.71

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25