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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25