1 /*
2 * Copyright © 1988-2004 Keith Packard and Bart Massey.
3 * All Rights Reserved. See the file COPYING in this directory
4 * for licensing information.
5 */
6
7 /*
8 * rational.c
9 *
10 * operationalns on rationals
11 */
12
13 #include "nickle.h"
14 #include <math.h>
15
16 int
RationalInit(void)17 RationalInit (void)
18 {
19 return 1;
20 }
21
22 #if 0
23 static Value
24 natural_to_rational (Natural *n)
25 {
26 ENTER ();
27 RETURN (NewRational (Positive, n, one_natural));
28 }
29 #endif
30
31 static Value
RationalPlusHelper(Sign sign,Rational * a,Rational * b)32 RationalPlusHelper (Sign sign, Rational *a, Rational *b)
33 {
34 ENTER ();
35 RETURN (NewRational (sign,
36 NaturalPlus (NaturalTimes (a->num, b->den),
37 NaturalTimes (b->num, a->den)),
38 NaturalTimes (a->den, b->den)));
39 }
40
41 static Value
RationalMinusHelper(Rational * a,Rational * b)42 RationalMinusHelper (Rational *a, Rational *b)
43 {
44 ENTER ();
45 Natural *ra, *rb, *t;
46 Sign sign = Positive;
47
48 ra = NaturalTimes (a->num, b->den);
49 rb = NaturalTimes (b->num, a->den);
50 if (NaturalLess (ra, rb))
51 {
52 sign = Negative;
53 t = ra;
54 ra = rb;
55 rb = t;
56 }
57 RETURN (NewRational (sign, NaturalMinus (ra, rb),
58 NaturalTimes (a->den, b->den)));
59 }
60
61 static Value
RationalPlus(Value av,Value bv,int expandOk)62 RationalPlus (Value av, Value bv, int expandOk)
63 {
64 ENTER ();
65 Rational *a = &av->rational, *b = &bv->rational;
66 Value ret;
67
68 switch (catagorize_signs(a->sign, b->sign)) {
69 case BothPositive:
70 case BothNegative:
71 ret = RationalPlusHelper (a->sign, a, b);
72 break;
73 case FirstPositive:
74 ret = RationalMinusHelper (a, b);
75 break;
76 case SecondPositive:
77 ret = RationalMinusHelper (b, a);
78 break;
79 default:
80 abort();
81 }
82 RETURN (ret);
83 }
84
85 static Value
RationalMinus(Value av,Value bv,int expandOk)86 RationalMinus (Value av, Value bv, int expandOk)
87 {
88 ENTER ();
89 Rational *a = &av->rational, *b = &bv->rational;
90 Value ret;
91
92 switch (catagorize_signs(a->sign, b->sign)) {
93 case BothPositive:
94 ret = RationalMinusHelper (a, b);
95 break;
96 case FirstPositive:
97 case SecondPositive:
98 ret = RationalPlusHelper (a->sign, a, b);
99 break;
100 case BothNegative:
101 ret = RationalMinusHelper (b, a);
102 break;
103 default:
104 abort();
105 }
106 RETURN (ret);
107 }
108
109 static Value
RationalTimes(Value av,Value bv,int expandOk)110 RationalTimes (Value av, Value bv, int expandOk)
111 {
112 ENTER ();
113 Rational *a = &av->rational, *b = &bv->rational;
114 Sign sign;
115
116 sign = Positive;
117 if (a->sign != b->sign)
118 sign = Negative;
119 RETURN (NewRational (sign,
120 NaturalTimes (a->num, b->num),
121 NaturalTimes (a->den, b->den)));
122 }
123
124 static Value
RationalDivide(Value av,Value bv,int expandOk)125 RationalDivide (Value av, Value bv, int expandOk)
126 {
127 ENTER ();
128 Rational *a = &av->rational, *b = &bv->rational;
129 Sign sign;
130
131 if (NaturalZero (b->num))
132 {
133 RaiseStandardException (exception_divide_by_zero, 2,
134 av, bv);
135 RETURN (Void);
136 }
137 sign = Positive;
138 if (a->sign != b->sign)
139 sign = Negative;
140 RETURN (NewRational (sign,
141 NaturalTimes (a->num, b->den),
142 NaturalTimes (a->den, b->num)));
143 }
144
145 /*
146 * Modulus for rational values.
147 *
148 * Sorta like for integers:
149 *
150 * c/d * (a/b | c/d) + a/b % c/d = a/b
151 *
152 * 0 <= a/b % c/d < abs (c/d)
153 * a/b | c/d is an integer
154 *
155 * To calculate modulus (e/f):
156 *
157 * c/d * n + e/f = a/b
158 * e/f = a/b - c/d * n
159 * (e * b * d) / f = a * d - c * b * n
160 *
161 * therefore (e * b * d) / f is integer
162 *
163 * c * b * n + (e * b * d) / f = a * d
164 * (e * b * d) / f = (a * d) % (c * b)
165 * e / f = ((a * d) % (c * b)) / (b * d)
166 */
167
168 static Value
RationalMod(Value av,Value bv,int expandOk)169 RationalMod (Value av, Value bv, int expandOk)
170 {
171 ENTER ();
172 Rational *a = &av->rational, *b = &bv->rational;
173 Natural *rem, *div;
174
175 if (NaturalZero (b->num))
176 {
177 RaiseStandardException (exception_divide_by_zero, 2,
178 av, bv);
179 RETURN (Void);
180 }
181 div = NaturalTimes (b->num, a->den);
182 (void) NaturalDivide (NaturalTimes (a->num, b->den), div, &rem);
183 if (a->sign == Negative && !NaturalZero (rem))
184 rem = NaturalMinus (div, rem);
185 RETURN (NewRational (Positive, rem, NaturalTimes (a->den, b->den)));
186 }
187
188
189 static Value
RationalLess(Value av,Value bv,int expandOk)190 RationalLess (Value av, Value bv, int expandOk)
191 {
192 ENTER ();
193 Rational *a = &av->rational, *b = &bv->rational;
194 Rational *t;
195 int ret;
196
197 switch (catagorize_signs (a->sign, b->sign)) {
198 case BothNegative:
199 t = a;
200 a = b;
201 b = t;
202 case BothPositive:
203 if (!NaturalEqual (a->den, b->den))
204 ret = NaturalLess (NaturalTimes (a->num, b->den),
205 NaturalTimes (b->num, a->den));
206 else
207 ret = NaturalLess (a->num, b->num);
208 break;
209 case FirstPositive:
210 ret = 0;
211 break;
212 case SecondPositive:
213 ret = 1;
214 break;
215 default:
216 abort();
217 }
218 RETURN (ret ? TrueVal : FalseVal);
219 }
220
221 static Value
RationalEqual(Value av,Value bv,int expandOk)222 RationalEqual (Value av, Value bv, int expandOk)
223 {
224 Rational *a = &av->rational, *b = &bv->rational;
225
226 if (a->sign == b->sign &&
227 NaturalEqual (a->num, b->num) &&
228 NaturalEqual (a->den, b->den))
229 {
230 return TrueVal;
231 }
232 return FalseVal;
233 }
234
235 static Value
RationalNegate(Value av,int expandOk)236 RationalNegate (Value av, int expandOk)
237 {
238 ENTER ();
239 Rational *a = &av->rational;
240
241 RETURN (NewRational (SignNegate (a->sign), a->num, a->den));
242 }
243
244 static Value
RationalFloor(Value av,int expandOk)245 RationalFloor (Value av, int expandOk)
246 {
247 ENTER ();
248 Rational *a = &av->rational;
249 Natural *quo, *rem;
250
251 quo = NaturalDivide (a->num, a->den, &rem);
252 if (!NaturalZero (rem) && a->sign == Negative)
253 quo = NaturalPlus (quo, one_natural);
254 RETURN (NewInteger (a->sign, quo));
255 }
256
257 static Value
RationalCeil(Value av,int expandOk)258 RationalCeil (Value av, int expandOk)
259 {
260 ENTER ();
261 Rational *a = &av->rational;
262 Natural *quo, *rem;
263
264 quo = NaturalDivide (a->num, a->den, &rem);
265 if (!NaturalZero (rem) && a->sign == Positive)
266 quo = NaturalPlus (quo, one_natural);
267 RETURN (NewInteger (a->sign, quo));
268 }
269
270 static Value
RationalPromote(Value av,Value bv)271 RationalPromote (Value av, Value bv)
272 {
273 ENTER ();
274
275 switch (ValueTag(av)) {
276 case rep_int:
277 av = NewIntRational (ValueInt(av));
278 break;
279 case rep_integer:
280 av = NewIntegerRational (&av->integer);
281 break;
282 default:
283 break;
284 }
285 RETURN (av);
286 }
287
288 static Value
RationalReduce(Value av)289 RationalReduce (Value av)
290 {
291 ENTER ();
292 Rational *a = &av->rational;
293
294 if (NaturalEqual (a->den, one_natural))
295 av = Reduce (NewInteger (a->sign, a->num));
296 RETURN (av);
297 }
298
299 static HashValue
RationalHash(Value av)300 RationalHash (Value av)
301 {
302 Rational *a = &av->rational;
303
304 return NaturalHash (a->den) ^ NaturalHash(a->num) ^ a->sign;
305 }
306
307 extern ValueRep IntegerRep;
308
309 extern Natural *NaturalFactor (Natural *, Natural *);
310 extern Natural *NaturalSqrt (Natural *);
311 extern Natural *NaturalIntPow (Natural *, int);
312 extern Natural *NaturalPow (Natural *, Natural *);
313 extern Natural *NaturalPowMod (Natural *, Natural *, Natural *);
314 extern Natural *two_natural;
315
316 static Natural *
NaturalPsi(Natural * a,Natural * max)317 NaturalPsi(Natural *a, Natural *max)
318 {
319 ENTER ();
320 Natural *p;
321 int n;
322 Natural *ret;
323 Natural *rem;
324 Natural *next;
325 Natural *pow;
326 Natural *fact;
327
328 ret = one_natural;
329 while (!NaturalEqual (a, one_natural))
330 {
331 p = NaturalFactor (a, max);
332 if (!p)
333 {
334 ret = 0;
335 break;
336 }
337 n = 0;
338 for (;;)
339 {
340 next = NaturalDivide (a, p, &rem);
341 if (!NaturalZero (rem))
342 break;
343 a = next;
344 n++;
345 }
346 pow = NaturalIntPow (p, n-1);
347 fact = NaturalMinus (NaturalTimes (pow, p), pow);
348 ret = NaturalTimes (ret, fact);
349 if (max && NaturalLess (max, fact))
350 break;
351 }
352 RETURN (ret);
353 }
354
355 #if 0
356 static int
357 IntSqrt (int a)
358 {
359 int l, h, m;
360
361 l = 2;
362 h = a/2;
363 while ((h-l) > 1)
364 {
365 m = (h+l) >> 1;
366 if (m * m < a)
367 l = m;
368 else
369 h = m;
370 }
371 return h;
372 }
373
374 static int
375 IntFactor (int a)
376 {
377 int v, lim;
378
379 if (!a)
380 return 0;
381 if ((a & 1) == 0)
382 return 2;
383 lim = IntSqrt (a);
384 for (v = 3; v <= lim; v += 2)
385 {
386 if (a % v == 0)
387 return v;
388 }
389 return a;
390 }
391
392 static int
393 IntPow (int a, int p)
394 {
395 int result;
396
397 result = 1;
398 while (p)
399 {
400 if (p & 1)
401 result = result * a;
402 p >>= 1;
403 if (p)
404 a = a * a;
405 }
406 return result;
407 }
408
409 static int
410 IntPowMod (int a, int p, int m)
411 {
412 int result;
413
414 if (m >= 32767)
415 {
416 #if DIGITBITS == 32
417 signed_digit la = a, lm = m, lr;
418 lr = 1;
419 while (p)
420 {
421 if (p & 1)
422 lr = (lr * la) % lm;
423 p >>= 1;
424 if (p)
425 la = (la * la) % lm;
426 }
427 result = (int) lr;
428 #else
429 ENTER ();
430 result = NaturalToInt (NaturalPowMod (NewNatural (a),
431 NewNatural (p),
432 NewNatural (m)));
433 EXIT ();
434 #endif
435 }
436 else
437 {
438 result = 1;
439 while (p)
440 {
441 if (p & 1)
442 result = (result * a) % m;
443 p >>= 1;
444 if (p)
445 a = (a * a) % m;
446 }
447 }
448 return result;
449 }
450
451 static int
452 IntPsi (int a)
453 {
454 int p;
455 int n;
456 int ret;
457
458 ret = 1;
459 while (a != 1)
460 {
461 p = IntFactor (a);
462 n = 0;
463 do
464 {
465 n++;
466 a /= p;
467 } while (a % p == 0);
468 ret = ret * (IntPow (p, n-1) * (p - 1));
469 }
470 return ret;
471 }
472 #endif
473
474 typedef struct _partial {
475 DataType *data;
476 struct _partial *down;
477 Natural *partial;
478 int power;
479 } Partial, *PartialPtr;
480
PartialMark(void * object)481 static void PartialMark (void *object)
482 {
483 PartialPtr p = object;
484
485 MemReference (p->partial);
486 MemReference (p->down);
487 }
488
489 DataType PartialType = { PartialMark, 0, "PartialType" };
490
491 static PartialPtr
NewPartial(Natural * partial)492 NewPartial (Natural *partial)
493 {
494 ENTER ();
495 PartialPtr p;
496
497 if (!partial)
498 RETURN (0);
499 p = ALLOCATE (&PartialType, sizeof (Partial));
500 p->down = 0;
501 p->partial = partial;
502 p->power = 0;
503 RETURN (p);
504 }
505
506 typedef struct _factor {
507 DataType *data;
508 struct _factor *next;
509 Natural *prime;
510 int power;
511 PartialPtr partials;
512 } Factor, *FactorPtr;
513
FactorMark(void * object)514 static void FactorMark (void *object)
515 {
516 FactorPtr f = object;
517
518 MemReference (f->prime);
519 MemReference (f->next);
520 MemReference (f->partials);
521 }
522
523 DataType FactorType = { FactorMark, 0, "FactorType" };
524
525 static FactorPtr
NewFactor(Natural * prime,int power,FactorPtr next)526 NewFactor (Natural *prime, int power, FactorPtr next)
527 {
528 ENTER ();
529 FactorPtr f;
530
531 f = ALLOCATE (&FactorType, sizeof (Factor));
532 f->next = next;
533 f->prime = prime;
534 f->power = power;
535 f->partials = 0;
536 f->partials = NewPartial (prime);
537 f->partials->power = 0;
538 RETURN (f);
539 }
540
541 static FactorPtr
GenerateFactors(Natural * n,Natural * max)542 GenerateFactors (Natural *n, Natural *max)
543 {
544 ENTER ();
545 FactorPtr f = 0;
546 Natural *p;
547 Natural *largest;
548 Natural *d, *rem;
549
550 p = 0;
551 largest = NaturalSqrt (n);
552 while (!NaturalEqual (n, one_natural))
553 {
554 int power = 1;
555 for (;;)
556 {
557 if (!p)
558 p = two_natural;
559 else if (NaturalEqual (p, two_natural))
560 p = NewNatural (3);
561 else
562 p = NaturalPlus (p, two_natural);
563
564 d = NaturalDivide (n, p, &rem);
565 if (NaturalZero (rem))
566 break;
567 if (max && NaturalLess (max, p))
568 RETURN(f);
569 if (NaturalLess (largest, p))
570 RETURN (NewFactor (n, 1, f));
571 }
572 n = d;
573 for (;;)
574 {
575 d = NaturalDivide (n, p, &rem);
576 if (!NaturalZero (rem))
577 break;
578 n = d;
579 power++;
580 }
581 f = NewFactor (p, power, f);
582 largest = NaturalSqrt (n);
583 }
584 RETURN (f);
585 }
586
587 static Natural *
FactorBump(FactorPtr f)588 FactorBump (FactorPtr f)
589 {
590 PartialPtr p, minp;
591 Natural *factor;
592
593 ENTER ();
594 if (!f)
595 RETURN(0);
596 p = f->partials;
597 if (!p)
598 RETURN(0);
599 minp = p;
600 while (p->power)
601 {
602 if (!p->down)
603 p->down = NewPartial (FactorBump (f->next));
604 p = p->down;
605 if (!p)
606 break;
607 if (NaturalLess (p->partial, minp->partial))
608 minp = p;
609 }
610 if (!minp)
611 RETURN(0);
612 factor = minp->partial;
613 if (minp->power < f->power)
614 {
615 minp->partial = NaturalTimes (minp->partial, f->prime);
616 minp->power++;
617 }
618 else
619 {
620 f->partials = minp->down;
621 }
622 RETURN (factor);
623 }
624
625 static int
RationalRepeatLength(int prec,Natural * nden,int ibase)626 RationalRepeatLength (int prec, Natural *nden, int ibase)
627 {
628 ENTER ();
629 Natural *nbase;
630 Natural *ndigits;
631 FactorPtr factors;
632 Natural *factor;
633 int digits;
634 Natural *max = 0;
635
636 if (NaturalEqual (nden, one_natural))
637 return 0;
638 if (prec > 0)
639 max = NewNatural (prec);
640 nbase = NewNatural (ibase);
641 ndigits = NaturalPsi (nden, max);
642 if (!ndigits)
643 {
644 factor = one_natural;
645 for (factor = one_natural;;
646 factor = NaturalPlus (factor, one_natural))
647 {
648 if (NaturalEqual (NaturalPowMod (nbase, factor, nden),
649 one_natural))
650 break;
651 if (aborting)
652 break;
653 if (NaturalLess (max, factor))
654 {
655 EXIT ();
656 return -1;
657 }
658 }
659 }
660 else
661 {
662 factors = GenerateFactors (ndigits, max);
663 if (aborting)
664 return 0;
665 factor = one_natural;
666 while (factor)
667 {
668 if (NaturalEqual (NaturalPowMod (nbase, factor, nden),
669 one_natural))
670 break;
671 if (aborting)
672 break;
673 factor = FactorBump (factors);
674 if (max && factor && NaturalLess (max, factor))
675 {
676 EXIT ();
677 return -1;
678 }
679 }
680 }
681 if (!factor)
682 factor = ndigits;
683 if (NaturalLess (max_int_natural, factor))
684 factor = max_int_natural;
685 digits = NaturalToInt (factor);
686 EXIT ();
687 return digits;
688 }
689
690 static void
CheckDecimalLength(int prec,Natural * nden,int ibase,int * initial,int * repeat)691 CheckDecimalLength (int prec, Natural *nden, int ibase, int *initial, int *repeat)
692 {
693 ENTER ();
694 Natural *rem;
695 Natural *nbase;
696 Natural *g;
697 int offset;
698 int rep;
699
700 nbase = NewNatural (ibase);
701 offset = 0;
702 while (!NaturalEqual ((g = NaturalGcd (nden, nbase)), one_natural))
703 {
704 if (aborting)
705 {
706 EXIT ();
707 return;
708 }
709 offset++;
710 if (prec >= 0 && offset > prec)
711 break;
712 nden = NaturalDivide (nden, g, &rem);
713 }
714 if (prec >= 0 && offset >= prec)
715 {
716 if (offset > prec)
717 offset = -prec;
718 else
719 offset = prec;
720 rep = 0;
721 }
722 else if (NaturalEqual (nden, one_natural))
723 {
724 rep = 0;
725 }
726 else
727 {
728 if (prec >= 0)
729 prec -= offset;
730 rep = RationalRepeatLength (prec, nden, ibase);
731 }
732 *initial = offset;
733 *repeat = rep;
734 EXIT ();
735 }
736
737 static Bool
RationalDecimalPrint(Value f,Value rv,char format,int base,int width,int prec,int fill)738 RationalDecimalPrint (Value f, Value rv, char format, int base, int width, int prec, int fill)
739 {
740 ENTER ();
741 Rational *r = &rv->rational;
742 Natural *quo;
743 Natural *partial;
744 Natural *rep, *init;
745 Natural *dig;
746 int exponent = 0;
747 int exponent_width = 0;
748 char *initial = 0, *in;
749 char *repeat = 0, *re;
750 char *whole;
751 int initial_width, repeat_width = 0;
752 int frac_width;
753 int rep_width, brace_width = 0, dot_width = 0;
754 int whole_width;
755 int fraction_width;
756 int print_width;
757 int min_prec;
758 Bool use_braces = True;
759
760 min_prec = 0;
761 if (format == 'f' || format == 'e' || format == 'g')
762 {
763 min_prec = prec;
764 use_braces = False;
765 }
766 if (prec == DEFAULT_OUTPUT_PRECISION)
767 {
768 min_prec = 0;
769 prec = 15;
770 }
771 else if (prec == INFINITE_OUTPUT_PRECISION)
772 prec = -1;
773 dig = NewNatural (base);
774 /*
775 * Check for small numbers for 'e' format
776 */
777 if (NaturalLess (r->num, r->den) && !NaturalZero(r->num))
778 {
779 Natural *quo, *rem;
780 Natural *mag;
781 int bits;
782
783 if (format == 'e' || (format == 'g' && prec > 0))
784 {
785 quo = NaturalDivide (r->den, r->num, &rem);
786 bits = NaturalWidth (quo);
787 exponent = (int) ((double) bits / (log ((double) base) / log (2.0)));
788 if (exponent < 0)
789 exponent = 0;
790 mag = NaturalIntPow (dig, exponent);
791 while (NaturalLess (mag, quo))
792 {
793 mag = NaturalTimes (mag, dig);
794 exponent++;
795 }
796 if (format == 'g' && prec > 0)
797 if (prec - exponent < 3)
798 format = 'e';
799 if (format == 'e')
800 {
801 int ev;
802 rv = RationalTimes (rv, NewRational (Positive, mag, one_natural), True);
803 r = &rv->rational;
804 exponent_width = 3;
805 ev = exponent;
806 while (ev >= base)
807 {
808 exponent_width++;
809 ev /= base;
810 }
811 exponent = -exponent;
812 }
813 else
814 exponent = 0;
815 }
816 else
817 exponent = 0;
818 }
819 CheckDecimalLength (prec, r->den, base, &initial_width, &repeat_width);
820 if (aborting)
821 {
822 EXIT ();
823 return False;
824 }
825 if ((rep_width = repeat_width))
826 {
827 /*
828 * When using %f format, just fill the
829 * result with digits
830 */
831 if (!use_braces && prec != -1)
832 {
833 initial_width = -prec;
834 repeat_width = 0;
835 rep_width = 0;
836 }
837 else
838 {
839 if (repeat_width < 0)
840 rep_width = prec - initial_width;
841 }
842 }
843 if (initial_width)
844 {
845 Natural *half_digit;
846 if (initial_width < 0)
847 {
848 initial_width = -initial_width;
849 half_digit = NaturalTimes (NaturalIntPow (dig, initial_width),
850 two_natural);
851 rv = RationalPlusHelper (r->sign,
852 r,
853 &NewRational (Positive,
854 one_natural,
855 half_digit)->rational);
856 r = &rv->rational;
857 }
858 else
859 {
860 if (!repeat_width && initial_width < min_prec)
861 initial_width = min_prec;
862 }
863 initial = malloc (initial_width + 1);
864 if (!initial)
865 {
866 EXIT ();
867 return False;
868 }
869 }
870 quo = NaturalDivide (r->num, r->den, &partial);
871 whole = NaturalSprint (0, quo, base, &whole_width);
872 brace_width = 0;
873 if (repeat_width)
874 {
875 brace_width++;
876 if (repeat_width > 0)
877 brace_width++;
878 }
879 dot_width = 0;
880 if (initial_width + rep_width)
881 dot_width = 1;
882 /*
883 * Compute how much space is available for the fractional part
884 */
885 if (width)
886 {
887 if (width < 0)
888 fraction_width = -width;
889 else
890 fraction_width = width;
891 fraction_width = fraction_width - (whole_width + exponent_width);
892 if (fraction_width < 0)
893 fraction_width = 0;
894 if (prec > 0 && fraction_width > prec + dot_width)
895 fraction_width = prec + dot_width;
896 }
897 else if (prec > 0)
898 fraction_width = prec + dot_width;
899 else
900 fraction_width = -1;
901 /*
902 * Start paring down parts of the output to fit the desired size
903 */
904 while (fraction_width >= 0 &&
905 (frac_width = dot_width + initial_width + rep_width + brace_width)
906 && frac_width > fraction_width)
907 {
908 if (rep_width)
909 {
910 if (brace_width > 1)
911 {
912 brace_width = 1;
913 repeat_width = -repeat_width;
914 }
915 rep_width = fraction_width - (dot_width + initial_width +
916 brace_width);
917 if (rep_width < 0)
918 {
919 rep_width = 0;
920 }
921 }
922 else if (brace_width)
923 brace_width = 0;
924 else if (initial_width)
925 {
926 initial_width = fraction_width - dot_width;
927 if (initial_width < 0)
928 initial_width = 0;
929 }
930 else
931 dot_width = 0;
932 }
933 if (initial_width)
934 {
935 init = NaturalDivide (NaturalTimes (partial,
936 NaturalIntPow (dig, initial_width)),
937 r->den,
938 &partial);
939 if (aborting)
940 {
941 free (initial);
942 EXIT ();
943 return False;
944 }
945 in = NaturalSprint (initial + initial_width + 1,
946 init, base, &initial_width);
947 if (!in)
948 {
949 free (initial);
950 EXIT ();
951 return False;
952 }
953 while (in > initial)
954 {
955 *--in = '0';
956 ++initial_width;
957 }
958 }
959 if (rep_width)
960 {
961 #define MAX_SENSIBLE 10000000
962 if (rep_width > MAX_SENSIBLE)
963 {
964 repeat_width = -1;
965 rep_width = MAX_SENSIBLE;
966 }
967 /*
968 * allocate the output buffer; keep trying until this works
969 */
970 while (!(repeat = malloc (rep_width + 1)))
971 {
972 repeat_width = -1;
973 rep_width >>= 1;
974 }
975 rep = NaturalDivide (NaturalTimes (partial,
976 NaturalIntPow (dig, rep_width)),
977 r->den,
978 &partial);
979 if (aborting)
980 {
981 free (initial);
982 free (repeat);
983 EXIT ();
984 return False;
985 }
986 re = NaturalSprint (repeat + rep_width + 1,
987 rep, base, &rep_width);
988 if (!re)
989 {
990 free (initial);
991 free (repeat);
992 EXIT ();
993 return False;
994 }
995 while (re > repeat)
996 {
997 *--re = '0';
998 ++rep_width;
999 }
1000 if (use_braces)
1001 {
1002 rep_width++; /* open { */
1003 if (repeat_width > 0)
1004 rep_width++; /* close } */
1005 }
1006 }
1007 fraction_width = initial_width + rep_width;
1008 print_width = whole_width + 1 + fraction_width + exponent_width;
1009 if (r->sign == Negative)
1010 print_width = print_width + 1;
1011 while (width > print_width)
1012 {
1013 FileOutchar (f, fill);
1014 width--;
1015 }
1016 if (r->sign == Negative)
1017 FileOutput (f, '-');
1018 FilePuts (f, whole);
1019 FileOutput (f, '.');
1020 if (initial_width)
1021 {
1022 FilePuts (f, initial);
1023 free (initial);
1024 }
1025 if (rep_width)
1026 {
1027 if (use_braces)
1028 FileOutput (f, '{');
1029 FilePuts (f, repeat);
1030 free (repeat);
1031 if (use_braces && repeat_width > 0)
1032 FileOutput (f, '}');
1033 }
1034 if (exponent)
1035 {
1036 FilePrintf (f, "e%d", exponent);
1037 }
1038 while (-width > print_width)
1039 {
1040 FileOutchar (f, fill);
1041 width++;
1042 }
1043 EXIT ();
1044 return True;
1045 }
1046
1047 static Bool
RationalPrint(Value f,Value rv,char format,int base,int width,int prec,int fill)1048 RationalPrint (Value f, Value rv, char format, int base, int width, int prec, int fill)
1049 {
1050 Rational *r = &rv->rational;
1051 char *num, *num_base, *den, *den_base;
1052 int num_width, den_width;
1053 int print_width;
1054 Bool ret = True;
1055
1056 if (base == 0)
1057 base = 10;
1058 switch (format) {
1059 case 'v':
1060 num_width = NaturalEstimateLength (r->num, base);
1061 num_base = malloc (num_width);
1062 num = NaturalSprint (num_base + num_width, r->num, base, &num_width);
1063 if (!num)
1064 {
1065 free (num_base);
1066 ret = False;
1067 break;
1068 }
1069 den_width = NaturalEstimateLength (r->den, base);
1070 den_base = malloc (den_width);
1071 den = NaturalSprint (den_base + den_width, r->den, base, &den_width);
1072 if (!den)
1073 {
1074 free (num_base);
1075 free (den_base);
1076 ret = False;
1077 break;
1078 }
1079 print_width = 1 + num_width + 1 + den_width + 1;
1080 if (r->sign == Negative)
1081 print_width++;
1082 while (width > print_width)
1083 {
1084 FileOutchar (f, fill);
1085 width--;
1086 }
1087 FileOutput (f, '(');
1088 if (r->sign == Negative)
1089 FileOutput (f, '-');
1090 FilePuts (f, num);
1091 FileOutput (f, '/');
1092 FilePuts (f, den);
1093 FileOutput (f, ')');
1094 free (num_base);
1095 free (den_base);
1096 while (-width > print_width)
1097 {
1098 FileOutchar (f, fill);
1099 width++;
1100 }
1101 break;
1102 default:
1103 ret = RationalDecimalPrint (f, rv, format, base, width, prec, fill);
1104 break;
1105 }
1106 return ret;
1107 }
1108
1109 static void
RationalMark(void * object)1110 RationalMark (void *object)
1111 {
1112 Rational *rational = object;
1113
1114 MemReference (rational->num);
1115 MemReference (rational->den);
1116 }
1117
1118 ValueRep RationalRep = {
1119 { RationalMark, 0, "RationalRep" }, /* base */
1120 rep_rational, /* tag */
1121 { /* binary */
1122 RationalPlus,
1123 RationalMinus,
1124 RationalTimes,
1125 RationalDivide,
1126 NumericDiv,
1127 RationalMod,
1128 RationalLess,
1129 RationalEqual,
1130 0,
1131 0,
1132 },
1133 { /* unary */
1134 RationalNegate,
1135 RationalFloor,
1136 RationalCeil,
1137 },
1138 RationalPromote,
1139 RationalReduce,
1140 RationalPrint,
1141 0,
1142 RationalHash,
1143 };
1144
1145 Value
NewRational(Sign sign,Natural * num,Natural * den)1146 NewRational (Sign sign, Natural *num, Natural *den)
1147 {
1148 ENTER ();
1149 Value ret;
1150 Natural *g;
1151 Natural *rem;
1152
1153 if (NaturalZero (num))
1154 den = one_natural;
1155 else
1156 {
1157 if (NaturalLength(den) != 1 || NaturalDigits(den)[0] != 1)
1158 {
1159 g = NaturalGcd (num, den);
1160 if (NaturalLength (g) != 1 || NaturalDigits(g)[0] != 1)
1161 {
1162 num = NaturalDivide (num, g, &rem);
1163 den = NaturalDivide (den, g, &rem);
1164 }
1165 }
1166 }
1167 ret = ALLOCATE (&RationalRep.data, sizeof (Rational));
1168 ret->rational.sign = sign;
1169 ret->rational.num = num;
1170 ret->rational.den = den;
1171 RETURN (ret);
1172 }
1173
1174 Value
NewIntRational(int i)1175 NewIntRational (int i)
1176 {
1177 ENTER ();
1178 if (i < 0)
1179 RETURN (NewRational (Negative, NewNatural ((unsigned) -i), one_natural));
1180 else
1181 RETURN (NewRational (Positive, NewNatural ((unsigned) i), one_natural));
1182 }
1183
1184 Value
NewIntegerRational(Integer * i)1185 NewIntegerRational (Integer *i)
1186 {
1187 ENTER ();
1188 RETURN (NewRational (IntegerSign((Value) i), IntegerMag((Value) i),
1189 one_natural));
1190 }
1191