Ruby  1.9.3p551(2014-11-13revision48407)
rational.c
Go to the documentation of this file.
1 /*
2  rational.c: Coded by Tadayoshi Funaba 2008-2011
3 
4  This implementation is based on Keiju Ishitsuka's Rational library
5  which is written in ruby.
6 */
7 
8 #include "ruby.h"
9 #include "internal.h"
10 #include <math.h>
11 #include <float.h>
12 
13 #ifdef HAVE_IEEEFP_H
14 #include <ieeefp.h>
15 #endif
16 
17 #define NDEBUG
18 #include <assert.h>
19 
20 #define ZERO INT2FIX(0)
21 #define ONE INT2FIX(1)
22 #define TWO INT2FIX(2)
23 
25 
29 
30 #define f_boolcast(x) ((x) ? Qtrue : Qfalse)
31 
32 #define binop(n,op) \
33 inline static VALUE \
34 f_##n(VALUE x, VALUE y)\
35 {\
36  return rb_funcall(x, (op), 1, y);\
37 }
38 
39 #define fun1(n) \
40 inline static VALUE \
41 f_##n(VALUE x)\
42 {\
43  return rb_funcall(x, id_##n, 0);\
44 }
45 
46 #define fun2(n) \
47 inline static VALUE \
48 f_##n(VALUE x, VALUE y)\
49 {\
50  return rb_funcall(x, id_##n, 1, y);\
51 }
52 
53 inline static VALUE
55 {
56  if (FIXNUM_P(y) && FIX2LONG(y) == 0)
57  return x;
58  else if (FIXNUM_P(x) && FIX2LONG(x) == 0)
59  return y;
60  return rb_funcall(x, '+', 1, y);
61 }
62 
63 inline static VALUE
65 {
66  if (FIXNUM_P(x) && FIXNUM_P(y)) {
67  long c = FIX2LONG(x) - FIX2LONG(y);
68  if (c > 0)
69  c = 1;
70  else if (c < 0)
71  c = -1;
72  return INT2FIX(c);
73  }
74  return rb_funcall(x, id_cmp, 1, y);
75 }
76 
77 inline static VALUE
79 {
80  if (FIXNUM_P(y) && FIX2LONG(y) == 1)
81  return x;
82  return rb_funcall(x, '/', 1, y);
83 }
84 
85 inline static VALUE
87 {
88  if (FIXNUM_P(x) && FIXNUM_P(y))
89  return f_boolcast(FIX2LONG(x) > FIX2LONG(y));
90  return rb_funcall(x, '>', 1, y);
91 }
92 
93 inline static VALUE
95 {
96  if (FIXNUM_P(x) && FIXNUM_P(y))
97  return f_boolcast(FIX2LONG(x) < FIX2LONG(y));
98  return rb_funcall(x, '<', 1, y);
99 }
100 
101 binop(mod, '%')
102 
103 inline static VALUE
104 f_mul(VALUE x, VALUE y)
105 {
106  if (FIXNUM_P(y)) {
107  long iy = FIX2LONG(y);
108  if (iy == 0) {
109  if (FIXNUM_P(x) || TYPE(x) == T_BIGNUM)
110  return ZERO;
111  }
112  else if (iy == 1)
113  return x;
114  }
115  else if (FIXNUM_P(x)) {
116  long ix = FIX2LONG(x);
117  if (ix == 0) {
118  if (FIXNUM_P(y) || TYPE(y) == T_BIGNUM)
119  return ZERO;
120  }
121  else if (ix == 1)
122  return y;
123  }
124  return rb_funcall(x, '*', 1, y);
125 }
126 
127 inline static VALUE
129 {
130  if (FIXNUM_P(y) && FIX2LONG(y) == 0)
131  return x;
132  return rb_funcall(x, '-', 1, y);
133 }
134 
135 fun1(abs)
136 fun1(floor)
137 fun1(inspect)
138 fun1(integer_p)
139 fun1(negate)
140 
141 inline static VALUE
142 f_to_i(VALUE x)
143 {
144  if (TYPE(x) == T_STRING)
145  return rb_str_to_inum(x, 10, 0);
146  return rb_funcall(x, id_to_i, 0);
147 }
148 inline static VALUE
150 {
151  if (TYPE(x) == T_STRING)
152  return DBL2NUM(rb_str_to_dbl(x, 0));
153  return rb_funcall(x, id_to_f, 0);
154 }
155 
156 fun1(to_s)
157 fun1(truncate)
158 
159 inline static VALUE
160 f_eqeq_p(VALUE x, VALUE y)
161 {
162  if (FIXNUM_P(x) && FIXNUM_P(y))
163  return f_boolcast(FIX2LONG(x) == FIX2LONG(y));
164  return rb_funcall(x, id_eqeq_p, 1, y);
165 }
166 
167 fun2(expt)
168 fun2(fdiv)
169 fun2(idiv)
170 
171 #define f_expt10(x) f_expt(INT2FIX(10), x)
172 
173 inline static VALUE
175 {
176  if (FIXNUM_P(x))
177  return f_boolcast(FIX2LONG(x) < 0);
178  return rb_funcall(x, '<', 1, ZERO);
179 }
180 
181 #define f_positive_p(x) (!f_negative_p(x))
182 
183 inline static VALUE
185 {
186  switch (TYPE(x)) {
187  case T_FIXNUM:
188  return f_boolcast(FIX2LONG(x) == 0);
189  case T_BIGNUM:
190  return Qfalse;
191  case T_RATIONAL:
192  {
193  VALUE num = RRATIONAL(x)->num;
194 
195  return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 0);
196  }
197  }
198  return rb_funcall(x, id_eqeq_p, 1, ZERO);
199 }
200 
201 #define f_nonzero_p(x) (!f_zero_p(x))
202 
203 inline static VALUE
205 {
206  switch (TYPE(x)) {
207  case T_FIXNUM:
208  return f_boolcast(FIX2LONG(x) == 1);
209  case T_BIGNUM:
210  return Qfalse;
211  case T_RATIONAL:
212  {
213  VALUE num = RRATIONAL(x)->num;
214  VALUE den = RRATIONAL(x)->den;
215 
216  return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 1 &&
217  FIXNUM_P(den) && FIX2LONG(den) == 1);
218  }
219  }
220  return rb_funcall(x, id_eqeq_p, 1, ONE);
221 }
222 
223 inline static VALUE
225 {
226  return rb_obj_is_kind_of(x, c);
227 }
228 
229 inline static VALUE
231 {
232  return f_kind_of_p(x, rb_cNumeric);
233 }
234 
235 inline static VALUE
237 {
238  return f_kind_of_p(x, rb_cInteger);
239 }
240 
241 inline static VALUE
243 {
244  return f_kind_of_p(x, rb_cFloat);
245 }
246 
247 inline static VALUE
249 {
250  return f_kind_of_p(x, rb_cRational);
251 }
252 
253 #define k_exact_p(x) (!k_float_p(x))
254 #define k_inexact_p(x) k_float_p(x)
255 
256 #define k_exact_zero_p(x) (k_exact_p(x) && f_zero_p(x))
257 #define k_exact_one_p(x) (k_exact_p(x) && f_one_p(x))
258 
259 #ifndef NDEBUG
260 #define f_gcd f_gcd_orig
261 #endif
262 
263 inline static long
264 i_gcd(long x, long y)
265 {
266  if (x < 0)
267  x = -x;
268  if (y < 0)
269  y = -y;
270 
271  if (x == 0)
272  return y;
273  if (y == 0)
274  return x;
275 
276  while (x > 0) {
277  long t = x;
278  x = y % x;
279  y = t;
280  }
281  return y;
282 }
283 
284 inline static VALUE
286 {
287  VALUE z;
288 
289  if (FIXNUM_P(x) && FIXNUM_P(y))
290  return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
291 
292  if (f_negative_p(x))
293  x = f_negate(x);
294  if (f_negative_p(y))
295  y = f_negate(y);
296 
297  if (f_zero_p(x))
298  return y;
299  if (f_zero_p(y))
300  return x;
301 
302  for (;;) {
303  if (FIXNUM_P(x)) {
304  if (FIX2LONG(x) == 0)
305  return y;
306  if (FIXNUM_P(y))
307  return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
308  }
309  z = x;
310  x = f_mod(y, x);
311  y = z;
312  }
313  /* NOTREACHED */
314 }
315 
316 #ifndef NDEBUG
317 #undef f_gcd
318 
319 inline static VALUE
320 f_gcd(VALUE x, VALUE y)
321 {
322  VALUE r = f_gcd_orig(x, y);
323  if (f_nonzero_p(r)) {
324  assert(f_zero_p(f_mod(x, r)));
325  assert(f_zero_p(f_mod(y, r)));
326  }
327  return r;
328 }
329 #endif
330 
331 inline static VALUE
333 {
334  if (f_zero_p(x) || f_zero_p(y))
335  return ZERO;
336  return f_abs(f_mul(f_div(x, f_gcd(x, y)), y));
337 }
338 
339 #define get_dat1(x) \
340  struct RRational *dat;\
341  dat = ((struct RRational *)(x))
342 
343 #define get_dat2(x,y) \
344  struct RRational *adat, *bdat;\
345  adat = ((struct RRational *)(x));\
346  bdat = ((struct RRational *)(y))
347 
348 inline static VALUE
350 {
351  NEWOBJ(obj, struct RRational);
352  OBJSETUP(obj, klass, T_RATIONAL);
353 
354  obj->num = num;
355  obj->den = den;
356 
357  return (VALUE)obj;
358 }
359 
360 static VALUE
362 {
363  return nurat_s_new_internal(klass, ZERO, ONE);
364 }
365 
366 #define rb_raise_zerodiv() rb_raise(rb_eZeroDivError, "divided by 0")
367 
368 #if 0
369 static VALUE
370 nurat_s_new_bang(int argc, VALUE *argv, VALUE klass)
371 {
372  VALUE num, den;
373 
374  switch (rb_scan_args(argc, argv, "11", &num, &den)) {
375  case 1:
376  if (!k_integer_p(num))
377  num = f_to_i(num);
378  den = ONE;
379  break;
380  default:
381  if (!k_integer_p(num))
382  num = f_to_i(num);
383  if (!k_integer_p(den))
384  den = f_to_i(den);
385 
386  switch (FIX2INT(f_cmp(den, ZERO))) {
387  case -1:
388  num = f_negate(num);
389  den = f_negate(den);
390  break;
391  case 0:
393  break;
394  }
395  break;
396  }
397 
398  return nurat_s_new_internal(klass, num, den);
399 }
400 #endif
401 
402 inline static VALUE
404 {
405  return nurat_s_new_internal(klass, x, ONE);
406 }
407 
408 inline static VALUE
410 {
411  assert(f_positive_p(y));
412  assert(f_nonzero_p(y));
413  return nurat_s_new_internal(klass, x, y);
414 }
415 
416 #ifdef CANONICALIZATION_FOR_MATHN
417 #define CANON
418 #endif
419 
420 #ifdef CANON
421 static int canonicalization = 0;
422 
425 {
426  canonicalization = f;
427 }
428 #endif
429 
430 inline static void
432 {
433  switch (TYPE(num)) {
434  case T_FIXNUM:
435  case T_BIGNUM:
436  break;
437  default:
438  if (!k_numeric_p(num) || !f_integer_p(num))
439  rb_raise(rb_eTypeError, "not an integer");
440  }
441 }
442 
443 inline static VALUE
445 {
446  nurat_int_check(num);
447  if (!k_integer_p(num))
448  num = f_to_i(num);
449  return num;
450 }
451 
452 inline static VALUE
454 {
455  VALUE gcd;
456 
457  switch (FIX2INT(f_cmp(den, ZERO))) {
458  case -1:
459  num = f_negate(num);
460  den = f_negate(den);
461  break;
462  case 0:
464  break;
465  }
466 
467  gcd = f_gcd(num, den);
468  num = f_idiv(num, gcd);
469  den = f_idiv(den, gcd);
470 
471 #ifdef CANON
472  if (f_one_p(den) && canonicalization)
473  return num;
474 #endif
475  return nurat_s_new_internal(klass, num, den);
476 }
477 
478 inline static VALUE
480 {
481  switch (FIX2INT(f_cmp(den, ZERO))) {
482  case -1:
483  num = f_negate(num);
484  den = f_negate(den);
485  break;
486  case 0:
488  break;
489  }
490 
491 #ifdef CANON
492  if (f_one_p(den) && canonicalization)
493  return num;
494 #endif
495  return nurat_s_new_internal(klass, num, den);
496 }
497 
498 static VALUE
500 {
501  VALUE num, den;
502 
503  switch (rb_scan_args(argc, argv, "11", &num, &den)) {
504  case 1:
505  num = nurat_int_value(num);
506  den = ONE;
507  break;
508  default:
509  num = nurat_int_value(num);
510  den = nurat_int_value(den);
511  break;
512  }
513 
514  return nurat_s_canonicalize_internal(klass, num, den);
515 }
516 
517 inline static VALUE
519 {
520  assert(!k_rational_p(x));
521  return nurat_s_canonicalize_internal(klass, x, ONE);
522 }
523 
524 inline static VALUE
526 {
527  assert(!k_rational_p(x));
528  assert(!k_rational_p(y));
529  return nurat_s_canonicalize_internal(klass, x, y);
530 }
531 
532 inline static VALUE
534 {
535  assert(!k_rational_p(x));
537 }
538 
539 inline static VALUE
541 {
542  assert(!k_rational_p(x));
543  assert(!k_rational_p(y));
544  return nurat_s_canonicalize_internal_no_reduce(klass, x, y);
545 }
546 
547 /*
548  * call-seq:
549  * Rational(x[, y]) -> numeric
550  *
551  * Returns x/y;
552  */
553 static VALUE
555 {
556  return rb_funcall2(rb_cRational, id_convert, argc, argv);
557 }
558 
559 /*
560  * call-seq:
561  * rat.numerator -> integer
562  *
563  * Returns the numerator.
564  *
565  * For example:
566  *
567  * Rational(7).numerator #=> 7
568  * Rational(7, 1).numerator #=> 7
569  * Rational(9, -4).numerator #=> -9
570  * Rational(-2, -10).numerator #=> 1
571  */
572 static VALUE
574 {
575  get_dat1(self);
576  return dat->num;
577 }
578 
579 /*
580  * call-seq:
581  * rat.denominator -> integer
582  *
583  * Returns the denominator (always positive).
584  *
585  * For example:
586  *
587  * Rational(7).denominator #=> 1
588  * Rational(7, 1).denominator #=> 1
589  * Rational(9, -4).denominator #=> 4
590  * Rational(-2, -10).denominator #=> 5
591  * rat.numerator.gcd(rat.denominator) #=> 1
592  */
593 static VALUE
595 {
596  get_dat1(self);
597  return dat->den;
598 }
599 
600 #ifndef NDEBUG
601 #define f_imul f_imul_orig
602 #endif
603 
604 inline static VALUE
605 f_imul(long a, long b)
606 {
607  VALUE r;
608  volatile long c;
609 
610  if (a == 0 || b == 0)
611  return ZERO;
612  else if (a == 1)
613  return LONG2NUM(b);
614  else if (b == 1)
615  return LONG2NUM(a);
616 
617  c = a * b;
618  r = LONG2NUM(c);
619  if (NUM2LONG(r) != c || (c / a) != b)
620  r = rb_big_mul(rb_int2big(a), rb_int2big(b));
621  return r;
622 }
623 
624 #ifndef NDEBUG
625 #undef f_imul
626 
627 inline static VALUE
628 f_imul(long x, long y)
629 {
630  VALUE r = f_imul_orig(x, y);
631  assert(f_eqeq_p(r, f_mul(LONG2NUM(x), LONG2NUM(y))));
632  return r;
633 }
634 #endif
635 
636 inline static VALUE
637 f_addsub(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
638 {
639  VALUE num, den;
640 
641  if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
642  FIXNUM_P(bnum) && FIXNUM_P(bden)) {
643  long an = FIX2LONG(anum);
644  long ad = FIX2LONG(aden);
645  long bn = FIX2LONG(bnum);
646  long bd = FIX2LONG(bden);
647  long ig = i_gcd(ad, bd);
648 
649  VALUE g = LONG2NUM(ig);
650  VALUE a = f_imul(an, bd / ig);
651  VALUE b = f_imul(bn, ad / ig);
652  VALUE c;
653 
654  if (k == '+')
655  c = f_add(a, b);
656  else
657  c = f_sub(a, b);
658 
659  b = f_idiv(aden, g);
660  g = f_gcd(c, g);
661  num = f_idiv(c, g);
662  a = f_idiv(bden, g);
663  den = f_mul(a, b);
664  }
665  else {
666  VALUE g = f_gcd(aden, bden);
667  VALUE a = f_mul(anum, f_idiv(bden, g));
668  VALUE b = f_mul(bnum, f_idiv(aden, g));
669  VALUE c;
670 
671  if (k == '+')
672  c = f_add(a, b);
673  else
674  c = f_sub(a, b);
675 
676  b = f_idiv(aden, g);
677  g = f_gcd(c, g);
678  num = f_idiv(c, g);
679  a = f_idiv(bden, g);
680  den = f_mul(a, b);
681  }
682  return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
683 }
684 
685 /*
686  * call-seq:
687  * rat + numeric -> numeric
688  *
689  * Performs addition.
690  *
691  * For example:
692  *
693  * Rational(2, 3) + Rational(2, 3) #=> (4/3)
694  * Rational(900) + Rational(1) #=> (900/1)
695  * Rational(-2, 9) + Rational(-9, 2) #=> (-85/18)
696  * Rational(9, 8) + 4 #=> (41/8)
697  * Rational(20, 9) + 9.8 #=> 12.022222222222222
698  */
699 static VALUE
700 nurat_add(VALUE self, VALUE other)
701 {
702  switch (TYPE(other)) {
703  case T_FIXNUM:
704  case T_BIGNUM:
705  {
706  get_dat1(self);
707 
708  return f_addsub(self,
709  dat->num, dat->den,
710  other, ONE, '+');
711  }
712  case T_FLOAT:
713  return f_add(f_to_f(self), other);
714  case T_RATIONAL:
715  {
716  get_dat2(self, other);
717 
718  return f_addsub(self,
719  adat->num, adat->den,
720  bdat->num, bdat->den, '+');
721  }
722  default:
723  return rb_num_coerce_bin(self, other, '+');
724  }
725 }
726 
727 /*
728  * call-seq:
729  * rat - numeric -> numeric
730  *
731  * Performs subtraction.
732  *
733  * For example:
734  *
735  * Rational(2, 3) - Rational(2, 3) #=> (0/1)
736  * Rational(900) - Rational(1) #=> (899/1)
737  * Rational(-2, 9) - Rational(-9, 2) #=> (77/18)
738  * Rational(9, 8) - 4 #=> (23/8)
739  * Rational(20, 9) - 9.8 #=> -7.577777777777778
740  */
741 static VALUE
742 nurat_sub(VALUE self, VALUE other)
743 {
744  switch (TYPE(other)) {
745  case T_FIXNUM:
746  case T_BIGNUM:
747  {
748  get_dat1(self);
749 
750  return f_addsub(self,
751  dat->num, dat->den,
752  other, ONE, '-');
753  }
754  case T_FLOAT:
755  return f_sub(f_to_f(self), other);
756  case T_RATIONAL:
757  {
758  get_dat2(self, other);
759 
760  return f_addsub(self,
761  adat->num, adat->den,
762  bdat->num, bdat->den, '-');
763  }
764  default:
765  return rb_num_coerce_bin(self, other, '-');
766  }
767 }
768 
769 inline static VALUE
770 f_muldiv(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
771 {
772  VALUE num, den;
773 
774  if (k == '/') {
775  VALUE t;
776 
777  if (f_negative_p(bnum)) {
778  anum = f_negate(anum);
779  bnum = f_negate(bnum);
780  }
781  t = bnum;
782  bnum = bden;
783  bden = t;
784  }
785 
786  if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
787  FIXNUM_P(bnum) && FIXNUM_P(bden)) {
788  long an = FIX2LONG(anum);
789  long ad = FIX2LONG(aden);
790  long bn = FIX2LONG(bnum);
791  long bd = FIX2LONG(bden);
792  long g1 = i_gcd(an, bd);
793  long g2 = i_gcd(ad, bn);
794 
795  num = f_imul(an / g1, bn / g2);
796  den = f_imul(ad / g2, bd / g1);
797  }
798  else {
799  VALUE g1 = f_gcd(anum, bden);
800  VALUE g2 = f_gcd(aden, bnum);
801 
802  num = f_mul(f_idiv(anum, g1), f_idiv(bnum, g2));
803  den = f_mul(f_idiv(aden, g2), f_idiv(bden, g1));
804  }
805  return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
806 }
807 
808 /*
809  * call-seq:
810  * rat * numeric -> numeric
811  *
812  * Performs multiplication.
813  *
814  * For example:
815  *
816  * Rational(2, 3) * Rational(2, 3) #=> (4/9)
817  * Rational(900) * Rational(1) #=> (900/1)
818  * Rational(-2, 9) * Rational(-9, 2) #=> (1/1)
819  * Rational(9, 8) * 4 #=> (9/2)
820  * Rational(20, 9) * 9.8 #=> 21.77777777777778
821  */
822 static VALUE
823 nurat_mul(VALUE self, VALUE other)
824 {
825  switch (TYPE(other)) {
826  case T_FIXNUM:
827  case T_BIGNUM:
828  {
829  get_dat1(self);
830 
831  return f_muldiv(self,
832  dat->num, dat->den,
833  other, ONE, '*');
834  }
835  case T_FLOAT:
836  return f_mul(f_to_f(self), other);
837  case T_RATIONAL:
838  {
839  get_dat2(self, other);
840 
841  return f_muldiv(self,
842  adat->num, adat->den,
843  bdat->num, bdat->den, '*');
844  }
845  default:
846  return rb_num_coerce_bin(self, other, '*');
847  }
848 }
849 
850 /*
851  * call-seq:
852  * rat / numeric -> numeric
853  * rat.quo(numeric) -> numeric
854  *
855  * Performs division.
856  *
857  * For example:
858  *
859  * Rational(2, 3) / Rational(2, 3) #=> (1/1)
860  * Rational(900) / Rational(1) #=> (900/1)
861  * Rational(-2, 9) / Rational(-9, 2) #=> (4/81)
862  * Rational(9, 8) / 4 #=> (9/32)
863  * Rational(20, 9) / 9.8 #=> 0.22675736961451246
864  */
865 static VALUE
866 nurat_div(VALUE self, VALUE other)
867 {
868  switch (TYPE(other)) {
869  case T_FIXNUM:
870  case T_BIGNUM:
871  if (f_zero_p(other))
873  {
874  get_dat1(self);
875 
876  return f_muldiv(self,
877  dat->num, dat->den,
878  other, ONE, '/');
879  }
880  case T_FLOAT:
881  {
882  double x = RFLOAT_VALUE(other), den;
883  get_dat1(self);
884 
885  if (isnan(x)) return DBL2NUM(NAN);
886  if (isinf(x)) return INT2FIX(0);
887  if (x != 0.0 && modf(x, &den) == 0.0) {
888  return rb_rational_raw2(dat->num, f_mul(rb_dbl2big(den), dat->den));
889  }
890  }
891  return rb_funcall(f_to_f(self), '/', 1, other);
892  case T_RATIONAL:
893  if (f_zero_p(other))
895  {
896  get_dat2(self, other);
897 
898  if (f_one_p(self))
899  return f_rational_new_no_reduce2(CLASS_OF(self),
900  bdat->den, bdat->num);
901 
902  return f_muldiv(self,
903  adat->num, adat->den,
904  bdat->num, bdat->den, '/');
905  }
906  default:
907  return rb_num_coerce_bin(self, other, '/');
908  }
909 }
910 
911 /*
912  * call-seq:
913  * rat.fdiv(numeric) -> float
914  *
915  * Performs division and returns the value as a float.
916  *
917  * For example:
918  *
919  * Rational(2, 3).fdiv(1) #=> 0.6666666666666666
920  * Rational(2, 3).fdiv(0.5) #=> 1.3333333333333333
921  * Rational(2).fdiv(3) #=> 0.6666666666666666
922  */
923 static VALUE
924 nurat_fdiv(VALUE self, VALUE other)
925 {
926  if (f_zero_p(other))
927  return f_div(self, f_to_f(other));
928  return f_to_f(f_div(self, other));
929 }
930 
931 /*
932  * call-seq:
933  * rat ** numeric -> numeric
934  *
935  * Performs exponentiation.
936  *
937  * For example:
938  *
939  * Rational(2) ** Rational(3) #=> (8/1)
940  * Rational(10) ** -2 #=> (1/100)
941  * Rational(10) ** -2.0 #=> 0.01
942  * Rational(-4) ** Rational(1,2) #=> (1.2246063538223773e-16+2.0i)
943  * Rational(1, 2) ** 0 #=> (1/1)
944  * Rational(1, 2) ** 0.0 #=> 1.0
945  */
946 static VALUE
947 nurat_expt(VALUE self, VALUE other)
948 {
949  if (k_numeric_p(other) && k_exact_zero_p(other))
950  return f_rational_new_bang1(CLASS_OF(self), ONE);
951 
952  if (k_rational_p(other)) {
953  get_dat1(other);
954 
955  if (f_one_p(dat->den))
956  other = dat->num; /* c14n */
957  }
958 
959  switch (TYPE(other)) {
960  case T_FIXNUM:
961  {
962  VALUE num, den;
963 
964  get_dat1(self);
965 
966  switch (FIX2INT(f_cmp(other, ZERO))) {
967  case 1:
968  num = f_expt(dat->num, other);
969  den = f_expt(dat->den, other);
970  break;
971  case -1:
972  num = f_expt(dat->den, f_negate(other));
973  den = f_expt(dat->num, f_negate(other));
974  break;
975  default:
976  num = ONE;
977  den = ONE;
978  break;
979  }
980  return f_rational_new2(CLASS_OF(self), num, den);
981  }
982  case T_BIGNUM:
983  rb_warn("in a**b, b may be too big");
984  /* fall through */
985  case T_FLOAT:
986  case T_RATIONAL:
987  return f_expt(f_to_f(self), other);
988  default:
989  return rb_num_coerce_bin(self, other, id_expt);
990  }
991 }
992 
993 /*
994  * call-seq:
995  * rat <=> numeric -> -1, 0, +1 or nil
996  *
997  * Performs comparison and returns -1, 0, or +1.
998  *
999  * For example:
1000  *
1001  * Rational(2, 3) <=> Rational(2, 3) #=> 0
1002  * Rational(5) <=> 5 #=> 0
1003  * Rational(2,3) <=> Rational(1,3) #=> 1
1004  * Rational(1,3) <=> 1 #=> -1
1005  * Rational(1,3) <=> 0.3 #=> 1
1006  */
1007 static VALUE
1008 nurat_cmp(VALUE self, VALUE other)
1009 {
1010  switch (TYPE(other)) {
1011  case T_FIXNUM:
1012  case T_BIGNUM:
1013  {
1014  get_dat1(self);
1015 
1016  if (FIXNUM_P(dat->den) && FIX2LONG(dat->den) == 1)
1017  return f_cmp(dat->num, other); /* c14n */
1018  return f_cmp(self, f_rational_new_bang1(CLASS_OF(self), other));
1019  }
1020  case T_FLOAT:
1021  return f_cmp(f_to_f(self), other);
1022  case T_RATIONAL:
1023  {
1024  VALUE num1, num2;
1025 
1026  get_dat2(self, other);
1027 
1028  if (FIXNUM_P(adat->num) && FIXNUM_P(adat->den) &&
1029  FIXNUM_P(bdat->num) && FIXNUM_P(bdat->den)) {
1030  num1 = f_imul(FIX2LONG(adat->num), FIX2LONG(bdat->den));
1031  num2 = f_imul(FIX2LONG(bdat->num), FIX2LONG(adat->den));
1032  }
1033  else {
1034  num1 = f_mul(adat->num, bdat->den);
1035  num2 = f_mul(bdat->num, adat->den);
1036  }
1037  return f_cmp(f_sub(num1, num2), ZERO);
1038  }
1039  default:
1040  return rb_num_coerce_cmp(self, other, id_cmp);
1041  }
1042 }
1043 
1044 /*
1045  * call-seq:
1046  * rat == object -> true or false
1047  *
1048  * Returns true if rat equals object numerically.
1049  *
1050  * For example:
1051  *
1052  * Rational(2, 3) == Rational(2, 3) #=> true
1053  * Rational(5) == 5 #=> true
1054  * Rational(0) == 0.0 #=> true
1055  * Rational('1/3') == 0.33 #=> false
1056  * Rational('1/2') == '1/2' #=> false
1057  */
1058 static VALUE
1060 {
1061  switch (TYPE(other)) {
1062  case T_FIXNUM:
1063  case T_BIGNUM:
1064  {
1065  get_dat1(self);
1066 
1067  if (f_zero_p(dat->num) && f_zero_p(other))
1068  return Qtrue;
1069 
1070  if (!FIXNUM_P(dat->den))
1071  return Qfalse;
1072  if (FIX2LONG(dat->den) != 1)
1073  return Qfalse;
1074  if (f_eqeq_p(dat->num, other))
1075  return Qtrue;
1076  return Qfalse;
1077  }
1078  case T_FLOAT:
1079  return f_eqeq_p(f_to_f(self), other);
1080  case T_RATIONAL:
1081  {
1082  get_dat2(self, other);
1083 
1084  if (f_zero_p(adat->num) && f_zero_p(bdat->num))
1085  return Qtrue;
1086 
1087  return f_boolcast(f_eqeq_p(adat->num, bdat->num) &&
1088  f_eqeq_p(adat->den, bdat->den));
1089  }
1090  default:
1091  return f_eqeq_p(other, self);
1092  }
1093 }
1094 
1095 /* :nodoc: */
1096 static VALUE
1098 {
1099  switch (TYPE(other)) {
1100  case T_FIXNUM:
1101  case T_BIGNUM:
1102  return rb_assoc_new(f_rational_new_bang1(CLASS_OF(self), other), self);
1103  case T_FLOAT:
1104  return rb_assoc_new(other, f_to_f(self));
1105  case T_RATIONAL:
1106  return rb_assoc_new(other, self);
1107  case T_COMPLEX:
1108  if (k_exact_zero_p(RCOMPLEX(other)->imag))
1110  (CLASS_OF(self), RCOMPLEX(other)->real), self);
1111  else
1112  return rb_assoc_new(other, rb_Complex(self, INT2FIX(0)));
1113  }
1114 
1115  rb_raise(rb_eTypeError, "%s can't be coerced into %s",
1116  rb_obj_classname(other), rb_obj_classname(self));
1117  return Qnil;
1118 }
1119 
1120 #if 0
1121 /* :nodoc: */
1122 static VALUE
1123 nurat_idiv(VALUE self, VALUE other)
1124 {
1125  return f_idiv(self, other);
1126 }
1127 
1128 /* :nodoc: */
1129 static VALUE
1130 nurat_quot(VALUE self, VALUE other)
1131 {
1132  return f_truncate(f_div(self, other));
1133 }
1134 
1135 /* :nodoc: */
1136 static VALUE
1137 nurat_quotrem(VALUE self, VALUE other)
1138 {
1139  VALUE val = f_truncate(f_div(self, other));
1140  return rb_assoc_new(val, f_sub(self, f_mul(other, val)));
1141 }
1142 #endif
1143 
1144 #if 0
1145 /* :nodoc: */
1146 static VALUE
1147 nurat_true(VALUE self)
1148 {
1149  return Qtrue;
1150 }
1151 #endif
1152 
1153 static VALUE
1155 {
1156  get_dat1(self);
1157  return f_idiv(dat->num, dat->den);
1158 }
1159 
1160 static VALUE
1162 {
1163  get_dat1(self);
1164  return f_negate(f_idiv(f_negate(dat->num), dat->den));
1165 }
1166 
1167 /*
1168  * call-seq:
1169  * rat.to_i -> integer
1170  *
1171  * Returns the truncated value as an integer.
1172  *
1173  * Equivalent to
1174  * rat.truncate.
1175  *
1176  * For example:
1177  *
1178  * Rational(2, 3).to_i #=> 0
1179  * Rational(3).to_i #=> 3
1180  * Rational(300.6).to_i #=> 300
1181  * Rational(98,71).to_i #=> 1
1182  * Rational(-30,2).to_i #=> -15
1183  */
1184 static VALUE
1186 {
1187  get_dat1(self);
1188  if (f_negative_p(dat->num))
1189  return f_negate(f_idiv(f_negate(dat->num), dat->den));
1190  return f_idiv(dat->num, dat->den);
1191 }
1192 
1193 static VALUE
1195 {
1196  VALUE num, den, neg;
1197 
1198  get_dat1(self);
1199 
1200  num = dat->num;
1201  den = dat->den;
1202  neg = f_negative_p(num);
1203 
1204  if (neg)
1205  num = f_negate(num);
1206 
1207  num = f_add(f_mul(num, TWO), den);
1208  den = f_mul(den, TWO);
1209  num = f_idiv(num, den);
1210 
1211  if (neg)
1212  num = f_negate(num);
1213 
1214  return num;
1215 }
1216 
1217 static VALUE
1219 {
1220  VALUE n, b, s;
1221 
1222  if (argc == 0)
1223  return (*func)(self);
1224 
1225  rb_scan_args(argc, argv, "01", &n);
1226 
1227  if (!k_integer_p(n))
1228  rb_raise(rb_eTypeError, "not an integer");
1229 
1230  b = f_expt10(n);
1231  s = f_mul(self, b);
1232 
1233  if (!k_rational_p(s)) {
1234  s = f_rational_new_bang1(CLASS_OF(self), s);
1235  }
1236 
1237  s = (*func)(s);
1238 
1239  s = f_div(f_rational_new_bang1(CLASS_OF(self), s), b);
1240 
1241  if (f_lt_p(n, ONE))
1242  s = f_to_i(s);
1243 
1244  return s;
1245 }
1246 
1247 /*
1248  * call-seq:
1249  * rat.floor -> integer
1250  * rat.floor(precision=0) -> rational
1251  *
1252  * Returns the truncated value (toward negative infinity).
1253  *
1254  * For example:
1255  *
1256  * Rational(3).floor #=> 3
1257  * Rational(2, 3).floor #=> 0
1258  * Rational(-3, 2).floor #=> -1
1259  *
1260  * decimal - 1 2 3 . 4 5 6
1261  * ^ ^ ^ ^ ^ ^
1262  * precision -3 -2 -1 0 +1 +2
1263  *
1264  * '%f' % Rational('-123.456').floor(+1) #=> "-123.500000"
1265  * '%f' % Rational('-123.456').floor(-1) #=> "-130.000000"
1266  */
1267 static VALUE
1268 nurat_floor_n(int argc, VALUE *argv, VALUE self)
1269 {
1270  return f_round_common(argc, argv, self, nurat_floor);
1271 }
1272 
1273 /*
1274  * call-seq:
1275  * rat.ceil -> integer
1276  * rat.ceil(precision=0) -> rational
1277  *
1278  * Returns the truncated value (toward positive infinity).
1279  *
1280  * For example:
1281  *
1282  * Rational(3).ceil #=> 3
1283  * Rational(2, 3).ceil #=> 1
1284  * Rational(-3, 2).ceil #=> -1
1285  *
1286  * decimal - 1 2 3 . 4 5 6
1287  * ^ ^ ^ ^ ^ ^
1288  * precision -3 -2 -1 0 +1 +2
1289  *
1290  * '%f' % Rational('-123.456').ceil(+1) #=> "-123.400000"
1291  * '%f' % Rational('-123.456').ceil(-1) #=> "-120.000000"
1292  */
1293 static VALUE
1294 nurat_ceil_n(int argc, VALUE *argv, VALUE self)
1295 {
1296  return f_round_common(argc, argv, self, nurat_ceil);
1297 }
1298 
1299 /*
1300  * call-seq:
1301  * rat.truncate -> integer
1302  * rat.truncate(precision=0) -> rational
1303  *
1304  * Returns the truncated value (toward zero).
1305  *
1306  * For example:
1307  *
1308  * Rational(3).truncate #=> 3
1309  * Rational(2, 3).truncate #=> 0
1310  * Rational(-3, 2).truncate #=> -1
1311  *
1312  * decimal - 1 2 3 . 4 5 6
1313  * ^ ^ ^ ^ ^ ^
1314  * precision -3 -2 -1 0 +1 +2
1315  *
1316  * '%f' % Rational('-123.456').truncate(+1) #=> "-123.400000"
1317  * '%f' % Rational('-123.456').truncate(-1) #=> "-120.000000"
1318  */
1319 static VALUE
1320 nurat_truncate_n(int argc, VALUE *argv, VALUE self)
1321 {
1322  return f_round_common(argc, argv, self, nurat_truncate);
1323 }
1324 
1325 /*
1326  * call-seq:
1327  * rat.round -> integer
1328  * rat.round(precision=0) -> rational
1329  *
1330  * Returns the truncated value (toward the nearest integer;
1331  * 0.5 => 1; -0.5 => -1).
1332  *
1333  * For example:
1334  *
1335  * Rational(3).round #=> 3
1336  * Rational(2, 3).round #=> 1
1337  * Rational(-3, 2).round #=> -2
1338  *
1339  * decimal - 1 2 3 . 4 5 6
1340  * ^ ^ ^ ^ ^ ^
1341  * precision -3 -2 -1 0 +1 +2
1342  *
1343  * '%f' % Rational('-123.456').round(+1) #=> "-123.500000"
1344  * '%f' % Rational('-123.456').round(-1) #=> "-120.000000"
1345  */
1346 static VALUE
1347 nurat_round_n(int argc, VALUE *argv, VALUE self)
1348 {
1349  return f_round_common(argc, argv, self, nurat_round);
1350 }
1351 
1352 /*
1353  * call-seq:
1354  * rat.to_f -> float
1355  *
1356  * Return the value as a float.
1357  *
1358  * For example:
1359  *
1360  * Rational(2).to_f #=> 2.0
1361  * Rational(9, 4).to_f #=> 2.25
1362  * Rational(-3, 4).to_f #=> -0.75
1363  * Rational(20, 3).to_f #=> 6.666666666666667
1364  */
1365 static VALUE
1367 {
1368  get_dat1(self);
1369  return f_fdiv(dat->num, dat->den);
1370 }
1371 
1372 /*
1373  * call-seq:
1374  * rat.to_r -> self
1375  *
1376  * Returns self.
1377  *
1378  * For example:
1379  *
1380  * Rational(2).to_r #=> (2/1)
1381  * Rational(-8, 6).to_r #=> (-4/3)
1382  */
1383 static VALUE
1385 {
1386  return self;
1387 }
1388 
1389 #define id_ceil rb_intern("ceil")
1390 #define f_ceil(x) rb_funcall((x), id_ceil, 0)
1391 
1392 #define id_quo rb_intern("quo")
1393 #define f_quo(x,y) rb_funcall((x), id_quo, 1, (y))
1394 
1395 #define f_reciprocal(x) f_quo(ONE, (x))
1396 
1397 /*
1398  The algorithm here is the method described in CLISP. Bruno Haible has
1399  graciously given permission to use this algorithm. He says, "You can use
1400  it, if you present the following explanation of the algorithm."
1401 
1402  Algorithm (recursively presented):
1403  If x is a rational number, return x.
1404  If x = 0.0, return 0.
1405  If x < 0.0, return (- (rationalize (- x))).
1406  If x > 0.0:
1407  Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
1408  exponent, sign).
1409  If m = 0 or e >= 0: return x = m*2^e.
1410  Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
1411  with smallest possible numerator and denominator.
1412  Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
1413  But in this case the result will be x itself anyway, regardless of
1414  the choice of a. Therefore we can simply ignore this case.
1415  Note 2: At first, we need to consider the closed interval [a,b].
1416  but since a and b have the denominator 2^(|e|+1) whereas x itself
1417  has a denominator <= 2^|e|, we can restrict the search to the open
1418  interval (a,b).
1419  So, for given a and b (0 < a < b) we are searching a rational number
1420  y with a <= y <= b.
1421  Recursive algorithm fraction_between(a,b):
1422  c := (ceiling a)
1423  if c < b
1424  then return c ; because a <= c < b, c integer
1425  else
1426  ; a is not integer (otherwise we would have had c = a < b)
1427  k := c-1 ; k = floor(a), k < a < b <= k+1
1428  return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
1429  ; note 1 <= 1/(b-k) < 1/(a-k)
1430 
1431  You can see that we are actually computing a continued fraction expansion.
1432 
1433  Algorithm (iterative):
1434  If x is rational, return x.
1435  Call (integer-decode-float x). It returns a m,e,s (mantissa,
1436  exponent, sign).
1437  If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
1438  Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
1439  (positive and already in lowest terms because the denominator is a
1440  power of two and the numerator is odd).
1441  Start a continued fraction expansion
1442  p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
1443  Loop
1444  c := (ceiling a)
1445  if c >= b
1446  then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
1447  goto Loop
1448  finally partial_quotient(c).
1449  Here partial_quotient(c) denotes the iteration
1450  i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
1451  At the end, return s * (p[i]/q[i]).
1452  This rational number is already in lowest terms because
1453  p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
1454 */
1455 
1456 static void
1458 {
1459  VALUE c, k, t, p0, p1, p2, q0, q1, q2;
1460 
1461  p0 = ZERO;
1462  p1 = ONE;
1463  q0 = ONE;
1464  q1 = ZERO;
1465 
1466  while (1) {
1467  c = f_ceil(a);
1468  if (f_lt_p(c, b))
1469  break;
1470  k = f_sub(c, ONE);
1471  p2 = f_add(f_mul(k, p1), p0);
1472  q2 = f_add(f_mul(k, q1), q0);
1473  t = f_reciprocal(f_sub(b, k));
1474  b = f_reciprocal(f_sub(a, k));
1475  a = t;
1476  p0 = p1;
1477  q0 = q1;
1478  p1 = p2;
1479  q1 = q2;
1480  }
1481  *p = f_add(f_mul(c, p1), p0);
1482  *q = f_add(f_mul(c, q1), q0);
1483 }
1484 
1485 /*
1486  * call-seq:
1487  * rat.rationalize -> self
1488  * rat.rationalize(eps) -> rational
1489  *
1490  * Returns a simpler approximation of the value if an optional
1491  * argument eps is given (rat-|eps| <= result <= rat+|eps|), self
1492  * otherwise.
1493  *
1494  * For example:
1495  *
1496  * r = Rational(5033165, 16777216)
1497  * r.rationalize #=> (5033165/16777216)
1498  * r.rationalize(Rational('0.01')) #=> (3/10)
1499  * r.rationalize(Rational('0.1')) #=> (1/3)
1500  */
1501 static VALUE
1502 nurat_rationalize(int argc, VALUE *argv, VALUE self)
1503 {
1504  VALUE e, a, b, p, q;
1505 
1506  if (argc == 0)
1507  return self;
1508 
1509  if (f_negative_p(self))
1510  return f_negate(nurat_rationalize(argc, argv, f_abs(self)));
1511 
1512  rb_scan_args(argc, argv, "01", &e);
1513  e = f_abs(e);
1514  a = f_sub(self, e);
1515  b = f_add(self, e);
1516 
1517  if (f_eqeq_p(a, b))
1518  return self;
1519 
1520  nurat_rationalize_internal(a, b, &p, &q);
1521  return f_rational_new2(CLASS_OF(self), p, q);
1522 }
1523 
1524 /* :nodoc: */
1525 static VALUE
1527 {
1528  st_index_t v, h[2];
1529  VALUE n;
1530 
1531  get_dat1(self);
1532  n = rb_hash(dat->num);
1533  h[0] = NUM2LONG(n);
1534  n = rb_hash(dat->den);
1535  h[1] = NUM2LONG(n);
1536  v = rb_memhash(h, sizeof(h));
1537  return LONG2FIX(v);
1538 }
1539 
1540 static VALUE
1542 {
1543  VALUE s;
1544  get_dat1(self);
1545 
1546  s = (*func)(dat->num);
1547  rb_str_cat2(s, "/");
1548  rb_str_concat(s, (*func)(dat->den));
1549 
1550  return s;
1551 }
1552 
1553 /*
1554  * call-seq:
1555  * rat.to_s -> string
1556  *
1557  * Returns the value as a string.
1558  *
1559  * For example:
1560  *
1561  * Rational(2).to_s #=> "2/1"
1562  * Rational(-8, 6).to_s #=> "-4/3"
1563  * Rational('0.5').to_s #=> "1/2"
1564  */
1565 static VALUE
1567 {
1568  return f_format(self, f_to_s);
1569 }
1570 
1571 /*
1572  * call-seq:
1573  * rat.inspect -> string
1574  *
1575  * Returns the value as a string for inspection.
1576  *
1577  * For example:
1578  *
1579  * Rational(2).inspect #=> "(2/1)"
1580  * Rational(-8, 6).inspect #=> "(-4/3)"
1581  * Rational('0.5').inspect #=> "(1/2)"
1582  */
1583 static VALUE
1585 {
1586  VALUE s;
1587 
1588  s = rb_usascii_str_new2("(");
1589  rb_str_concat(s, f_format(self, f_inspect));
1590  rb_str_cat2(s, ")");
1591 
1592  return s;
1593 }
1594 
1595 /* :nodoc: */
1596 static VALUE
1598 {
1599  VALUE a;
1600  get_dat1(self);
1601 
1602  a = rb_assoc_new(dat->num, dat->den);
1603  rb_copy_generic_ivar(a, self);
1604  return a;
1605 }
1606 
1607 /* :nodoc: */
1608 static VALUE
1610 {
1611  get_dat1(self);
1612  Check_Type(a, T_ARRAY);
1613  if (RARRAY_LEN(a) != 2)
1614  rb_raise(rb_eArgError, "marshaled rational must have an array whose length is 2 but %ld", RARRAY_LEN(a));
1615  dat->num = RARRAY_PTR(a)[0];
1616  dat->den = RARRAY_PTR(a)[1];
1617  rb_copy_generic_ivar(self, a);
1618 
1619  if (f_zero_p(dat->den))
1620  rb_raise_zerodiv();
1621 
1622  return self;
1623 }
1624 
1625 /* --- */
1626 
1627 VALUE
1629 {
1630  get_dat1(x);
1631  return f_rational_new_no_reduce2(CLASS_OF(x), dat->den, dat->num);
1632 }
1633 
1634 /*
1635  * call-seq:
1636  * int.gcd(int2) -> integer
1637  *
1638  * Returns the greatest common divisor (always positive). 0.gcd(x)
1639  * and x.gcd(0) return abs(x).
1640  *
1641  * For example:
1642  *
1643  * 2.gcd(2) #=> 2
1644  * 3.gcd(-7) #=> 1
1645  * ((1<<31)-1).gcd((1<<61)-1) #=> 1
1646  */
1647 VALUE
1648 rb_gcd(VALUE self, VALUE other)
1649 {
1650  other = nurat_int_value(other);
1651  return f_gcd(self, other);
1652 }
1653 
1654 /*
1655  * call-seq:
1656  * int.lcm(int2) -> integer
1657  *
1658  * Returns the least common multiple (always positive). 0.lcm(x) and
1659  * x.lcm(0) return zero.
1660  *
1661  * For example:
1662  *
1663  * 2.lcm(2) #=> 2
1664  * 3.lcm(-7) #=> 21
1665  * ((1<<31)-1).lcm((1<<61)-1) #=> 4951760154835678088235319297
1666  */
1667 VALUE
1668 rb_lcm(VALUE self, VALUE other)
1669 {
1670  other = nurat_int_value(other);
1671  return f_lcm(self, other);
1672 }
1673 
1674 /*
1675  * call-seq:
1676  * int.gcdlcm(int2) -> array
1677  *
1678  * Returns an array; [int.gcd(int2), int.lcm(int2)].
1679  *
1680  * For example:
1681  *
1682  * 2.gcdlcm(2) #=> [2, 2]
1683  * 3.gcdlcm(-7) #=> [1, 21]
1684  * ((1<<31)-1).gcdlcm((1<<61)-1) #=> [1, 4951760154835678088235319297]
1685  */
1686 VALUE
1687 rb_gcdlcm(VALUE self, VALUE other)
1688 {
1689  other = nurat_int_value(other);
1690  return rb_assoc_new(f_gcd(self, other), f_lcm(self, other));
1691 }
1692 
1693 VALUE
1695 {
1696  return nurat_s_new_internal(rb_cRational, x, y);
1697 }
1698 
1699 VALUE
1701 {
1703 }
1704 
1705 static VALUE nurat_s_convert(int argc, VALUE *argv, VALUE klass);
1706 
1707 VALUE
1709 {
1710  VALUE a[2];
1711  a[0] = x;
1712  a[1] = y;
1713  return nurat_s_convert(2, a, rb_cRational);
1714 }
1715 
1716 #define id_numerator rb_intern("numerator")
1717 #define f_numerator(x) rb_funcall((x), id_numerator, 0)
1718 
1719 #define id_denominator rb_intern("denominator")
1720 #define f_denominator(x) rb_funcall((x), id_denominator, 0)
1721 
1722 #define id_to_r rb_intern("to_r")
1723 #define f_to_r(x) rb_funcall((x), id_to_r, 0)
1724 
1725 /*
1726  * call-seq:
1727  * num.numerator -> integer
1728  *
1729  * Returns the numerator.
1730  */
1731 static VALUE
1733 {
1734  return f_numerator(f_to_r(self));
1735 }
1736 
1737 /*
1738  * call-seq:
1739  * num.denominator -> integer
1740  *
1741  * Returns the denominator (always positive).
1742  */
1743 static VALUE
1745 {
1746  return f_denominator(f_to_r(self));
1747 }
1748 
1749 /*
1750  * call-seq:
1751  * int.numerator -> self
1752  *
1753  * Returns self.
1754  */
1755 static VALUE
1757 {
1758  return self;
1759 }
1760 
1761 /*
1762  * call-seq:
1763  * int.denominator -> 1
1764  *
1765  * Returns 1.
1766  */
1767 static VALUE
1769 {
1770  return INT2FIX(1);
1771 }
1772 
1773 /*
1774  * call-seq:
1775  * flo.numerator -> integer
1776  *
1777  * Returns the numerator. The result is machine dependent.
1778  *
1779  * For example:
1780  *
1781  * n = 0.3.numerator #=> 5404319552844595
1782  * d = 0.3.denominator #=> 18014398509481984
1783  * n.fdiv(d) #=> 0.3
1784  */
1785 static VALUE
1787 {
1788  double d = RFLOAT_VALUE(self);
1789  if (isinf(d) || isnan(d))
1790  return self;
1791  return rb_call_super(0, 0);
1792 }
1793 
1794 /*
1795  * call-seq:
1796  * flo.denominator -> integer
1797  *
1798  * Returns the denominator (always positive). The result is machine
1799  * dependent.
1800  *
1801  * See numerator.
1802  */
1803 static VALUE
1805 {
1806  double d = RFLOAT_VALUE(self);
1807  if (isinf(d) || isnan(d))
1808  return INT2FIX(1);
1809  return rb_call_super(0, 0);
1810 }
1811 
1812 /*
1813  * call-seq:
1814  * nil.to_r -> (0/1)
1815  *
1816  * Returns zero as a rational.
1817  */
1818 static VALUE
1820 {
1821  return rb_rational_new1(INT2FIX(0));
1822 }
1823 
1824 /*
1825  * call-seq:
1826  * nil.rationalize([eps]) -> (0/1)
1827  *
1828  * Returns zero as a rational. An optional argument eps is always
1829  * ignored.
1830  */
1831 static VALUE
1832 nilclass_rationalize(int argc, VALUE *argv, VALUE self)
1833 {
1834  rb_scan_args(argc, argv, "01", NULL);
1835  return nilclass_to_r(self);
1836 }
1837 
1838 /*
1839  * call-seq:
1840  * int.to_r -> rational
1841  *
1842  * Returns the value as a rational.
1843  *
1844  * For example:
1845  *
1846  * 1.to_r #=> (1/1)
1847  * (1<<64).to_r #=> (18446744073709551616/1)
1848  */
1849 static VALUE
1851 {
1852  return rb_rational_new1(self);
1853 }
1854 
1855 /*
1856  * call-seq:
1857  * int.rationalize([eps]) -> rational
1858  *
1859  * Returns the value as a rational. An optional argument eps is
1860  * always ignored.
1861  */
1862 static VALUE
1863 integer_rationalize(int argc, VALUE *argv, VALUE self)
1864 {
1865  rb_scan_args(argc, argv, "01", NULL);
1866  return integer_to_r(self);
1867 }
1868 
1869 static void
1871 {
1872  double f;
1873  int n;
1874 
1875  f = frexp(RFLOAT_VALUE(self), &n);
1876  f = ldexp(f, DBL_MANT_DIG);
1877  n -= DBL_MANT_DIG;
1878  *rf = rb_dbl2big(f);
1879  *rn = INT2FIX(n);
1880 }
1881 
1882 #if 0
1883 static VALUE
1884 float_decode(VALUE self)
1885 {
1886  VALUE f, n;
1887 
1888  float_decode_internal(self, &f, &n);
1889  return rb_assoc_new(f, n);
1890 }
1891 #endif
1892 
1893 #define id_lshift rb_intern("<<")
1894 #define f_lshift(x,n) rb_funcall((x), id_lshift, 1, (n))
1895 
1896 /*
1897  * call-seq:
1898  * flt.to_r -> rational
1899  *
1900  * Returns the value as a rational.
1901  *
1902  * NOTE: 0.3.to_r isn't the same as '0.3'.to_r. The latter is
1903  * equivalent to '3/10'.to_r, but the former isn't so.
1904  *
1905  * For example:
1906  *
1907  * 2.0.to_r #=> (2/1)
1908  * 2.5.to_r #=> (5/2)
1909  * -0.75.to_r #=> (-3/4)
1910  * 0.0.to_r #=> (0/1)
1911  */
1912 static VALUE
1914 {
1915  VALUE f, n;
1916 
1917  float_decode_internal(self, &f, &n);
1918 #if FLT_RADIX == 2
1919  {
1920  long ln = FIX2LONG(n);
1921 
1922  if (ln == 0)
1923  return f_to_r(f);
1924  if (ln > 0)
1925  return f_to_r(f_lshift(f, n));
1926  ln = -ln;
1927  return rb_rational_new2(f, f_lshift(ONE, INT2FIX(ln)));
1928  }
1929 #else
1930  return f_to_r(f_mul(f, f_expt(INT2FIX(FLT_RADIX), n)));
1931 #endif
1932 }
1933 
1934 /*
1935  * call-seq:
1936  * flt.rationalize([eps]) -> rational
1937  *
1938  * Returns a simpler approximation of the value (flt-|eps| <= result
1939  * <= flt+|eps|). if eps is not given, it will be chosen
1940  * automatically.
1941  *
1942  * For example:
1943  *
1944  * 0.3.rationalize #=> (3/10)
1945  * 1.333.rationalize #=> (1333/1000)
1946  * 1.333.rationalize(0.01) #=> (4/3)
1947  */
1948 static VALUE
1949 float_rationalize(int argc, VALUE *argv, VALUE self)
1950 {
1951  VALUE e, a, b, p, q;
1952 
1953  if (f_negative_p(self))
1954  return f_negate(float_rationalize(argc, argv, f_abs(self)));
1955 
1956  rb_scan_args(argc, argv, "01", &e);
1957 
1958  if (argc != 0) {
1959  e = f_abs(e);
1960  a = f_sub(self, e);
1961  b = f_add(self, e);
1962  }
1963  else {
1964  VALUE f, n;
1965 
1966  float_decode_internal(self, &f, &n);
1967  if (f_zero_p(f) || f_positive_p(n))
1968  return rb_rational_new1(f_lshift(f, n));
1969 
1970 #if FLT_RADIX == 2
1971  a = rb_rational_new2(f_sub(f_mul(TWO, f), ONE),
1972  f_lshift(ONE, f_sub(ONE, n)));
1973  b = rb_rational_new2(f_add(f_mul(TWO, f), ONE),
1974  f_lshift(ONE, f_sub(ONE, n)));
1975 #else
1977  INT2FIX(FLT_RADIX - 1)),
1978  f_expt(INT2FIX(FLT_RADIX), f_sub(ONE, n)));
1980  INT2FIX(FLT_RADIX - 1)),
1981  f_expt(INT2FIX(FLT_RADIX), f_sub(ONE, n)));
1982 #endif
1983  }
1984 
1985  if (f_eqeq_p(a, b))
1986  return f_to_r(self);
1987 
1988  nurat_rationalize_internal(a, b, &p, &q);
1989  return rb_rational_new2(p, q);
1990 }
1991 
1993 
1994 #define WS "\\s*"
1995 #define DIGITS "(?:[0-9](?:_[0-9]|[0-9])*)"
1996 #define NUMERATOR "(?:" DIGITS "?\\.)?" DIGITS "(?:[eE][-+]?" DIGITS ")?"
1997 #define DENOMINATOR DIGITS
1998 #define PATTERN "\\A" WS "([-+])?(" NUMERATOR ")(?:\\/(" DENOMINATOR "))?" WS
1999 
2000 static void
2002 {
2003  static const char rat_pat_source[] = PATTERN;
2004  static const char an_e_pat_source[] = "[eE]";
2005  static const char a_dot_pat_source[] = "\\.";
2006  static const char underscores_pat_source[] = "_+";
2007 
2008  if (rat_pat) return;
2009 
2010  rat_pat = rb_reg_new(rat_pat_source, sizeof rat_pat_source - 1, 0);
2011  rb_gc_register_mark_object(rat_pat);
2012 
2013  an_e_pat = rb_reg_new(an_e_pat_source, sizeof an_e_pat_source - 1, 0);
2014  rb_gc_register_mark_object(an_e_pat);
2015 
2016  a_dot_pat = rb_reg_new(a_dot_pat_source, sizeof a_dot_pat_source - 1, 0);
2017  rb_gc_register_mark_object(a_dot_pat);
2018 
2019  underscores_pat = rb_reg_new(underscores_pat_source,
2020  sizeof underscores_pat_source - 1, 0);
2021  rb_gc_register_mark_object(underscores_pat);
2022 
2023  an_underscore = rb_usascii_str_new2("_");
2024  rb_gc_register_mark_object(an_underscore);
2025 }
2026 
2027 #define id_match rb_intern("match")
2028 #define f_match(x,y) rb_funcall((x), id_match, 1, (y))
2029 
2030 #define id_split rb_intern("split")
2031 #define f_split(x,y) rb_funcall((x), id_split, 1, (y))
2032 
2033 #include <ctype.h>
2034 
2035 static VALUE
2037 {
2038  VALUE s, m;
2039 
2040  s = self;
2041 
2042  if (RSTRING_LEN(s) == 0)
2043  return rb_assoc_new(Qnil, self);
2044 
2045  m = f_match(rat_pat, s);
2046 
2047  if (!NIL_P(m)) {
2048  VALUE v, ifp, exp, ip, fp;
2049  VALUE si = rb_reg_nth_match(1, m);
2050  VALUE nu = rb_reg_nth_match(2, m);
2051  VALUE de = rb_reg_nth_match(3, m);
2052  VALUE re = rb_reg_match_post(m);
2053 
2054  {
2055  VALUE a;
2056 
2057  if (!strpbrk(RSTRING_PTR(nu), "eE")) {
2058  ifp = nu; /* not a copy */
2059  exp = Qnil;
2060  }
2061  else {
2062  a = f_split(nu, an_e_pat);
2063  ifp = RARRAY_PTR(a)[0];
2064  if (RARRAY_LEN(a) != 2)
2065  exp = Qnil;
2066  else
2067  exp = RARRAY_PTR(a)[1];
2068  }
2069 
2070  if (!strchr(RSTRING_PTR(ifp), '.')) {
2071  ip = ifp; /* not a copy */
2072  fp = Qnil;
2073  }
2074  else {
2075  a = f_split(ifp, a_dot_pat);
2076  ip = RARRAY_PTR(a)[0];
2077  if (RARRAY_LEN(a) != 2)
2078  fp = Qnil;
2079  else
2080  fp = RARRAY_PTR(a)[1];
2081  }
2082  }
2083 
2084  v = rb_rational_new1(f_to_i(ip));
2085 
2086  if (!NIL_P(fp)) {
2087  char *p = RSTRING_PTR(fp);
2088  long count = 0;
2089  VALUE l;
2090 
2091  while (*p) {
2092  if (rb_isdigit(*p))
2093  count++;
2094  p++;
2095  }
2096  l = f_expt10(LONG2NUM(count));
2097  v = f_mul(v, l);
2098  v = f_add(v, f_to_i(fp));
2099  v = f_div(v, l);
2100  }
2101  if (!NIL_P(si) && *RSTRING_PTR(si) == '-')
2102  v = f_negate(v);
2103  if (!NIL_P(exp))
2104  v = f_mul(v, f_expt10(f_to_i(exp)));
2105 #if 0
2106  if (!NIL_P(de) && (!NIL_P(fp) || !NIL_P(exp)))
2107  return rb_assoc_new(v, rb_usascii_str_new2("dummy"));
2108 #endif
2109  if (!NIL_P(de))
2110  v = f_div(v, f_to_i(de));
2111 
2112  return rb_assoc_new(v, re);
2113  }
2114  return rb_assoc_new(Qnil, self);
2115 }
2116 
2117 static VALUE
2119 {
2120  VALUE a = string_to_r_internal(self);
2121  if (NIL_P(RARRAY_PTR(a)[0]) || RSTRING_LEN(RARRAY_PTR(a)[1]) > 0) {
2122  VALUE s = f_inspect(self);
2123  rb_raise(rb_eArgError, "invalid value for convert(): %s",
2124  StringValuePtr(s));
2125  }
2126  return RARRAY_PTR(a)[0];
2127 }
2128 
2129 #define id_gsub rb_intern("gsub")
2130 #define f_gsub(x,y,z) rb_funcall((x), id_gsub, 2, (y), (z))
2131 
2132 /*
2133  * call-seq:
2134  * str.to_r -> rational
2135  *
2136  * Returns a rational which denotes the string form. The parser
2137  * ignores leading whitespaces and trailing garbage. Any digit
2138  * sequences can be separated by an underscore. Returns zero for null
2139  * or garbage string.
2140  *
2141  * NOTE: '0.3'.to_r isn't the same as 0.3.to_r. The former is
2142  * equivalent to '3/10'.to_r, but the latter isn't so.
2143  *
2144  * For example:
2145  *
2146  * ' 2 '.to_r #=> (2/1)
2147  * '300/2'.to_r #=> (150/1)
2148  * '-9.2'.to_r #=> (-46/5)
2149  * '-9.2e2'.to_r #=> (-920/1)
2150  * '1_234_567'.to_r #=> (1234567/1)
2151  * '21 june 09'.to_r #=> (21/1)
2152  * '21/06/09'.to_r #=> (7/2)
2153  * 'bwv 1079'.to_r #=> (0/1)
2154  */
2155 static VALUE
2157 {
2158  VALUE s, a, a1, backref;
2159 
2160  backref = rb_backref_get();
2161  rb_match_busy(backref);
2162 
2163  s = f_gsub(self, underscores_pat, an_underscore);
2164  a = string_to_r_internal(s);
2165 
2166  rb_backref_set(backref);
2167 
2168  a1 = RARRAY_PTR(a)[0];
2169  if (!NIL_P(a1)) {
2170  if (TYPE(a1) == T_FLOAT)
2171  rb_raise(rb_eFloatDomainError, "Infinity");
2172  return a1;
2173  }
2174  return rb_rational_new1(INT2FIX(0));
2175 }
2176 
2177 #define id_to_r rb_intern("to_r")
2178 #define f_to_r(x) rb_funcall((x), id_to_r, 0)
2179 
2180 static VALUE
2181 nurat_s_convert(int argc, VALUE *argv, VALUE klass)
2182 {
2183  VALUE a1, a2, backref;
2184 
2185  rb_scan_args(argc, argv, "11", &a1, &a2);
2186 
2187  if (NIL_P(a1) || (argc == 2 && NIL_P(a2)))
2188  rb_raise(rb_eTypeError, "can't convert nil into Rational");
2189 
2190  switch (TYPE(a1)) {
2191  case T_COMPLEX:
2192  if (k_exact_zero_p(RCOMPLEX(a1)->imag))
2193  a1 = RCOMPLEX(a1)->real;
2194  }
2195 
2196  switch (TYPE(a2)) {
2197  case T_COMPLEX:
2198  if (k_exact_zero_p(RCOMPLEX(a2)->imag))
2199  a2 = RCOMPLEX(a2)->real;
2200  }
2201 
2202  backref = rb_backref_get();
2203  rb_match_busy(backref);
2204 
2205  switch (TYPE(a1)) {
2206  case T_FIXNUM:
2207  case T_BIGNUM:
2208  break;
2209  case T_FLOAT:
2210  a1 = f_to_r(a1);
2211  break;
2212  case T_STRING:
2213  a1 = string_to_r_strict(a1);
2214  break;
2215  }
2216 
2217  switch (TYPE(a2)) {
2218  case T_FIXNUM:
2219  case T_BIGNUM:
2220  break;
2221  case T_FLOAT:
2222  a2 = f_to_r(a2);
2223  break;
2224  case T_STRING:
2225  a2 = string_to_r_strict(a2);
2226  break;
2227  }
2228 
2229  rb_backref_set(backref);
2230 
2231  switch (TYPE(a1)) {
2232  case T_RATIONAL:
2233  if (argc == 1 || (k_exact_one_p(a2)))
2234  return a1;
2235  }
2236 
2237  if (argc == 1) {
2238  if (!(k_numeric_p(a1) && k_integer_p(a1)))
2239  return rb_convert_type(a1, T_RATIONAL, "Rational", "to_r");
2240  }
2241  else {
2242  if ((k_numeric_p(a1) && k_numeric_p(a2)) &&
2243  (!f_integer_p(a1) || !f_integer_p(a2)))
2244  return f_div(a1, a2);
2245  }
2246 
2247  {
2248  VALUE argv2[2];
2249  argv2[0] = a1;
2250  argv2[1] = a2;
2251  return nurat_s_new(argc, argv2, klass);
2252  }
2253 }
2254 
2255 /*
2256  * A rational number can be represented as a paired integer number;
2257  * a/b (b>0). Where a is numerator and b is denominator. Integer a
2258  * equals rational a/1 mathematically.
2259  *
2260  * In ruby, you can create rational object with Rational, to_r or
2261  * rationalize method. The return values will be irreducible.
2262  *
2263  * Rational(1) #=> (1/1)
2264  * Rational(2, 3) #=> (2/3)
2265  * Rational(4, -6) #=> (-2/3)
2266  * 3.to_r #=> (3/1)
2267  *
2268  * You can also create rational object from floating-point numbers or
2269  * strings.
2270  *
2271  * Rational(0.3) #=> (5404319552844595/18014398509481984)
2272  * Rational('0.3') #=> (3/10)
2273  * Rational('2/3') #=> (2/3)
2274  *
2275  * 0.3.to_r #=> (5404319552844595/18014398509481984)
2276  * '0.3'.to_r #=> (3/10)
2277  * '2/3'.to_r #=> (2/3)
2278  * 0.3.rationalize #=> (3/10)
2279  *
2280  * A rational object is an exact number, which helps you to write
2281  * program without any rounding errors.
2282  *
2283  * 10.times.inject(0){|t,| t + 0.1} #=> 0.9999999999999999
2284  * 10.times.inject(0){|t,| t + Rational('0.1')} #=> (1/1)
2285  *
2286  * However, when an expression has inexact factor (numerical value or
2287  * operation), will produce an inexact result.
2288  *
2289  * Rational(10) / 3 #=> (10/3)
2290  * Rational(10) / 3.0 #=> 3.3333333333333335
2291  *
2292  * Rational(-8) ** Rational(1, 3)
2293  * #=> (1.0000000000000002+1.7320508075688772i)
2294  */
2295 void
2297 {
2298 #undef rb_intern
2299 #define rb_intern(str) rb_intern_const(str)
2300 
2301  assert(fprintf(stderr, "assert() is now active\n"));
2302 
2303  id_abs = rb_intern("abs");
2304  id_cmp = rb_intern("<=>");
2305  id_convert = rb_intern("convert");
2306  id_eqeq_p = rb_intern("==");
2307  id_expt = rb_intern("**");
2308  id_fdiv = rb_intern("fdiv");
2309  id_floor = rb_intern("floor");
2310  id_idiv = rb_intern("div");
2311  id_inspect = rb_intern("inspect");
2312  id_integer_p = rb_intern("integer?");
2313  id_negate = rb_intern("-@");
2314  id_to_f = rb_intern("to_f");
2315  id_to_i = rb_intern("to_i");
2316  id_to_s = rb_intern("to_s");
2317  id_truncate = rb_intern("truncate");
2318 
2319  rb_cRational = rb_define_class("Rational", rb_cNumeric);
2320 
2322  rb_undef_method(CLASS_OF(rb_cRational), "allocate");
2323 
2324 #if 0
2325  rb_define_private_method(CLASS_OF(rb_cRational), "new!", nurat_s_new_bang, -1);
2327 #else
2329 #endif
2330 
2332 
2333  rb_define_method(rb_cRational, "numerator", nurat_numerator, 0);
2334  rb_define_method(rb_cRational, "denominator", nurat_denominator, 0);
2335 
2343 
2347 
2348 #if 0 /* NUBY */
2349  rb_define_method(rb_cRational, "//", nurat_idiv, 1);
2350 #endif
2351 
2352 #if 0
2353  rb_define_method(rb_cRational, "quot", nurat_quot, 1);
2354  rb_define_method(rb_cRational, "quotrem", nurat_quotrem, 1);
2355 #endif
2356 
2357 #if 0
2358  rb_define_method(rb_cRational, "rational?", nurat_true, 0);
2359  rb_define_method(rb_cRational, "exact?", nurat_true, 0);
2360 #endif
2361 
2366 
2370  rb_define_method(rb_cRational, "rationalize", nurat_rationalize, -1);
2371 
2373 
2376 
2377  rb_define_method(rb_cRational, "marshal_dump", nurat_marshal_dump, 0);
2378  rb_define_method(rb_cRational, "marshal_load", nurat_marshal_load, 1);
2379 
2380  /* --- */
2381 
2382  rb_define_method(rb_cInteger, "gcd", rb_gcd, 1);
2383  rb_define_method(rb_cInteger, "lcm", rb_lcm, 1);
2384  rb_define_method(rb_cInteger, "gcdlcm", rb_gcdlcm, 1);
2385 
2387  rb_define_method(rb_cNumeric, "denominator", numeric_denominator, 0);
2388 
2390  rb_define_method(rb_cInteger, "denominator", integer_denominator, 0);
2391 
2392  rb_define_method(rb_cFloat, "numerator", float_numerator, 0);
2393  rb_define_method(rb_cFloat, "denominator", float_denominator, 0);
2394 
2396  rb_define_method(rb_cNilClass, "rationalize", nilclass_rationalize, -1);
2398  rb_define_method(rb_cInteger, "rationalize", integer_rationalize, -1);
2399  rb_define_method(rb_cFloat, "to_r", float_to_r, 0);
2400  rb_define_method(rb_cFloat, "rationalize", float_rationalize, -1);
2401 
2402  make_patterns();
2403 
2405 
2407 }
2408 
2409 /*
2410 Local variables:
2411 c-file-style: "ruby"
2412 End:
2413 */
2414