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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 56 - (show annotations) (download)
Sat Oct 29 01:53:01 2016 UTC (6 years, 3 months ago) by dashley
File MIME type: text/plain
File size: 33222 byte(s)
License and property (keyword) changes.
1 //$Header$
2 //-------------------------------------------------------------------------------------------------
3 //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 //-------------------------------------------------------------------------------------------------
6 //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 //
16 //The above copyright notice and this permission notice shall be included in all
17 //copies or substantial portions of the Software.
18 //
19 //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 //-------------------------------------------------------------------------------------------------
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 return("$Header$");
1037 }
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 //End of gmp_rats.c.

Properties

Name Value
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25