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("ient); |
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("ient, &remainder, |
620 |
&(rn->num), &gcd); |
621 |
GMP_INTS_mpz_copy(&(rn->num), "ient); |
622 |
|
623 |
//Divide the denominator by the GCD and store it |
624 |
//back. |
625 |
GMP_INTS_mpz_tdiv_qr("ient, &remainder, |
626 |
&(rn->den), &gcd); |
627 |
GMP_INTS_mpz_copy(&(rn->den), "ient); |
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("ient); |
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. |