add a comment.
[ruby-svn.git] / rational.c
blob48ed7eae1772e1966e664a98284a906e56a09a6e
1 /*
2 rational.c: Coded by Tadayoshi Funaba 2008
4 This implementation is based on Keiju Ishitsuka's Rational library
5 which is written in ruby.
6 */
8 #include "ruby.h"
9 #include <math.h>
10 #include <float.h>
12 #define NDEBUG
13 #include <assert.h>
15 #ifdef HAVE_IEEEFP_H
16 #include <ieeefp.h>
17 #endif
19 #ifndef RATIONAL_NAME
20 #define RATIONAL_NAME "Rational"
21 #endif
23 #define ZERO INT2FIX(0)
24 #define ONE INT2FIX(1)
25 #define TWO INT2FIX(2)
27 VALUE rb_cRational;
29 static ID id_Unify, id_abs, id_cmp, id_convert, id_equal_p,
30 id_expt, id_floor, id_format, id_idiv, id_inspect, id_negate, id_new,
31 id_new_bang, id_to_f, id_to_i, id_to_s, id_truncate;
33 #define f_boolcast(x) ((x) ? Qtrue : Qfalse)
35 #define binop(n,op) \
36 inline static VALUE \
37 f_##n(VALUE x, VALUE y)\
39 return rb_funcall(x, op, 1, y);\
42 #define fun1(n) \
43 inline static VALUE \
44 f_##n(VALUE x)\
46 return rb_funcall(x, id_##n, 0);\
49 #define fun2(n) \
50 inline static VALUE \
51 f_##n(VALUE x, VALUE y)\
53 return rb_funcall(x, id_##n, 1, y);\
56 inline static VALUE
57 f_add(VALUE x, VALUE y)
59 if (FIXNUM_P(y)) {
60 if (FIX2LONG(y) == 0)
61 return x;
63 else if (FIXNUM_P(x)) {
64 if (FIX2LONG(x) == 0)
65 return y;
67 return rb_funcall(x, '+', 1, y);
70 inline static VALUE
71 f_cmp(VALUE x, VALUE y)
73 if (FIXNUM_P(x) && FIXNUM_P(y)) {
74 long c = FIX2LONG(x) - FIX2LONG(y);
75 if (c > 0)
76 c = 1;
77 else if (c < 0)
78 c = -1;
79 return INT2FIX(c);
81 return rb_funcall(x, id_cmp, 1, y);
84 inline static VALUE
85 f_div(VALUE x, VALUE y)
87 if (FIXNUM_P(y) && FIX2LONG(y) == 1)
88 return x;
89 return rb_funcall(x, '/', 1, y);
92 inline static VALUE
93 f_gt_p(VALUE x, VALUE y)
95 if (FIXNUM_P(x) && FIXNUM_P(y))
96 return f_boolcast(FIX2LONG(x) > FIX2LONG(y));
97 return rb_funcall(x, '>', 1, y);
100 inline static VALUE
101 f_lt_p(VALUE x, VALUE y)
103 if (FIXNUM_P(x) && FIXNUM_P(y))
104 return f_boolcast(FIX2LONG(x) < FIX2LONG(y));
105 return rb_funcall(x, '<', 1, y);
108 binop(mod, '%')
110 inline static VALUE
111 f_mul(VALUE x, VALUE y)
113 if (FIXNUM_P(y)) {
114 long _iy = FIX2LONG(y);
115 if (_iy == 0) {
116 if (TYPE(x) == T_FLOAT)
117 return rb_float_new(0.0);
118 else
119 return ZERO;
121 else if (_iy == 1)
122 return x;
124 else if (FIXNUM_P(x)) {
125 long _ix = FIX2LONG(x);
126 if (_ix == 0) {
127 if (TYPE(y) == T_FLOAT)
128 return rb_float_new(0.0);
129 else
130 return ZERO;
132 else if (_ix == 1)
133 return y;
135 return rb_funcall(x, '*', 1, y);
138 inline static VALUE
139 f_sub(VALUE x, VALUE y)
141 if (FIXNUM_P(y))
142 if (FIX2LONG(y) == 0)
143 return x;
144 return rb_funcall(x, '-', 1, y);
147 binop(xor, '^')
149 fun1(abs)
150 fun1(floor)
151 fun1(inspect)
152 fun1(negate)
153 fun1(to_f)
154 fun1(to_i)
155 fun1(to_s)
156 fun1(truncate)
158 inline static VALUE
159 f_equal_p(VALUE x, VALUE y)
161 if (FIXNUM_P(x) && FIXNUM_P(y))
162 return f_boolcast(FIX2LONG(x) == FIX2LONG(y));
163 return rb_funcall(x, id_equal_p, 1, y);
166 fun2(expt)
167 fun2(idiv)
169 inline static VALUE
170 f_negative_p(VALUE x)
172 if (FIXNUM_P(x))
173 return f_boolcast(FIX2LONG(x) < 0);
174 return rb_funcall(x, '<', 1, ZERO);
177 inline static VALUE
178 f_zero_p(VALUE x)
180 if (FIXNUM_P(x))
181 return f_boolcast(FIX2LONG(x) == 0);
182 return rb_funcall(x, id_equal_p, 1, ZERO);
185 inline static VALUE
186 f_one_p(VALUE x)
188 if (FIXNUM_P(x))
189 return f_boolcast(FIX2LONG(x) == 1);
190 return rb_funcall(x, id_equal_p, 1, ONE);
193 inline static VALUE
194 f_kind_of_p(VALUE x, VALUE c)
196 return rb_obj_is_kind_of(x, c);
199 inline static VALUE
200 k_numeric_p(VALUE x)
202 return f_kind_of_p(x, rb_cNumeric);
205 inline static VALUE
206 k_integer_p(VALUE x)
208 return f_kind_of_p(x, rb_cInteger);
211 inline static VALUE
212 k_float_p(VALUE x)
214 return f_kind_of_p(x, rb_cFloat);
217 inline static VALUE
218 k_rational_p(VALUE x)
220 return f_kind_of_p(x, rb_cRational);
223 #ifndef NDEBUG
224 #define f_gcd f_gcd_orig
225 #endif
227 inline static long
228 i_gcd(long x, long y)
230 long b;
232 if (x < 0)
233 x = -x;
234 if (y < 0)
235 y = -y;
237 if (x == 0)
238 return y;
239 if (y == 0)
240 return x;
242 b = 0;
243 while ((x & 1) == 0 && (y & 1) == 0) {
244 b += 1;
245 x >>= 1;
246 y >>= 1;
249 while ((x & 1) == 0)
250 x >>= 1;
252 while ((y & 1) == 0)
253 y >>= 1;
255 while (x != y) {
256 if (y > x) {
257 long t;
258 t = x;
259 x = y;
260 y = t;
262 x -= y;
263 while ((x & 1) == 0)
264 x >>= 1;
267 return x << b;
270 inline static VALUE
271 f_gcd(VALUE x, VALUE y)
273 VALUE z;
275 if (FIXNUM_P(x) && FIXNUM_P(y))
276 return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
278 if (f_negative_p(x))
279 x = f_negate(x);
280 if (f_negative_p(y))
281 y = f_negate(y);
283 if (f_zero_p(x))
284 return y;
285 if (f_zero_p(y))
286 return x;
288 for (;;) {
289 if (FIXNUM_P(x)) {
290 if (FIX2LONG(x) == 0)
291 return y;
292 if (FIXNUM_P(y))
293 return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
295 z = x;
296 x = f_mod(y, x);
297 y = z;
299 /* NOTREACHED */
302 #ifndef NDEBUG
303 #undef f_gcd
305 inline static VALUE
306 f_gcd(VALUE x, VALUE y)
308 VALUE r = f_gcd_orig(x, y);
309 if (!f_zero_p(r)) {
310 assert(f_zero_p(f_mod(x, r)));
311 assert(f_zero_p(f_mod(y, r)));
313 return r;
315 #endif
317 inline static VALUE
318 f_lcm(VALUE x, VALUE y)
320 if (f_zero_p(x) || f_zero_p(y))
321 return ZERO;
322 return f_abs(f_mul(f_div(x, f_gcd(x, y)), y));
325 #define get_dat1(x) \
326 struct RRational *dat;\
327 dat = ((struct RRational *)(x))
329 #define get_dat2(x,y) \
330 struct RRational *adat, *bdat;\
331 adat = ((struct RRational *)(x));\
332 bdat = ((struct RRational *)(y))
334 inline static VALUE
335 nurat_s_new_internal(VALUE klass, VALUE num, VALUE den)
337 NEWOBJ(obj, struct RRational);
338 OBJSETUP(obj, klass, T_RATIONAL);
340 obj->num = num;
341 obj->den = den;
343 return (VALUE)obj;
346 static VALUE
347 nurat_s_alloc(VALUE klass)
349 return nurat_s_new_internal(klass, ZERO, ONE);
352 static VALUE
353 nurat_s_new_bang(int argc, VALUE *argv, VALUE klass)
355 VALUE num, den;
357 switch (rb_scan_args(argc, argv, "11", &num, &den)) {
358 case 1:
359 if (!k_integer_p(num))
360 num = f_to_i(num);
361 den = ONE;
362 break;
363 default:
364 if (!k_integer_p(num))
365 num = f_to_i(num);
366 if (!k_integer_p(den))
367 den = f_to_i(den);
369 switch (FIX2INT(f_cmp(den, ZERO))) {
370 case -1:
371 num = f_negate(num);
372 den = f_negate(den);
373 break;
374 case 0:
375 rb_raise(rb_eZeroDivError, "devided by zero");
376 break;
378 break;
381 return nurat_s_new_internal(klass, num, den);
384 inline static VALUE
385 f_rational_new_bang1(VALUE klass, VALUE x)
387 return nurat_s_new_internal(klass, x, ONE);
390 inline static VALUE
391 f_rational_new_bang2(VALUE klass, VALUE x, VALUE y)
393 assert(!f_negative_p(y));
394 assert(!f_zero_p(y));
395 return nurat_s_new_internal(klass, x, y);
398 #define f_unify_p(klass) rb_const_defined(klass, id_Unify)
400 inline static void
401 nurat_int_check(VALUE num)
403 switch (TYPE(num)) {
404 case T_FIXNUM:
405 case T_BIGNUM:
406 break;
407 default:
408 rb_raise(rb_eArgError, "not an integer");
412 inline static VALUE
413 nurat_s_canonicalize_internal(VALUE klass, VALUE num, VALUE den)
415 VALUE gcd;
417 switch (FIX2INT(f_cmp(den, ZERO))) {
418 case -1:
419 num = f_negate(num);
420 den = f_negate(den);
421 break;
422 case 0:
423 rb_raise(rb_eZeroDivError, "devided by zero");
424 break;
427 gcd = f_gcd(num, den);
428 num = f_idiv(num, gcd);
429 den = f_idiv(den, gcd);
431 if (f_one_p(den) && f_unify_p(klass))
432 return num;
433 return nurat_s_new_internal(klass, num, den);
436 inline static VALUE
437 nurat_s_canonicalize_internal_no_reduce(VALUE klass, VALUE num, VALUE den)
439 switch (FIX2INT(f_cmp(den, ZERO))) {
440 case -1:
441 num = f_negate(num);
442 den = f_negate(den);
443 break;
444 case 0:
445 rb_raise(rb_eZeroDivError, "devided by zero");
446 break;
449 if (f_equal_p(den, ONE) && f_unify_p(klass))
450 return num;
451 return nurat_s_new_internal(klass, num, den);
454 #if 0
455 static VALUE
456 nurat_s_canonicalize(int argc, VALUE *argv, VALUE klass)
458 VALUE num, den;
460 if (rb_scan_args(argc, argv, "11", &num, &den) == 1) {
461 den = ONE;
464 nurat_int_check(num);
465 nurat_int_check(den);
467 return nurat_s_canonicalize_internal(klass, num, den);
469 #endif
471 static VALUE
472 nurat_s_new(int argc, VALUE *argv, VALUE klass)
474 VALUE num, den;
476 if (rb_scan_args(argc, argv, "11", &num, &den) == 1)
477 den = ONE;
479 nurat_int_check(num);
480 nurat_int_check(den);
482 return nurat_s_canonicalize_internal(klass, num, den);
485 inline static VALUE
486 f_rational_new1(VALUE klass, VALUE x)
488 assert(!k_rational_p(x));
489 return nurat_s_canonicalize_internal(klass, x, ONE);
492 inline static VALUE
493 f_rational_new2(VALUE klass, VALUE x, VALUE y)
495 assert(!k_rational_p(x));
496 assert(!k_rational_p(y));
497 return nurat_s_canonicalize_internal(klass, x, y);
500 inline static VALUE
501 f_rational_new_no_reduce1(VALUE klass, VALUE x)
503 assert(!k_rational_p(x));
504 return nurat_s_canonicalize_internal_no_reduce(klass, x, ONE);
507 inline static VALUE
508 f_rational_new_no_reduce2(VALUE klass, VALUE x, VALUE y)
510 assert(!k_rational_p(x));
511 assert(!k_rational_p(y));
512 return nurat_s_canonicalize_internal_no_reduce(klass, x, y);
515 static VALUE
516 nurat_f_rational(int argc, VALUE *argv, VALUE klass)
518 return rb_funcall2(rb_cRational, id_convert, argc, argv);
521 static VALUE
522 nurat_numerator(VALUE self)
524 get_dat1(self);
525 return dat->num;
528 static VALUE
529 nurat_denominator(VALUE self)
531 get_dat1(self);
532 return dat->den;
535 #ifndef NDEBUG
536 #define f_imul f_imul_orig
537 #endif
539 inline static VALUE
540 f_imul(long a, long b)
542 VALUE r;
543 long c;
545 if (a == 0 || b == 0)
546 return ZERO;
547 else if (a == 1)
548 return LONG2NUM(b);
549 else if (b == 1)
550 return LONG2NUM(a);
552 c = a * b;
553 r = LONG2NUM(c);
554 if (NUM2LONG(r) != c || (c / a) != b)
555 r = rb_big_mul(rb_int2big(a), rb_int2big(b));
556 return r;
559 #ifndef NDEBUG
560 #undef f_imul
562 inline static VALUE
563 f_imul(long x, long y)
565 VALUE r = f_imul_orig(x, y);
566 assert(f_equal_p(r, f_mul(LONG2NUM(x), LONG2NUM(y))));
567 return r;
569 #endif
571 inline static VALUE
572 f_addsub(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
574 VALUE num, den;
576 if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
577 FIXNUM_P(bnum) && FIXNUM_P(bden)) {
578 long an = FIX2LONG(anum);
579 long ad = FIX2LONG(aden);
580 long bn = FIX2LONG(bnum);
581 long bd = FIX2LONG(bden);
582 long ig = i_gcd(ad, bd);
584 VALUE g = LONG2NUM(ig);
585 VALUE a = f_imul(an, bd / ig);
586 VALUE b = f_imul(bn, ad / ig);
587 VALUE c;
589 if (k == '+')
590 c = f_add(a, b);
591 else
592 c = f_sub(a, b);
594 b = f_idiv(aden, g);
595 g = f_gcd(c, g);
596 num = f_idiv(c, g);
597 a = f_idiv(bden, g);
598 den = f_mul(a, b);
600 else {
601 VALUE g = f_gcd(aden, bden);
602 VALUE a = f_mul(anum, f_idiv(bden, g));
603 VALUE b = f_mul(bnum, f_idiv(aden, g));
604 VALUE c;
606 if (k == '+')
607 c = f_add(a, b);
608 else
609 c = f_sub(a, b);
611 b = f_idiv(aden, g);
612 g = f_gcd(c, g);
613 num = f_idiv(c, g);
614 a = f_idiv(bden, g);
615 den = f_mul(a, b);
617 return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
620 static VALUE
621 nurat_add(VALUE self, VALUE other)
623 switch (TYPE(other)) {
624 case T_FIXNUM:
625 case T_BIGNUM:
627 get_dat1(self);
629 return f_addsub(self,
630 dat->num, dat->den,
631 other, ONE, '+');
633 case T_FLOAT:
634 return f_add(f_to_f(self), other);
635 case T_RATIONAL:
637 get_dat2(self, other);
639 return f_addsub(self,
640 adat->num, adat->den,
641 bdat->num, bdat->den, '+');
643 default:
644 return rb_num_coerce_bin(self, other, '+');
648 static VALUE
649 nurat_sub(VALUE self, VALUE other)
651 switch (TYPE(other)) {
652 case T_FIXNUM:
653 case T_BIGNUM:
655 get_dat1(self);
657 return f_addsub(self,
658 dat->num, dat->den,
659 other, ONE, '-');
661 case T_FLOAT:
662 return f_sub(f_to_f(self), other);
663 case T_RATIONAL:
665 get_dat2(self, other);
667 return f_addsub(self,
668 adat->num, adat->den,
669 bdat->num, bdat->den, '-');
671 default:
672 return rb_num_coerce_bin(self, other, '-');
676 inline static VALUE
677 f_muldiv(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
679 VALUE num, den;
681 if (k == '/') {
682 VALUE t;
684 if (f_negative_p(bnum)) {
685 anum = f_negate(anum);
686 bnum = f_negate(bnum);
688 t = bnum;
689 bnum = bden;
690 bden = t;
693 if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
694 FIXNUM_P(bnum) && FIXNUM_P(bden)) {
695 long an = FIX2LONG(anum);
696 long ad = FIX2LONG(aden);
697 long bn = FIX2LONG(bnum);
698 long bd = FIX2LONG(bden);
699 long g1 = i_gcd(an, bd);
700 long g2 = i_gcd(ad, bn);
702 num = f_imul(an / g1, bn / g2);
703 den = f_imul(ad / g2, bd / g1);
705 else {
706 VALUE g1 = f_gcd(anum, bden);
707 VALUE g2 = f_gcd(aden, bnum);
709 num = f_mul(f_idiv(anum, g1), f_idiv(bnum, g2));
710 den = f_mul(f_idiv(aden, g2), f_idiv(bden, g1));
712 return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
715 static VALUE
716 nurat_mul(VALUE self, VALUE other)
718 switch (TYPE(other)) {
719 case T_FIXNUM:
720 case T_BIGNUM:
722 get_dat1(self);
724 return f_muldiv(self,
725 dat->num, dat->den,
726 other, ONE, '*');
728 case T_FLOAT:
729 return f_mul(f_to_f(self), other);
730 case T_RATIONAL:
732 get_dat2(self, other);
734 return f_muldiv(self,
735 adat->num, adat->den,
736 bdat->num, bdat->den, '*');
738 default:
739 return rb_num_coerce_bin(self, other, '*');
743 static VALUE
744 nurat_div(VALUE self, VALUE other)
746 switch (TYPE(other)) {
747 case T_FIXNUM:
748 case T_BIGNUM:
749 if (f_zero_p(other))
750 rb_raise(rb_eZeroDivError, "devided by zero");
752 get_dat1(self);
754 return f_muldiv(self,
755 dat->num, dat->den,
756 other, ONE, '/');
758 case T_FLOAT:
759 return rb_funcall(f_to_f(self), '/', 1, other);
760 case T_RATIONAL:
761 if (f_zero_p(other))
762 rb_raise(rb_eZeroDivError, "devided by zero");
764 get_dat2(self, other);
766 return f_muldiv(self,
767 adat->num, adat->den,
768 bdat->num, bdat->den, '/');
770 default:
771 return rb_num_coerce_bin(self, other, '/');
775 static VALUE
776 nurat_fdiv(VALUE self, VALUE other)
778 return f_div(f_to_f(self), other);
781 static VALUE
782 nurat_expt(VALUE self, VALUE other)
784 if (f_zero_p(other))
785 return f_rational_new_bang1(CLASS_OF(self), ONE);
787 if (k_rational_p(other)) {
788 get_dat1(other);
790 if (f_one_p(dat->den))
791 other = dat->num; /* good? */
794 switch (TYPE(other)) {
795 case T_FIXNUM:
796 case T_BIGNUM:
798 VALUE num, den;
800 get_dat1(self);
802 switch (FIX2INT(f_cmp(other, ZERO))) {
803 case 1:
804 num = f_expt(dat->num, other);
805 den = f_expt(dat->den, other);
806 break;
807 case -1:
808 num = f_expt(dat->den, f_negate(other));
809 den = f_expt(dat->num, f_negate(other));
810 break;
811 default:
812 num = ONE;
813 den = ONE;
814 break;
816 return f_rational_new2(CLASS_OF(self), num, den);
818 case T_FLOAT:
819 case T_RATIONAL:
820 return f_expt(f_to_f(self), other);
821 default:
822 return rb_num_coerce_bin(self, other, id_expt);
826 static VALUE
827 nurat_cmp(VALUE self, VALUE other)
829 switch (TYPE(other)) {
830 case T_FIXNUM:
831 case T_BIGNUM:
833 get_dat1(self);
835 if (FIXNUM_P(dat->den) && FIX2LONG(dat->den) == 1)
836 return f_cmp(dat->num, other);
837 return f_cmp(self, f_rational_new_bang1(CLASS_OF(self), other));
839 case T_FLOAT:
840 return f_cmp(f_to_f(self), other);
841 case T_RATIONAL:
843 VALUE num1, num2;
845 get_dat2(self, other);
847 if (FIXNUM_P(adat->num) && FIXNUM_P(adat->den) &&
848 FIXNUM_P(bdat->num) && FIXNUM_P(bdat->den)) {
849 num1 = f_imul(FIX2LONG(adat->num), FIX2LONG(bdat->den));
850 num2 = f_imul(FIX2LONG(bdat->num), FIX2LONG(adat->den));
852 else {
853 num1 = f_mul(adat->num, bdat->den);
854 num2 = f_mul(bdat->num, adat->den);
856 return f_cmp(f_sub(num1, num2), ZERO);
858 default:
859 return rb_num_coerce_bin(self, other, id_cmp);
863 static VALUE
864 nurat_equal_p(VALUE self, VALUE other)
866 switch (TYPE(other)) {
867 case T_FIXNUM:
868 case T_BIGNUM:
870 get_dat1(self);
872 if (f_zero_p(dat->num) && f_zero_p(other))
873 return Qtrue;
875 if (!FIXNUM_P(dat->den))
876 return Qfalse;
877 if (FIX2LONG(dat->den) != 1)
878 return Qfalse;
879 if (f_equal_p(dat->num, other))
880 return Qtrue;
881 return Qfalse;
883 case T_FLOAT:
884 return f_equal_p(f_to_f(self), other);
885 case T_RATIONAL:
887 get_dat2(self, other);
889 if (f_zero_p(adat->num) && f_zero_p(bdat->num))
890 return Qtrue;
892 return f_boolcast(f_equal_p(adat->num, bdat->num) &&
893 f_equal_p(adat->den, bdat->den));
895 default:
896 return f_equal_p(other, self);
900 static VALUE
901 nurat_coerce(VALUE self, VALUE other)
903 switch (TYPE(other)) {
904 case T_FIXNUM:
905 case T_BIGNUM:
906 return rb_assoc_new(f_rational_new_bang1(CLASS_OF(self), other), self);
907 case T_FLOAT:
908 return rb_assoc_new(other, f_to_f(self));
911 rb_raise(rb_eTypeError, "%s can't be coerced into %s",
912 rb_obj_classname(other), rb_obj_classname(self));
913 return Qnil;
916 static VALUE
917 nurat_idiv(VALUE self, VALUE other)
919 return f_floor(f_div(self, other));
922 static VALUE
923 nurat_mod(VALUE self, VALUE other)
925 VALUE val = f_floor(f_div(self, other));
926 return f_sub(self, f_mul(other, val));
929 static VALUE
930 nurat_divmod(VALUE self, VALUE other)
932 VALUE val = f_floor(f_div(self, other));
933 return rb_assoc_new(val, f_sub(self, f_mul(other, val)));
936 #if 0
937 static VALUE
938 nurat_quot(VALUE self, VALUE other)
940 return f_truncate(f_div(self, other));
942 #endif
944 static VALUE
945 nurat_rem(VALUE self, VALUE other)
947 VALUE val = f_truncate(f_div(self, other));
948 return f_sub(self, f_mul(other, val));
951 #if 0
952 static VALUE
953 nurat_quotrem(VALUE self, VALUE other)
955 VALUE val = f_truncate(f_div(self, other));
956 return rb_assoc_new(val, f_sub(self, f_mul(other, val)));
958 #endif
960 static VALUE
961 nurat_abs(VALUE self)
963 if (!f_negative_p(self))
964 return self;
965 return f_negate(self);
968 #if 0
969 static VALUE
970 nurat_true(VALUE self)
972 return Qtrue;
974 #endif
976 static VALUE
977 nurat_floor(VALUE self)
979 get_dat1(self);
980 return f_idiv(dat->num, dat->den);
983 static VALUE
984 nurat_ceil(VALUE self)
986 get_dat1(self);
987 return f_negate(f_idiv(f_negate(dat->num), dat->den));
990 static VALUE
991 nurat_truncate(VALUE self)
993 get_dat1(self);
994 if (f_negative_p(dat->num))
995 return f_negate(f_idiv(f_negate(dat->num), dat->den));
996 return f_idiv(dat->num, dat->den);
999 static VALUE
1000 nurat_round(VALUE self)
1002 get_dat1(self);
1004 if (f_negative_p(dat->num)) {
1005 VALUE num, den;
1007 num = f_negate(dat->num);
1008 num = f_add(f_mul(num, TWO), dat->den);
1009 den = f_mul(dat->den, TWO);
1010 return f_negate(f_idiv(num, den));
1012 else {
1013 VALUE num = f_add(f_mul(dat->num, TWO), dat->den);
1014 VALUE den = f_mul(dat->den, TWO);
1015 return f_idiv(num, den);
1019 #define f_size(x) rb_funcall(x, rb_intern("size"), 0)
1020 #define f_rshift(x,y) rb_funcall(x, rb_intern(">>"), 1, y)
1022 inline static long
1023 i_ilog2(VALUE x)
1025 long q, r, fx;
1027 assert(!f_lt_p(x, ONE));
1029 q = (NUM2LONG(f_size(x)) - sizeof(long)) * 8 + 1;
1031 if (q > 0)
1032 x = f_rshift(x, LONG2NUM(q));
1034 fx = NUM2LONG(x);
1036 r = -1;
1037 while (fx) {
1038 fx >>= 1;
1039 r += 1;
1042 return q + r;
1045 static VALUE
1046 nurat_to_f(VALUE self)
1048 VALUE num, den;
1049 int minus = 0;
1050 long nl, dl, ml, ne, de;
1051 int e;
1052 double f;
1055 get_dat1(self);
1057 if (f_zero_p(dat->num))
1058 return rb_float_new(0.0);
1060 num = dat->num;
1061 den = dat->den;
1064 if (f_negative_p(num)) {
1065 num = f_negate(num);
1066 minus = 1;
1069 nl = i_ilog2(num);
1070 dl = i_ilog2(den);
1071 ml = (long)(log(DBL_MAX) / log(2.0) - 1); /* should be a static */
1073 ne = 0;
1074 if (nl > ml) {
1075 ne = nl - ml;
1076 num = f_rshift(num, LONG2NUM(ne));
1079 de = 0;
1080 if (dl > ml) {
1081 de = dl - ml;
1082 den = f_rshift(den, LONG2NUM(de));
1085 e = (int)(ne - de);
1087 if ((e > DBL_MAX_EXP) || (e < DBL_MIN_EXP)) {
1088 rb_warning("%s out of Float range", rb_obj_classname(self));
1089 return rb_float_new(e > 0 ? HUGE_VAL : 0.0);
1092 f = NUM2DBL(num) / NUM2DBL(den);
1093 if (minus)
1094 f = -f;
1095 f = ldexp(f, e);
1097 if (isinf(f) || isnan(f))
1098 rb_warning("%s out of Float range", rb_obj_classname(self));
1100 return rb_float_new(f);
1103 static VALUE
1104 nurat_to_r(VALUE self)
1106 return self;
1109 static VALUE
1110 nurat_hash(VALUE self)
1112 get_dat1(self);
1113 return f_xor(dat->num, dat->den);
1116 static VALUE
1117 nurat_to_s(VALUE self)
1119 get_dat1(self);
1121 if (f_one_p(dat->den))
1122 return f_to_s(dat->num);
1123 return rb_funcall(rb_mKernel, id_format, 3,
1124 rb_str_new2("%d/%d"), dat->num, dat->den);
1127 static VALUE
1128 nurat_inspect(VALUE self)
1130 get_dat1(self);
1131 return rb_funcall(rb_mKernel, id_format, 3,
1132 rb_str_new2("Rational(%d, %d)"), dat->num, dat->den);
1135 static VALUE
1136 nurat_marshal_dump(VALUE self)
1138 get_dat1(self);
1139 return rb_assoc_new(dat->num, dat->den);
1142 static VALUE
1143 nurat_marshal_load(VALUE self, VALUE a)
1145 get_dat1(self);
1146 dat->num = RARRAY_PTR(a)[0];
1147 dat->den = RARRAY_PTR(a)[1];
1149 if (f_zero_p(dat->den))
1150 rb_raise(rb_eZeroDivError, "devided by zero");
1152 return self;
1155 /* --- */
1157 VALUE
1158 rb_gcd(VALUE self, VALUE other)
1160 nurat_int_check(other);
1161 return f_gcd(self, other);
1164 VALUE
1165 rb_lcm(VALUE self, VALUE other)
1167 nurat_int_check(other);
1168 return f_lcm(self, other);
1171 VALUE
1172 rb_gcdlcm(VALUE self, VALUE other)
1174 nurat_int_check(other);
1175 return rb_assoc_new(f_gcd(self, other), f_lcm(self, other));
1178 VALUE
1179 rb_rational_raw(VALUE x, VALUE y)
1181 return nurat_s_new_internal(rb_cRational, x, y);
1184 VALUE
1185 rb_rational_new(VALUE x, VALUE y)
1187 return nurat_s_canonicalize_internal(rb_cRational, x, y);
1190 static VALUE nurat_s_convert(int argc, VALUE *argv, VALUE klass);
1192 VALUE
1193 rb_Rational(VALUE x, VALUE y)
1195 VALUE a[2];
1196 a[0] = x;
1197 a[1] = y;
1198 return nurat_s_convert(2, a, rb_cRational);
1201 static VALUE
1202 nilclass_to_r(VALUE self)
1204 return rb_rational_new1(INT2FIX(0));
1207 static VALUE
1208 integer_to_r(VALUE self)
1210 return rb_rational_new1(self);
1213 static VALUE
1214 float_decode(VALUE self)
1216 double f;
1217 int n;
1219 f = frexp(RFLOAT_VALUE(self), &n);
1220 f = ldexp(f, DBL_MANT_DIG);
1221 n -= DBL_MANT_DIG;
1222 return rb_assoc_new(f_to_i(rb_float_new(f)), INT2FIX(n));
1225 static VALUE
1226 float_to_r(VALUE self)
1228 VALUE a = float_decode(self);
1229 return f_mul(RARRAY_PTR(a)[0],
1230 f_expt(INT2FIX(FLT_RADIX), RARRAY_PTR(a)[1]));
1233 static VALUE rat_pat, an_e_pat, a_dot_pat, underscores_pat, an_underscore;
1235 #define DIGITS "(?:\\d(?:_\\d|\\d)*)"
1236 #define NUMERATOR "(?:" DIGITS "?\\.)?" DIGITS "(?:[eE][-+]?" DIGITS ")?"
1237 #define DENOMINATOR "[-+]?" DIGITS
1238 #define PATTERN "\\A([-+])?(" NUMERATOR ")(?:\\/(" DENOMINATOR "))?"
1240 static void
1241 make_patterns(void)
1243 static const char rat_pat_source[] = PATTERN;
1244 static const char an_e_pat_source[] = "[eE]";
1245 static const char a_dot_pat_source[] = "\\.";
1246 static const char underscores_pat_source[] = "_+";
1248 if (rat_pat) return;
1250 rat_pat = rb_reg_new(rat_pat_source, sizeof rat_pat_source - 1, 0);
1251 rb_global_variable(&rat_pat);
1253 an_e_pat = rb_reg_new(an_e_pat_source, sizeof an_e_pat_source - 1, 0);
1254 rb_global_variable(&an_e_pat);
1256 a_dot_pat = rb_reg_new(a_dot_pat_source, sizeof a_dot_pat_source - 1, 0);
1257 rb_global_variable(&a_dot_pat);
1259 underscores_pat = rb_reg_new(underscores_pat_source,
1260 sizeof underscores_pat_source - 1, 0);
1261 rb_global_variable(&underscores_pat);
1263 an_underscore = rb_str_new2("_");
1264 rb_global_variable(&an_underscore);
1267 #define id_strip rb_intern("strip")
1268 #define f_strip(x) rb_funcall(x, id_strip, 0)
1270 #define id_match rb_intern("match")
1271 #define f_match(x,y) rb_funcall(x, id_match, 1, y)
1273 #define id_aref rb_intern("[]")
1274 #define f_aref(x,y) rb_funcall(x, id_aref, 1, y)
1276 #define id_post_match rb_intern("post_match")
1277 #define f_post_match(x) rb_funcall(x, id_post_match, 0)
1279 #define id_split rb_intern("split")
1280 #define f_split(x,y) rb_funcall(x, id_split, 1, y)
1282 #include <ctype.h>
1284 static VALUE
1285 string_to_r_internal(VALUE self)
1287 VALUE s, m;
1289 s = f_strip(self);
1291 if (RSTRING_LEN(s) == 0)
1292 return rb_assoc_new(Qnil, self);
1294 m = f_match(rat_pat, s);
1296 if (!NIL_P(m)) {
1297 VALUE v, ifp, exp, ip, fp;
1298 VALUE si = f_aref(m, INT2FIX(1));
1299 VALUE nu = f_aref(m, INT2FIX(2));
1300 VALUE de = f_aref(m, INT2FIX(3));
1301 VALUE re = f_post_match(m);
1304 VALUE a;
1306 a = f_split(nu, an_e_pat);
1307 ifp = RARRAY_PTR(a)[0];
1308 if (RARRAY_LEN(a) != 2)
1309 exp = Qnil;
1310 else
1311 exp = RARRAY_PTR(a)[1];
1313 a = f_split(ifp, a_dot_pat);
1314 ip = RARRAY_PTR(a)[0];
1315 if (RARRAY_LEN(a) != 2)
1316 fp = Qnil;
1317 else
1318 fp = RARRAY_PTR(a)[1];
1321 v = rb_rational_new1(f_to_i(ip));
1323 if (!NIL_P(fp)) {
1324 char *p = StringValuePtr(fp);
1325 long count = 0;
1326 VALUE l;
1328 while (*p) {
1329 if (rb_isdigit(*p))
1330 count++;
1331 p++;
1334 l = f_expt(INT2FIX(10), LONG2NUM(count));
1335 v = f_mul(v, l);
1336 v = f_add(v, f_to_i(fp));
1337 v = f_div(v, l);
1339 if (!NIL_P(exp))
1340 v = f_mul(v, f_expt(INT2FIX(10), f_to_i(exp)));
1341 if (!NIL_P(si) && *StringValuePtr(si) == '-')
1342 v = f_negate(v);
1343 if (!NIL_P(de))
1344 v = f_div(v, f_to_i(de));
1346 return rb_assoc_new(v, re);
1348 return rb_assoc_new(Qnil, self);
1351 static VALUE
1352 string_to_r_strict(VALUE self)
1354 VALUE a = string_to_r_internal(self);
1355 if (NIL_P(RARRAY_PTR(a)[0]) || RSTRING_LEN(RARRAY_PTR(a)[1]) > 0) {
1356 VALUE s = f_inspect(self);
1357 rb_raise(rb_eArgError, "invalid value for Rational: %s",
1358 StringValuePtr(s));
1360 return RARRAY_PTR(a)[0];
1363 #define id_gsub rb_intern("gsub")
1364 #define f_gsub(x,y,z) rb_funcall(x, id_gsub, 2, y, z)
1366 static VALUE
1367 string_to_r(VALUE self)
1369 VALUE s, a, backref;
1371 backref = rb_backref_get();
1372 rb_match_busy(backref);
1374 s = f_gsub(self, underscores_pat, an_underscore);
1375 a = string_to_r_internal(s);
1377 rb_backref_set(backref);
1379 if (!NIL_P(RARRAY_PTR(a)[0]))
1380 return RARRAY_PTR(a)[0];
1381 return rb_rational_new1(INT2FIX(0));
1384 #define id_to_r rb_intern("to_r")
1385 #define f_to_r(x) rb_funcall(x, id_to_r, 0)
1387 static VALUE
1388 nurat_s_convert(int argc, VALUE *argv, VALUE klass)
1390 VALUE a1, a2, backref;
1392 rb_scan_args(argc, argv, "02", &a1, &a2);
1394 switch (TYPE(a1)) {
1395 case T_COMPLEX:
1396 if (k_float_p(RCOMPLEX(a1)->image) || !f_zero_p(RCOMPLEX(a1)->image)) {
1397 VALUE s = f_to_s(a1);
1398 rb_raise(rb_eRangeError, "can't accept %s",
1399 StringValuePtr(s));
1401 a1 = RCOMPLEX(a1)->real;
1404 switch (TYPE(a2)) {
1405 case T_COMPLEX:
1406 if (k_float_p(RCOMPLEX(a2)->image) || !f_zero_p(RCOMPLEX(a2)->image)) {
1407 VALUE s = f_to_s(a2);
1408 rb_raise(rb_eRangeError, "can't accept %s",
1409 StringValuePtr(s));
1411 a2 = RCOMPLEX(a2)->real;
1414 backref = rb_backref_get();
1415 rb_match_busy(backref);
1417 switch (TYPE(a1)) {
1418 case T_FIXNUM:
1419 case T_BIGNUM:
1420 break;
1421 case T_FLOAT:
1422 a1 = f_to_r(a1);
1423 break;
1424 case T_STRING:
1425 a1 = string_to_r_strict(a1);
1426 break;
1429 switch (TYPE(a2)) {
1430 case T_FIXNUM:
1431 case T_BIGNUM:
1432 break;
1433 case T_FLOAT:
1434 a2 = f_to_r(a2);
1435 break;
1436 case T_STRING:
1437 a2 = string_to_r_strict(a2);
1438 break;
1441 rb_backref_set(backref);
1443 switch (TYPE(a1)) {
1444 case T_RATIONAL:
1445 if (NIL_P(a2) || f_zero_p(a2))
1446 return a1;
1447 return f_div(a1, a2);
1450 switch (TYPE(a2)) {
1451 case T_RATIONAL:
1452 return f_div(a1, a2);
1456 VALUE argv2[2];
1457 argv2[0] = a1;
1458 argv2[1] = a2;
1459 return nurat_s_new(argc, argv2, klass);
1463 static VALUE
1464 nurat_s_induced_from(VALUE klass, VALUE n)
1466 return f_to_r(n);
1469 void
1470 Init_Rational(void)
1472 #undef rb_intern
1474 assert(fprintf(stderr, "assert() is now active\n"));
1476 id_Unify = rb_intern("Unify");
1477 id_abs = rb_intern("abs");
1478 id_cmp = rb_intern("<=>");
1479 id_convert = rb_intern("convert");
1480 id_equal_p = rb_intern("==");
1481 id_expt = rb_intern("**");
1482 id_floor = rb_intern("floor");
1483 id_format = rb_intern("format");
1484 id_idiv = rb_intern("div");
1485 id_inspect = rb_intern("inspect");
1486 id_negate = rb_intern("-@");
1487 id_new = rb_intern("new");
1488 id_new_bang = rb_intern("new!");
1489 id_to_f = rb_intern("to_f");
1490 id_to_i = rb_intern("to_i");
1491 id_to_s = rb_intern("to_s");
1492 id_truncate = rb_intern("truncate");
1494 rb_cRational = rb_define_class(RATIONAL_NAME, rb_cNumeric);
1496 rb_define_alloc_func(rb_cRational, nurat_s_alloc);
1497 rb_funcall(rb_cRational, rb_intern("private_class_method"), 1,
1498 ID2SYM(rb_intern("allocate")));
1500 rb_define_singleton_method(rb_cRational, "new!", nurat_s_new_bang, -1);
1501 rb_funcall(rb_cRational, rb_intern("private_class_method"), 1,
1502 ID2SYM(rb_intern("new!")));
1504 rb_define_singleton_method(rb_cRational, "new", nurat_s_new, -1);
1505 rb_funcall(rb_cRational, rb_intern("private_class_method"), 1,
1506 ID2SYM(rb_intern("new")));
1508 rb_define_global_function(RATIONAL_NAME, nurat_f_rational, -1);
1510 rb_define_method(rb_cRational, "numerator", nurat_numerator, 0);
1511 rb_define_method(rb_cRational, "denominator", nurat_denominator, 0);
1513 rb_define_method(rb_cRational, "+", nurat_add, 1);
1514 rb_define_method(rb_cRational, "-", nurat_sub, 1);
1515 rb_define_method(rb_cRational, "*", nurat_mul, 1);
1516 rb_define_method(rb_cRational, "/", nurat_div, 1);
1517 rb_define_method(rb_cRational, "quo", nurat_div, 1);
1518 rb_define_method(rb_cRational, "fdiv", nurat_fdiv, 1);
1519 rb_define_method(rb_cRational, "**", nurat_expt, 1);
1521 rb_define_method(rb_cRational, "<=>", nurat_cmp, 1);
1522 rb_define_method(rb_cRational, "==", nurat_equal_p, 1);
1523 rb_define_method(rb_cRational, "coerce", nurat_coerce, 1);
1525 rb_define_method(rb_cRational, "div", nurat_idiv, 1);
1526 #if NUBY
1527 rb_define_method(rb_cRational, "//", nurat_idiv, 1);
1528 #endif
1529 rb_define_method(rb_cRational, "modulo", nurat_mod, 1);
1530 rb_define_method(rb_cRational, "%", nurat_mod, 1);
1531 rb_define_method(rb_cRational, "divmod", nurat_divmod, 1);
1533 #if 0
1534 rb_define_method(rb_cRational, "quot", nurat_quot, 1);
1535 #endif
1536 rb_define_method(rb_cRational, "remainder", nurat_rem, 1);
1537 #if 0
1538 rb_define_method(rb_cRational, "quotrem", nurat_quotrem, 1);
1539 #endif
1541 rb_define_method(rb_cRational, "abs", nurat_abs, 0);
1543 #if 0
1544 rb_define_method(rb_cRational, "rational?", nurat_true, 0);
1545 rb_define_method(rb_cRational, "exact?", nurat_true, 0);
1546 #endif
1548 rb_define_method(rb_cRational, "floor", nurat_floor, 0);
1549 rb_define_method(rb_cRational, "ceil", nurat_ceil, 0);
1550 rb_define_method(rb_cRational, "truncate", nurat_truncate, 0);
1551 rb_define_method(rb_cRational, "round", nurat_round, 0);
1553 rb_define_method(rb_cRational, "to_i", nurat_truncate, 0);
1554 rb_define_method(rb_cRational, "to_f", nurat_to_f, 0);
1555 rb_define_method(rb_cRational, "to_r", nurat_to_r, 0);
1557 rb_define_method(rb_cRational, "hash", nurat_hash, 0);
1559 rb_define_method(rb_cRational, "to_s", nurat_to_s, 0);
1560 rb_define_method(rb_cRational, "inspect", nurat_inspect, 0);
1562 rb_define_method(rb_cRational, "marshal_dump", nurat_marshal_dump, 0);
1563 rb_define_method(rb_cRational, "marshal_load", nurat_marshal_load, 1);
1565 /* --- */
1567 rb_define_method(rb_cInteger, "gcd", rb_gcd, 1);
1568 rb_define_method(rb_cInteger, "lcm", rb_lcm, 1);
1569 rb_define_method(rb_cInteger, "gcdlcm", rb_gcdlcm, 1);
1571 rb_define_method(rb_cNilClass, "to_r", nilclass_to_r, 0);
1572 rb_define_method(rb_cInteger, "to_r", integer_to_r, 0);
1573 rb_define_method(rb_cFloat, "to_r", float_to_r, 0);
1575 make_patterns();
1577 rb_define_method(rb_cString, "to_r", string_to_r, 0);
1579 rb_define_singleton_method(rb_cRational, "convert", nurat_s_convert, -1);
1580 rb_funcall(rb_cRational, rb_intern("private_class_method"), 1,
1581 ID2SYM(rb_intern("convert")));
1583 rb_include_module(rb_cRational, rb_mPrecision);
1584 rb_define_singleton_method(rb_cRational, "induced_from",
1585 nurat_s_induced_from, 1);