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