Ruby: rational.c Source File

00001 /*
00002   rational.c: Coded by Tadayoshi Funaba 2008-2012
00003 
00004   This implementation is based on Keiju Ishitsuka's Rational library
00005   which is written in ruby.
00006 */
00007 
00008 #include "ruby.h"
00009 #include "internal.h"
00010 #include <math.h>
00011 #include <float.h>
00012 
00013 #ifdef HAVE_IEEEFP_H
00014 #include <ieeefp.h>
00015 #endif
00016 
00017 #define NDEBUG
00018 #include <assert.h>
00019 
00020 #define ZERO INT2FIX(0)
00021 #define ONE INT2FIX(1)
00022 #define TWO INT2FIX(2)
00023 
00024 VALUE rb_cRational;
00025 
00026 static ID id_abs, id_cmp, id_convert, id_eqeq_p, id_expt, id_fdiv,
00027     id_floor, id_idiv, id_inspect, id_integer_p, id_negate, id_to_f,
00028     id_to_i, id_to_s, id_truncate, id_i_num, id_i_den;
00029 
00030 #define f_boolcast(x) ((x) ? Qtrue : Qfalse)
00031 
00032 #define binop(n,op) \
00033 inline static VALUE \
00034 f_##n(VALUE x, VALUE y)\
00035 {\
00036   return rb_funcall(x, (op), 1, y);\
00037 }
00038 
00039 #define fun1(n) \
00040 inline static VALUE \
00041 f_##n(VALUE x)\
00042 {\
00043     return rb_funcall(x, id_##n, 0);\
00044 }
00045 
00046 #define fun2(n) \
00047 inline static VALUE \
00048 f_##n(VALUE x, VALUE y)\
00049 {\
00050     return rb_funcall(x, id_##n, 1, y);\
00051 }
00052 
00053 inline static VALUE
00054 f_add(VALUE x, VALUE y)
00055 {
00056     if (FIXNUM_P(y) && FIX2LONG(y) == 0)
00057         return x;
00058     else if (FIXNUM_P(x) && FIX2LONG(x) == 0)
00059         return y;
00060     return rb_funcall(x, '+', 1, y);
00061 }
00062 
00063 inline static VALUE
00064 f_cmp(VALUE x, VALUE y)
00065 {
00066     if (FIXNUM_P(x) && FIXNUM_P(y)) {
00067         long c = FIX2LONG(x) - FIX2LONG(y);
00068         if (c > 0)
00069             c = 1;
00070         else if (c < 0)
00071             c = -1;
00072         return INT2FIX(c);
00073     }
00074     return rb_funcall(x, id_cmp, 1, y);
00075 }
00076 
00077 inline static VALUE
00078 f_div(VALUE x, VALUE y)
00079 {
00080     if (FIXNUM_P(y) && FIX2LONG(y) == 1)
00081         return x;
00082     return rb_funcall(x, '/', 1, y);
00083 }
00084 
00085 inline static VALUE
00086 f_gt_p(VALUE x, VALUE y)
00087 {
00088     if (FIXNUM_P(x) && FIXNUM_P(y))
00089         return f_boolcast(FIX2LONG(x) > FIX2LONG(y));
00090     return rb_funcall(x, '>', 1, y);
00091 }
00092 
00093 inline static VALUE
00094 f_lt_p(VALUE x, VALUE y)
00095 {
00096     if (FIXNUM_P(x) && FIXNUM_P(y))
00097         return f_boolcast(FIX2LONG(x) < FIX2LONG(y));
00098     return rb_funcall(x, '<', 1, y);
00099 }
00100 
00101 binop(mod, '%')
00102 
00103 inline static VALUE
00104 f_mul(VALUE x, VALUE y)
00105 {
00106     if (FIXNUM_P(y)) {
00107         long iy = FIX2LONG(y);
00108         if (iy == 0) {
00109             if (FIXNUM_P(x) || RB_TYPE_P(x, T_BIGNUM))
00110                 return ZERO;
00111         }
00112         else if (iy == 1)
00113             return x;
00114     }
00115     else if (FIXNUM_P(x)) {
00116         long ix = FIX2LONG(x);
00117         if (ix == 0) {
00118             if (FIXNUM_P(y) || RB_TYPE_P(y, T_BIGNUM))
00119                 return ZERO;
00120         }
00121         else if (ix == 1)
00122             return y;
00123     }
00124     return rb_funcall(x, '*', 1, y);
00125 }
00126 
00127 inline static VALUE
00128 f_sub(VALUE x, VALUE y)
00129 {
00130     if (FIXNUM_P(y) && FIX2LONG(y) == 0)
00131         return x;
00132     return rb_funcall(x, '-', 1, y);
00133 }
00134 
00135 fun1(abs)
00136 fun1(floor)
00137 fun1(inspect)
00138 fun1(integer_p)
00139 fun1(negate)
00140 
00141 inline static VALUE
00142 f_to_i(VALUE x)
00143 {
00144     if (RB_TYPE_P(x, T_STRING))
00145         return rb_str_to_inum(x, 10, 0);
00146     return rb_funcall(x, id_to_i, 0);
00147 }
00148 inline static VALUE
00149 f_to_f(VALUE x)
00150 {
00151     if (RB_TYPE_P(x, T_STRING))
00152         return DBL2NUM(rb_str_to_dbl(x, 0));
00153     return rb_funcall(x, id_to_f, 0);
00154 }
00155 
00156 fun1(to_s)
00157 fun1(truncate)
00158 
00159 inline static VALUE
00160 f_eqeq_p(VALUE x, VALUE y)
00161 {
00162     if (FIXNUM_P(x) && FIXNUM_P(y))
00163         return f_boolcast(FIX2LONG(x) == FIX2LONG(y));
00164     return rb_funcall(x, id_eqeq_p, 1, y);
00165 }
00166 
00167 fun2(expt)
00168 fun2(fdiv)
00169 fun2(idiv)
00170 
00171 #define f_expt10(x) f_expt(INT2FIX(10), x)
00172 
00173 inline static VALUE
00174 f_negative_p(VALUE x)
00175 {
00176     if (FIXNUM_P(x))
00177         return f_boolcast(FIX2LONG(x) < 0);
00178     return rb_funcall(x, '<', 1, ZERO);
00179 }
00180 
00181 #define f_positive_p(x) (!f_negative_p(x))
00182 
00183 inline static VALUE
00184 f_zero_p(VALUE x)
00185 {
00186     switch (TYPE(x)) {
00187       case T_FIXNUM:
00188         return f_boolcast(FIX2LONG(x) == 0);
00189       case T_BIGNUM:
00190         return Qfalse;
00191       case T_RATIONAL:
00192       {
00193           VALUE num = RRATIONAL(x)->num;
00194 
00195           return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 0);
00196       }
00197     }
00198     return rb_funcall(x, id_eqeq_p, 1, ZERO);
00199 }
00200 
00201 #define f_nonzero_p(x) (!f_zero_p(x))
00202 
00203 inline static VALUE
00204 f_one_p(VALUE x)
00205 {
00206     switch (TYPE(x)) {
00207       case T_FIXNUM:
00208         return f_boolcast(FIX2LONG(x) == 1);
00209       case T_BIGNUM:
00210         return Qfalse;
00211       case T_RATIONAL:
00212       {
00213           VALUE num = RRATIONAL(x)->num;
00214           VALUE den = RRATIONAL(x)->den;
00215 
00216           return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 1 &&
00217                             FIXNUM_P(den) && FIX2LONG(den) == 1);
00218       }
00219     }
00220     return rb_funcall(x, id_eqeq_p, 1, ONE);
00221 }
00222 
00223 inline static VALUE
00224 f_minus_one_p(VALUE x)
00225 {
00226     switch (TYPE(x)) {
00227       case T_FIXNUM:
00228         return f_boolcast(FIX2LONG(x) == -1);
00229       case T_BIGNUM:
00230         return Qfalse;
00231       case T_RATIONAL:
00232       {
00233           VALUE num = RRATIONAL(x)->num;
00234           VALUE den = RRATIONAL(x)->den;
00235 
00236           return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == -1 &&
00237                             FIXNUM_P(den) && FIX2LONG(den) == 1);
00238       }
00239     }
00240     return rb_funcall(x, id_eqeq_p, 1, INT2FIX(-1));
00241 }
00242 
00243 inline static VALUE
00244 f_kind_of_p(VALUE x, VALUE c)
00245 {
00246     return rb_obj_is_kind_of(x, c);
00247 }
00248 
00249 inline static VALUE
00250 k_numeric_p(VALUE x)
00251 {
00252     return f_kind_of_p(x, rb_cNumeric);
00253 }
00254 
00255 inline static VALUE
00256 k_integer_p(VALUE x)
00257 {
00258     return f_kind_of_p(x, rb_cInteger);
00259 }
00260 
00261 inline static VALUE
00262 k_float_p(VALUE x)
00263 {
00264     return f_kind_of_p(x, rb_cFloat);
00265 }
00266 
00267 inline static VALUE
00268 k_rational_p(VALUE x)
00269 {
00270     return f_kind_of_p(x, rb_cRational);
00271 }
00272 
00273 #define k_exact_p(x) (!k_float_p(x))
00274 #define k_inexact_p(x) k_float_p(x)
00275 
00276 #define k_exact_zero_p(x) (k_exact_p(x) && f_zero_p(x))
00277 #define k_exact_one_p(x) (k_exact_p(x) && f_one_p(x))
00278 
00279 #ifndef NDEBUG
00280 #define f_gcd f_gcd_orig
00281 #endif
00282 
00283 inline static long
00284 i_gcd(long x, long y)
00285 {
00286     if (x < 0)
00287         x = -x;
00288     if (y < 0)
00289         y = -y;
00290 
00291     if (x == 0)
00292         return y;
00293     if (y == 0)
00294         return x;
00295 
00296     while (x > 0) {
00297         long t = x;
00298         x = y % x;
00299         y = t;
00300     }
00301     return y;
00302 }
00303 
00304 inline static VALUE
00305 f_gcd(VALUE x, VALUE y)
00306 {
00307     VALUE z;
00308 
00309     if (FIXNUM_P(x) && FIXNUM_P(y))
00310         return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
00311 
00312     if (f_negative_p(x))
00313         x = f_negate(x);
00314     if (f_negative_p(y))
00315         y = f_negate(y);
00316 
00317     if (f_zero_p(x))
00318         return y;
00319     if (f_zero_p(y))
00320         return x;
00321 
00322     for (;;) {
00323         if (FIXNUM_P(x)) {
00324             if (FIX2LONG(x) == 0)
00325                 return y;
00326             if (FIXNUM_P(y))
00327                 return LONG2NUM(i_gcd(FIX2LONG(x), FIX2LONG(y)));
00328         }
00329         z = x;
00330         x = f_mod(y, x);
00331         y = z;
00332     }
00333     /* NOTREACHED */
00334 }
00335 
00336 #ifndef NDEBUG
00337 #undef f_gcd
00338 
00339 inline static VALUE
00340 f_gcd(VALUE x, VALUE y)
00341 {
00342     VALUE r = f_gcd_orig(x, y);
00343     if (f_nonzero_p(r)) {
00344         assert(f_zero_p(f_mod(x, r)));
00345         assert(f_zero_p(f_mod(y, r)));
00346     }
00347     return r;
00348 }
00349 #endif
00350 
00351 inline static VALUE
00352 f_lcm(VALUE x, VALUE y)
00353 {
00354     if (f_zero_p(x) || f_zero_p(y))
00355         return ZERO;
00356     return f_abs(f_mul(f_div(x, f_gcd(x, y)), y));
00357 }
00358 
00359 #define get_dat1(x) \
00360     struct RRational *dat;\
00361     dat = ((struct RRational *)(x))
00362 
00363 #define get_dat2(x,y) \
00364     struct RRational *adat, *bdat;\
00365     adat = ((struct RRational *)(x));\
00366     bdat = ((struct RRational *)(y))
00367 
00368 inline static VALUE
00369 nurat_s_new_internal(VALUE klass, VALUE num, VALUE den)
00370 {
00371     NEWOBJ_OF(obj, struct RRational, klass, T_RATIONAL);
00372 
00373     obj->num = num;
00374     obj->den = den;
00375 
00376     return (VALUE)obj;
00377 }
00378 
00379 static VALUE
00380 nurat_s_alloc(VALUE klass)
00381 {
00382     return nurat_s_new_internal(klass, ZERO, ONE);
00383 }
00384 
00385 #define rb_raise_zerodiv() rb_raise(rb_eZeroDivError, "divided by 0")
00386 
00387 #if 0
00388 static VALUE
00389 nurat_s_new_bang(int argc, VALUE *argv, VALUE klass)
00390 {
00391     VALUE num, den;
00392 
00393     switch (rb_scan_args(argc, argv, "11", &num, &den)) {
00394       case 1:
00395         if (!k_integer_p(num))
00396             num = f_to_i(num);
00397         den = ONE;
00398         break;
00399       default:
00400         if (!k_integer_p(num))
00401             num = f_to_i(num);
00402         if (!k_integer_p(den))
00403             den = f_to_i(den);
00404 
00405         switch (FIX2INT(f_cmp(den, ZERO))) {
00406           case -1:
00407             num = f_negate(num);
00408             den = f_negate(den);
00409             break;
00410           case 0:
00411             rb_raise_zerodiv();
00412             break;
00413         }
00414         break;
00415     }
00416 
00417     return nurat_s_new_internal(klass, num, den);
00418 }
00419 #endif
00420 
00421 inline static VALUE
00422 f_rational_new_bang1(VALUE klass, VALUE x)
00423 {
00424     return nurat_s_new_internal(klass, x, ONE);
00425 }
00426 
00427 inline static VALUE
00428 f_rational_new_bang2(VALUE klass, VALUE x, VALUE y)
00429 {
00430     assert(f_positive_p(y));
00431     assert(f_nonzero_p(y));
00432     return nurat_s_new_internal(klass, x, y);
00433 }
00434 
00435 #ifdef CANONICALIZATION_FOR_MATHN
00436 #define CANON
00437 #endif
00438 
00439 #ifdef CANON
00440 static int canonicalization = 0;
00441 
00442 RUBY_FUNC_EXPORTED void
00443 nurat_canonicalization(int f)
00444 {
00445     canonicalization = f;
00446 }
00447 #endif
00448 
00449 inline static void
00450 nurat_int_check(VALUE num)
00451 {
00452     switch (TYPE(num)) {
00453       case T_FIXNUM:
00454       case T_BIGNUM:
00455         break;
00456       default:
00457         if (!k_numeric_p(num) || !f_integer_p(num))
00458             rb_raise(rb_eTypeError, "not an integer");
00459     }
00460 }
00461 
00462 inline static VALUE
00463 nurat_int_value(VALUE num)
00464 {
00465     nurat_int_check(num);
00466     if (!k_integer_p(num))
00467         num = f_to_i(num);
00468     return num;
00469 }
00470 
00471 inline static VALUE
00472 nurat_s_canonicalize_internal(VALUE klass, VALUE num, VALUE den)
00473 {
00474     VALUE gcd;
00475 
00476     switch (FIX2INT(f_cmp(den, ZERO))) {
00477       case -1:
00478         num = f_negate(num);
00479         den = f_negate(den);
00480         break;
00481       case 0:
00482         rb_raise_zerodiv();
00483         break;
00484     }
00485 
00486     gcd = f_gcd(num, den);
00487     num = f_idiv(num, gcd);
00488     den = f_idiv(den, gcd);
00489 
00490 #ifdef CANON
00491     if (f_one_p(den) && canonicalization)
00492         return num;
00493 #endif
00494     return nurat_s_new_internal(klass, num, den);
00495 }
00496 
00497 inline static VALUE
00498 nurat_s_canonicalize_internal_no_reduce(VALUE klass, VALUE num, VALUE den)
00499 {
00500     switch (FIX2INT(f_cmp(den, ZERO))) {
00501       case -1:
00502         num = f_negate(num);
00503         den = f_negate(den);
00504         break;
00505       case 0:
00506         rb_raise_zerodiv();
00507         break;
00508     }
00509 
00510 #ifdef CANON
00511     if (f_one_p(den) && canonicalization)
00512         return num;
00513 #endif
00514     return nurat_s_new_internal(klass, num, den);
00515 }
00516 
00517 static VALUE
00518 nurat_s_new(int argc, VALUE *argv, VALUE klass)
00519 {
00520     VALUE num, den;
00521 
00522     switch (rb_scan_args(argc, argv, "11", &num, &den)) {
00523       case 1:
00524         num = nurat_int_value(num);
00525         den = ONE;
00526         break;
00527       default:
00528         num = nurat_int_value(num);
00529         den = nurat_int_value(den);
00530         break;
00531     }
00532 
00533     return nurat_s_canonicalize_internal(klass, num, den);
00534 }
00535 
00536 inline static VALUE
00537 f_rational_new1(VALUE klass, VALUE x)
00538 {
00539     assert(!k_rational_p(x));
00540     return nurat_s_canonicalize_internal(klass, x, ONE);
00541 }
00542 
00543 inline static VALUE
00544 f_rational_new2(VALUE klass, VALUE x, VALUE y)
00545 {
00546     assert(!k_rational_p(x));
00547     assert(!k_rational_p(y));
00548     return nurat_s_canonicalize_internal(klass, x, y);
00549 }
00550 
00551 inline static VALUE
00552 f_rational_new_no_reduce1(VALUE klass, VALUE x)
00553 {
00554     assert(!k_rational_p(x));
00555     return nurat_s_canonicalize_internal_no_reduce(klass, x, ONE);
00556 }
00557 
00558 inline static VALUE
00559 f_rational_new_no_reduce2(VALUE klass, VALUE x, VALUE y)
00560 {
00561     assert(!k_rational_p(x));
00562     assert(!k_rational_p(y));
00563     return nurat_s_canonicalize_internal_no_reduce(klass, x, y);
00564 }
00565 
00566 /*
00567  * call-seq:
00568  *    Rational(x[, y])  ->  numeric
00569  *
00570  * Returns x/y;
00571  *
00572  *    Rational(1, 2)   #=> (1/2)
00573  *    Rational('1/2')  #=> (1/2)
00574  *    Rational(nil)    #=> TypeError
00575  *    Rational(1, nil) #=> TypeError
00576  *
00577  * Syntax of string form:
00578  *
00579  *   string form = extra spaces , rational , extra spaces ;
00580  *   rational = [ sign ] , unsigned rational ;
00581  *   unsigned rational = numerator | numerator , "/" , denominator ;
00582  *   numerator = integer part | fractional part | integer part , fractional part ;
00583  *   denominator = digits ;
00584  *   integer part = digits ;
00585  *   fractional part = "." , digits , [ ( "e" | "E" ) , [ sign ] , digits ] ;
00586  *   sign = "-" | "+" ;
00587  *   digits = digit , { digit | "_" , digit } ;
00588  *   digit = "0" | "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9" ;
00589  *   extra spaces = ? \s* ? ;
00590  *
00591  * See String#to_r.
00592  */
00593 static VALUE
00594 nurat_f_rational(int argc, VALUE *argv, VALUE klass)
00595 {
00596     return rb_funcall2(rb_cRational, id_convert, argc, argv);
00597 }
00598 
00599 /*
00600  * call-seq:
00601  *    rat.numerator  ->  integer
00602  *
00603  * Returns the numerator.
00604  *
00605  *    Rational(7).numerator        #=> 7
00606  *    Rational(7, 1).numerator     #=> 7
00607  *    Rational(9, -4).numerator    #=> -9
00608  *    Rational(-2, -10).numerator  #=> 1
00609  */
00610 static VALUE
00611 nurat_numerator(VALUE self)
00612 {
00613     get_dat1(self);
00614     return dat->num;
00615 }
00616 
00617 /*
00618  * call-seq:
00619  *    rat.denominator  ->  integer
00620  *
00621  * Returns the denominator (always positive).
00622  *
00623  *    Rational(7).denominator             #=> 1
00624  *    Rational(7, 1).denominator          #=> 1
00625  *    Rational(9, -4).denominator         #=> 4
00626  *    Rational(-2, -10).denominator       #=> 5
00627  *    rat.numerator.gcd(rat.denominator)  #=> 1
00628  */
00629 static VALUE
00630 nurat_denominator(VALUE self)
00631 {
00632     get_dat1(self);
00633     return dat->den;
00634 }
00635 
00636 #ifndef NDEBUG
00637 #define f_imul f_imul_orig
00638 #endif
00639 
00640 inline static VALUE
00641 f_imul(long a, long b)
00642 {
00643     VALUE r;
00644 
00645     if (a == 0 || b == 0)
00646         return ZERO;
00647     else if (a == 1)
00648         return LONG2NUM(b);
00649     else if (b == 1)
00650         return LONG2NUM(a);
00651 
00652     if (MUL_OVERFLOW_LONG_P(a, b))
00653         r = rb_big_mul(rb_int2big(a), rb_int2big(b));
00654     else
00655         r = LONG2NUM(a * b);
00656     return r;
00657 }
00658 
00659 #ifndef NDEBUG
00660 #undef f_imul
00661 
00662 inline static VALUE
00663 f_imul(long x, long y)
00664 {
00665     VALUE r = f_imul_orig(x, y);
00666     assert(f_eqeq_p(r, f_mul(LONG2NUM(x), LONG2NUM(y))));
00667     return r;
00668 }
00669 #endif
00670 
00671 inline static VALUE
00672 f_addsub(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
00673 {
00674     VALUE num, den;
00675 
00676     if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
00677         FIXNUM_P(bnum) && FIXNUM_P(bden)) {
00678         long an = FIX2LONG(anum);
00679         long ad = FIX2LONG(aden);
00680         long bn = FIX2LONG(bnum);
00681         long bd = FIX2LONG(bden);
00682         long ig = i_gcd(ad, bd);
00683 
00684         VALUE g = LONG2NUM(ig);
00685         VALUE a = f_imul(an, bd / ig);
00686         VALUE b = f_imul(bn, ad / ig);
00687         VALUE c;
00688 
00689         if (k == '+')
00690             c = f_add(a, b);
00691         else
00692             c = f_sub(a, b);
00693 
00694         b = f_idiv(aden, g);
00695         g = f_gcd(c, g);
00696         num = f_idiv(c, g);
00697         a = f_idiv(bden, g);
00698         den = f_mul(a, b);
00699     }
00700     else {
00701         VALUE g = f_gcd(aden, bden);
00702         VALUE a = f_mul(anum, f_idiv(bden, g));
00703         VALUE b = f_mul(bnum, f_idiv(aden, g));
00704         VALUE c;
00705 
00706         if (k == '+')
00707             c = f_add(a, b);
00708         else
00709             c = f_sub(a, b);
00710 
00711         b = f_idiv(aden, g);
00712         g = f_gcd(c, g);
00713         num = f_idiv(c, g);
00714         a = f_idiv(bden, g);
00715         den = f_mul(a, b);
00716     }
00717     return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
00718 }
00719 
00720 /*
00721  * call-seq:
00722  *    rat + numeric  ->  numeric
00723  *
00724  * Performs addition.
00725  *
00726  *    Rational(2, 3)  + Rational(2, 3)   #=> (4/3)
00727  *    Rational(900)   + Rational(1)      #=> (900/1)
00728  *    Rational(-2, 9) + Rational(-9, 2)  #=> (-85/18)
00729  *    Rational(9, 8)  + 4                #=> (41/8)
00730  *    Rational(20, 9) + 9.8              #=> 12.022222222222222
00731  */
00732 static VALUE
00733 nurat_add(VALUE self, VALUE other)
00734 {
00735     switch (TYPE(other)) {
00736       case T_FIXNUM:
00737       case T_BIGNUM:
00738         {
00739             get_dat1(self);
00740 
00741             return f_addsub(self,
00742                             dat->num, dat->den,
00743                             other, ONE, '+');
00744         }
00745       case T_FLOAT:
00746         return f_add(f_to_f(self), other);
00747       case T_RATIONAL:
00748         {
00749             get_dat2(self, other);
00750 
00751             return f_addsub(self,
00752                             adat->num, adat->den,
00753                             bdat->num, bdat->den, '+');
00754         }
00755       default:
00756         return rb_num_coerce_bin(self, other, '+');
00757     }
00758 }
00759 
00760 /*
00761  * call-seq:
00762  *    rat - numeric  ->  numeric
00763  *
00764  * Performs subtraction.
00765  *
00766  *    Rational(2, 3)  - Rational(2, 3)   #=> (0/1)
00767  *    Rational(900)   - Rational(1)      #=> (899/1)
00768  *    Rational(-2, 9) - Rational(-9, 2)  #=> (77/18)
00769  *    Rational(9, 8)  - 4                #=> (23/8)
00770  *    Rational(20, 9) - 9.8              #=> -7.577777777777778
00771  */
00772 static VALUE
00773 nurat_sub(VALUE self, VALUE other)
00774 {
00775     switch (TYPE(other)) {
00776       case T_FIXNUM:
00777       case T_BIGNUM:
00778         {
00779             get_dat1(self);
00780 
00781             return f_addsub(self,
00782                             dat->num, dat->den,
00783                             other, ONE, '-');
00784         }
00785       case T_FLOAT:
00786         return f_sub(f_to_f(self), other);
00787       case T_RATIONAL:
00788         {
00789             get_dat2(self, other);
00790 
00791             return f_addsub(self,
00792                             adat->num, adat->den,
00793                             bdat->num, bdat->den, '-');
00794         }
00795       default:
00796         return rb_num_coerce_bin(self, other, '-');
00797     }
00798 }
00799 
00800 inline static VALUE
00801 f_muldiv(VALUE self, VALUE anum, VALUE aden, VALUE bnum, VALUE bden, int k)
00802 {
00803     VALUE num, den;
00804 
00805     if (k == '/') {
00806         VALUE t;
00807 
00808         if (f_negative_p(bnum)) {
00809             anum = f_negate(anum);
00810             bnum = f_negate(bnum);
00811         }
00812         t = bnum;
00813         bnum = bden;
00814         bden = t;
00815     }
00816 
00817     if (FIXNUM_P(anum) && FIXNUM_P(aden) &&
00818         FIXNUM_P(bnum) && FIXNUM_P(bden)) {
00819         long an = FIX2LONG(anum);
00820         long ad = FIX2LONG(aden);
00821         long bn = FIX2LONG(bnum);
00822         long bd = FIX2LONG(bden);
00823         long g1 = i_gcd(an, bd);
00824         long g2 = i_gcd(ad, bn);
00825 
00826         num = f_imul(an / g1, bn / g2);
00827         den = f_imul(ad / g2, bd / g1);
00828     }
00829     else {
00830         VALUE g1 = f_gcd(anum, bden);
00831         VALUE g2 = f_gcd(aden, bnum);
00832 
00833         num = f_mul(f_idiv(anum, g1), f_idiv(bnum, g2));
00834         den = f_mul(f_idiv(aden, g2), f_idiv(bden, g1));
00835     }
00836     return f_rational_new_no_reduce2(CLASS_OF(self), num, den);
00837 }
00838 
00839 /*
00840  * call-seq:
00841  *    rat * numeric  ->  numeric
00842  *
00843  * Performs multiplication.
00844  *
00845  *    Rational(2, 3)  * Rational(2, 3)   #=> (4/9)
00846  *    Rational(900)   * Rational(1)      #=> (900/1)
00847  *    Rational(-2, 9) * Rational(-9, 2)  #=> (1/1)
00848  *    Rational(9, 8)  * 4                #=> (9/2)
00849  *    Rational(20, 9) * 9.8              #=> 21.77777777777778
00850  */
00851 static VALUE
00852 nurat_mul(VALUE self, VALUE other)
00853 {
00854     switch (TYPE(other)) {
00855       case T_FIXNUM:
00856       case T_BIGNUM:
00857         {
00858             get_dat1(self);
00859 
00860             return f_muldiv(self,
00861                             dat->num, dat->den,
00862                             other, ONE, '*');
00863         }
00864       case T_FLOAT:
00865         return f_mul(f_to_f(self), other);
00866       case T_RATIONAL:
00867         {
00868             get_dat2(self, other);
00869 
00870             return f_muldiv(self,
00871                             adat->num, adat->den,
00872                             bdat->num, bdat->den, '*');
00873         }
00874       default:
00875         return rb_num_coerce_bin(self, other, '*');
00876     }
00877 }
00878 
00879 /*
00880  * call-seq:
00881  *    rat / numeric     ->  numeric
00882  *    rat.quo(numeric)  ->  numeric
00883  *
00884  * Performs division.
00885  *
00886  *    Rational(2, 3)  / Rational(2, 3)   #=> (1/1)
00887  *    Rational(900)   / Rational(1)      #=> (900/1)
00888  *    Rational(-2, 9) / Rational(-9, 2)  #=> (4/81)
00889  *    Rational(9, 8)  / 4                #=> (9/32)
00890  *    Rational(20, 9) / 9.8              #=> 0.22675736961451246
00891  */
00892 static VALUE
00893 nurat_div(VALUE self, VALUE other)
00894 {
00895     switch (TYPE(other)) {
00896       case T_FIXNUM:
00897       case T_BIGNUM:
00898         if (f_zero_p(other))
00899             rb_raise_zerodiv();
00900         {
00901             get_dat1(self);
00902 
00903             return f_muldiv(self,
00904                             dat->num, dat->den,
00905                             other, ONE, '/');
00906         }
00907       case T_FLOAT:
00908         {
00909             double x = RFLOAT_VALUE(other), den;
00910             get_dat1(self);
00911 
00912             if (isnan(x)) return DBL2NUM(NAN);
00913             if (isinf(x)) return INT2FIX(0);
00914             if (x != 0.0 && modf(x, &den) == 0.0) {
00915                 return rb_rational_raw2(dat->num, f_mul(rb_dbl2big(den), dat->den));
00916             }
00917         }
00918         return rb_funcall(f_to_f(self), '/', 1, other);
00919       case T_RATIONAL:
00920         if (f_zero_p(other))
00921             rb_raise_zerodiv();
00922         {
00923             get_dat2(self, other);
00924 
00925             if (f_one_p(self))
00926                 return f_rational_new_no_reduce2(CLASS_OF(self),
00927                                                  bdat->den, bdat->num);
00928 
00929             return f_muldiv(self,
00930                             adat->num, adat->den,
00931                             bdat->num, bdat->den, '/');
00932         }
00933       default:
00934         return rb_num_coerce_bin(self, other, '/');
00935     }
00936 }
00937 
00938 /*
00939  * call-seq:
00940  *    rat.fdiv(numeric)  ->  float
00941  *
00942  * Performs division and returns the value as a float.
00943  *
00944  *    Rational(2, 3).fdiv(1)       #=> 0.6666666666666666
00945  *    Rational(2, 3).fdiv(0.5)     #=> 1.3333333333333333
00946  *    Rational(2).fdiv(3)          #=> 0.6666666666666666
00947  */
00948 static VALUE
00949 nurat_fdiv(VALUE self, VALUE other)
00950 {
00951     if (f_zero_p(other))
00952         return f_div(self, f_to_f(other));
00953     return f_to_f(f_div(self, other));
00954 }
00955 
00956 inline static VALUE
00957 f_odd_p(VALUE integer)
00958 {
00959     if (rb_funcall(integer, '%', 1, INT2FIX(2)) != INT2FIX(0)) {
00960         return Qtrue;
00961     }
00962     return Qfalse;
00963 
00964 }
00965 
00966 /*
00967  * call-seq:
00968  *    rat ** numeric  ->  numeric
00969  *
00970  * Performs exponentiation.
00971  *
00972  *    Rational(2)    ** Rational(3)    #=> (8/1)
00973  *    Rational(10)   ** -2             #=> (1/100)
00974  *    Rational(10)   ** -2.0           #=> 0.01
00975  *    Rational(-4)   ** Rational(1,2)  #=> (1.2246063538223773e-16+2.0i)
00976  *    Rational(1, 2) ** 0              #=> (1/1)
00977  *    Rational(1, 2) ** 0.0            #=> 1.0
00978  */
00979 static VALUE
00980 nurat_expt(VALUE self, VALUE other)
00981 {
00982     if (k_numeric_p(other) && k_exact_zero_p(other))
00983         return f_rational_new_bang1(CLASS_OF(self), ONE);
00984 
00985     if (k_rational_p(other)) {
00986         get_dat1(other);
00987 
00988         if (f_one_p(dat->den))
00989             other = dat->num; /* c14n */
00990     }
00991 
00992     /* Deal with special cases of 0**n and 1**n */
00993     if (k_numeric_p(other) && k_exact_p(other)) {
00994         get_dat1(self);
00995         if (f_one_p(dat->den))
00996             if (f_one_p(dat->num))
00997                 return f_rational_new_bang1(CLASS_OF(self), ONE);
00998             else if (f_minus_one_p(dat->num) && k_integer_p(other))
00999                 return f_rational_new_bang1(CLASS_OF(self), INT2FIX(f_odd_p(other) ? -1 : 1));
01000             else if (f_zero_p(dat->num))
01001                 if (FIX2INT(f_cmp(other, ZERO)) == -1)
01002                     rb_raise_zerodiv();
01003                 else
01004                     return f_rational_new_bang1(CLASS_OF(self), ZERO);
01005     }
01006 
01007     /* General case */
01008     switch (TYPE(other)) {
01009       case T_FIXNUM:
01010         {
01011             VALUE num, den;
01012 
01013             get_dat1(self);
01014 
01015             switch (FIX2INT(f_cmp(other, ZERO))) {
01016               case 1:
01017                 num = f_expt(dat->num, other);
01018                 den = f_expt(dat->den, other);
01019                 break;
01020               case -1:
01021                 num = f_expt(dat->den, f_negate(other));
01022                 den = f_expt(dat->num, f_negate(other));
01023                 break;
01024               default:
01025                 num = ONE;
01026                 den = ONE;
01027                 break;
01028             }
01029             return f_rational_new2(CLASS_OF(self), num, den);
01030         }
01031       case T_BIGNUM:
01032         rb_warn("in a**b, b may be too big");
01033         /* fall through */
01034       case T_FLOAT:
01035       case T_RATIONAL:
01036         return f_expt(f_to_f(self), other);
01037       default:
01038         return rb_num_coerce_bin(self, other, id_expt);
01039     }
01040 }
01041 
01042 /*
01043  * call-seq:
01044  *    rational <=> numeric  ->  -1, 0, +1 or nil
01045  *
01046  * Performs comparison and returns -1, 0, or +1.
01047  *
01048  * +nil+ is returned if the two values are incomparable.
01049  *
01050  *    Rational(2, 3)  <=> Rational(2, 3)  #=> 0
01051  *    Rational(5)     <=> 5               #=> 0
01052  *    Rational(2,3)   <=> Rational(1,3)   #=> 1
01053  *    Rational(1,3)   <=> 1               #=> -1
01054  *    Rational(1,3)   <=> 0.3             #=> 1
01055  */
01056 static VALUE
01057 nurat_cmp(VALUE self, VALUE other)
01058 {
01059     switch (TYPE(other)) {
01060       case T_FIXNUM:
01061       case T_BIGNUM:
01062         {
01063             get_dat1(self);
01064 
01065             if (FIXNUM_P(dat->den) && FIX2LONG(dat->den) == 1)
01066                 return f_cmp(dat->num, other); /* c14n */
01067             return f_cmp(self, f_rational_new_bang1(CLASS_OF(self), other));
01068         }
01069       case T_FLOAT:
01070         return f_cmp(f_to_f(self), other);
01071       case T_RATIONAL:
01072         {
01073             VALUE num1, num2;
01074 
01075             get_dat2(self, other);
01076 
01077             if (FIXNUM_P(adat->num) && FIXNUM_P(adat->den) &&
01078                 FIXNUM_P(bdat->num) && FIXNUM_P(bdat->den)) {
01079                 num1 = f_imul(FIX2LONG(adat->num), FIX2LONG(bdat->den));
01080                 num2 = f_imul(FIX2LONG(bdat->num), FIX2LONG(adat->den));
01081             }
01082             else {
01083                 num1 = f_mul(adat->num, bdat->den);
01084                 num2 = f_mul(bdat->num, adat->den);
01085             }
01086             return f_cmp(f_sub(num1, num2), ZERO);
01087         }
01088       default:
01089         return rb_num_coerce_cmp(self, other, id_cmp);
01090     }
01091 }
01092 
01093 /*
01094  * call-seq:
01095  *    rat == object  ->  true or false
01096  *
01097  * Returns true if rat equals object numerically.
01098  *
01099  *    Rational(2, 3)  == Rational(2, 3)   #=> true
01100  *    Rational(5)     == 5                #=> true
01101  *    Rational(0)     == 0.0              #=> true
01102  *    Rational('1/3') == 0.33             #=> false
01103  *    Rational('1/2') == '1/2'            #=> false
01104  */
01105 static VALUE
01106 nurat_eqeq_p(VALUE self, VALUE other)
01107 {
01108     switch (TYPE(other)) {
01109       case T_FIXNUM:
01110       case T_BIGNUM:
01111         {
01112             get_dat1(self);
01113 
01114             if (f_zero_p(dat->num) && f_zero_p(other))
01115                 return Qtrue;
01116 
01117             if (!FIXNUM_P(dat->den))
01118                 return Qfalse;
01119             if (FIX2LONG(dat->den) != 1)
01120                 return Qfalse;
01121             if (f_eqeq_p(dat->num, other))
01122                 return Qtrue;
01123             return Qfalse;
01124         }
01125       case T_FLOAT:
01126         return f_eqeq_p(f_to_f(self), other);
01127       case T_RATIONAL:
01128         {
01129             get_dat2(self, other);
01130 
01131             if (f_zero_p(adat->num) && f_zero_p(bdat->num))
01132                 return Qtrue;
01133 
01134             return f_boolcast(f_eqeq_p(adat->num, bdat->num) &&
01135                               f_eqeq_p(adat->den, bdat->den));
01136         }
01137       default:
01138         return f_eqeq_p(other, self);
01139     }
01140 }
01141 
01142 /* :nodoc: */
01143 static VALUE
01144 nurat_coerce(VALUE self, VALUE other)
01145 {
01146     switch (TYPE(other)) {
01147       case T_FIXNUM:
01148       case T_BIGNUM:
01149         return rb_assoc_new(f_rational_new_bang1(CLASS_OF(self), other), self);
01150       case T_FLOAT:
01151         return rb_assoc_new(other, f_to_f(self));
01152       case T_RATIONAL:
01153         return rb_assoc_new(other, self);
01154       case T_COMPLEX:
01155         if (k_exact_zero_p(RCOMPLEX(other)->imag))
01156             return rb_assoc_new(f_rational_new_bang1
01157                                 (CLASS_OF(self), RCOMPLEX(other)->real), self);
01158         else
01159             return rb_assoc_new(other, rb_Complex(self, INT2FIX(0)));
01160     }
01161 
01162     rb_raise(rb_eTypeError, "%s can't be coerced into %s",
01163              rb_obj_classname(other), rb_obj_classname(self));
01164     return Qnil;
01165 }
01166 
01167 #if 0
01168 /* :nodoc: */
01169 static VALUE
01170 nurat_idiv(VALUE self, VALUE other)
01171 {
01172     return f_idiv(self, other);
01173 }
01174 
01175 /* :nodoc: */
01176 static VALUE
01177 nurat_quot(VALUE self, VALUE other)
01178 {
01179     return f_truncate(f_div(self, other));
01180 }
01181 
01182 /* :nodoc: */
01183 static VALUE
01184 nurat_quotrem(VALUE self, VALUE other)
01185 {
01186     VALUE val = f_truncate(f_div(self, other));
01187     return rb_assoc_new(val, f_sub(self, f_mul(other, val)));
01188 }
01189 #endif
01190 
01191 #if 0
01192 /* :nodoc: */
01193 static VALUE
01194 nurat_true(VALUE self)
01195 {
01196     return Qtrue;
01197 }
01198 #endif
01199 
01200 static VALUE
01201 nurat_floor(VALUE self)
01202 {
01203     get_dat1(self);
01204     return f_idiv(dat->num, dat->den);
01205 }
01206 
01207 static VALUE
01208 nurat_ceil(VALUE self)
01209 {
01210     get_dat1(self);
01211     return f_negate(f_idiv(f_negate(dat->num), dat->den));
01212 }
01213 
01214 /*
01215  * call-seq:
01216  *    rat.to_i  ->  integer
01217  *
01218  * Returns the truncated value as an integer.
01219  *
01220  * Equivalent to
01221  *    rat.truncate.
01222  *
01223  *    Rational(2, 3).to_i   #=> 0
01224  *    Rational(3).to_i      #=> 3
01225  *    Rational(300.6).to_i  #=> 300
01226  *    Rational(98,71).to_i  #=> 1
01227  *    Rational(-30,2).to_i  #=> -15
01228  */
01229 static VALUE
01230 nurat_truncate(VALUE self)
01231 {
01232     get_dat1(self);
01233     if (f_negative_p(dat->num))
01234         return f_negate(f_idiv(f_negate(dat->num), dat->den));
01235     return f_idiv(dat->num, dat->den);
01236 }
01237 
01238 static VALUE
01239 nurat_round(VALUE self)
01240 {
01241     VALUE num, den, neg;
01242 
01243     get_dat1(self);
01244 
01245     num = dat->num;
01246     den = dat->den;
01247     neg = f_negative_p(num);
01248 
01249     if (neg)
01250         num = f_negate(num);
01251 
01252     num = f_add(f_mul(num, TWO), den);
01253     den = f_mul(den, TWO);
01254     num = f_idiv(num, den);
01255 
01256     if (neg)
01257         num = f_negate(num);
01258 
01259     return num;
01260 }
01261 
01262 static VALUE
01263 f_round_common(int argc, VALUE *argv, VALUE self, VALUE (*func)(VALUE))
01264 {
01265     VALUE n, b, s;
01266 
01267     if (argc == 0)
01268         return (*func)(self);
01269 
01270     rb_scan_args(argc, argv, "01", &n);
01271 
01272     if (!k_integer_p(n))
01273         rb_raise(rb_eTypeError, "not an integer");
01274 
01275     b = f_expt10(n);
01276     s = f_mul(self, b);
01277 
01278     if (k_float_p(s)) {
01279         if (f_lt_p(n, ZERO))
01280             return ZERO;
01281         return self;
01282     }
01283 
01284     if (!k_rational_p(s)) {
01285         s = f_rational_new_bang1(CLASS_OF(self), s);
01286     }
01287 
01288     s = (*func)(s);
01289 
01290     s = f_div(f_rational_new_bang1(CLASS_OF(self), s), b);
01291 
01292     if (f_lt_p(n, ONE))
01293         s = f_to_i(s);
01294 
01295     return s;
01296 }
01297 
01298 /*
01299  * call-seq:
01300  *    rat.floor               ->  integer
01301  *    rat.floor(precision=0)  ->  rational
01302  *
01303  * Returns the truncated value (toward negative infinity).
01304  *
01305  *    Rational(3).floor      #=> 3
01306  *    Rational(2, 3).floor   #=> 0
01307  *    Rational(-3, 2).floor  #=> -1
01308  *
01309  *           decimal      -  1  2  3 . 4  5  6
01310  *                          ^  ^  ^  ^   ^  ^
01311  *          precision      -3 -2 -1  0  +1 +2
01312  *
01313  *    '%f' % Rational('-123.456').floor(+1)  #=> "-123.500000"
01314  *    '%f' % Rational('-123.456').floor(-1)  #=> "-130.000000"
01315  */
01316 static VALUE
01317 nurat_floor_n(int argc, VALUE *argv, VALUE self)
01318 {
01319     return f_round_common(argc, argv, self, nurat_floor);
01320 }
01321 
01322 /*
01323  * call-seq:
01324  *    rat.ceil               ->  integer
01325  *    rat.ceil(precision=0)  ->  rational
01326  *
01327  * Returns the truncated value (toward positive infinity).
01328  *
01329  *    Rational(3).ceil      #=> 3
01330  *    Rational(2, 3).ceil   #=> 1
01331  *    Rational(-3, 2).ceil  #=> -1
01332  *
01333  *           decimal      -  1  2  3 . 4  5  6
01334  *                          ^  ^  ^  ^   ^  ^
01335  *          precision      -3 -2 -1  0  +1 +2
01336  *
01337  *    '%f' % Rational('-123.456').ceil(+1)  #=> "-123.400000"
01338  *    '%f' % Rational('-123.456').ceil(-1)  #=> "-120.000000"
01339  */
01340 static VALUE
01341 nurat_ceil_n(int argc, VALUE *argv, VALUE self)
01342 {
01343     return f_round_common(argc, argv, self, nurat_ceil);
01344 }
01345 
01346 /*
01347  * call-seq:
01348  *    rat.truncate               ->  integer
01349  *    rat.truncate(precision=0)  ->  rational
01350  *
01351  * Returns the truncated value (toward zero).
01352  *
01353  *    Rational(3).truncate      #=> 3
01354  *    Rational(2, 3).truncate   #=> 0
01355  *    Rational(-3, 2).truncate  #=> -1
01356  *
01357  *           decimal      -  1  2  3 . 4  5  6
01358  *                          ^  ^  ^  ^   ^  ^
01359  *          precision      -3 -2 -1  0  +1 +2
01360  *
01361  *    '%f' % Rational('-123.456').truncate(+1)  #=>  "-123.400000"
01362  *    '%f' % Rational('-123.456').truncate(-1)  #=>  "-120.000000"
01363  */
01364 static VALUE
01365 nurat_truncate_n(int argc, VALUE *argv, VALUE self)
01366 {
01367     return f_round_common(argc, argv, self, nurat_truncate);
01368 }
01369 
01370 /*
01371  * call-seq:
01372  *    rat.round               ->  integer
01373  *    rat.round(precision=0)  ->  rational
01374  *
01375  * Returns the truncated value (toward the nearest integer;
01376  * 0.5 => 1; -0.5 => -1).
01377  *
01378  *    Rational(3).round      #=> 3
01379  *    Rational(2, 3).round   #=> 1
01380  *    Rational(-3, 2).round  #=> -2
01381  *
01382  *           decimal      -  1  2  3 . 4  5  6
01383  *                          ^  ^  ^  ^   ^  ^
01384  *          precision      -3 -2 -1  0  +1 +2
01385  *
01386  *    '%f' % Rational('-123.456').round(+1)  #=> "-123.500000"
01387  *    '%f' % Rational('-123.456').round(-1)  #=> "-120.000000"
01388  */
01389 static VALUE
01390 nurat_round_n(int argc, VALUE *argv, VALUE self)
01391 {
01392     return f_round_common(argc, argv, self, nurat_round);
01393 }
01394 
01395 /*
01396  * call-seq:
01397  *    rat.to_f  ->  float
01398  *
01399  * Return the value as a float.
01400  *
01401  *    Rational(2).to_f      #=> 2.0
01402  *    Rational(9, 4).to_f   #=> 2.25
01403  *    Rational(-3, 4).to_f  #=> -0.75
01404  *    Rational(20, 3).to_f  #=> 6.666666666666667
01405  */
01406 static VALUE
01407 nurat_to_f(VALUE self)
01408 {
01409     get_dat1(self);
01410     return f_fdiv(dat->num, dat->den);
01411 }
01412 
01413 /*
01414  * call-seq:
01415  *    rat.to_r  ->  self
01416  *
01417  * Returns self.
01418  *
01419  *    Rational(2).to_r      #=> (2/1)
01420  *    Rational(-8, 6).to_r  #=> (-4/3)
01421  */
01422 static VALUE
01423 nurat_to_r(VALUE self)
01424 {
01425     return self;
01426 }
01427 
01428 #define id_ceil rb_intern("ceil")
01429 #define f_ceil(x) rb_funcall((x), id_ceil, 0)
01430 
01431 #define id_quo rb_intern("quo")
01432 #define f_quo(x,y) rb_funcall((x), id_quo, 1, (y))
01433 
01434 #define f_reciprocal(x) f_quo(ONE, (x))
01435 
01436 /*
01437   The algorithm here is the method described in CLISP.  Bruno Haible has
01438   graciously given permission to use this algorithm.  He says, "You can use
01439   it, if you present the following explanation of the algorithm."
01440 
01441   Algorithm (recursively presented):
01442     If x is a rational number, return x.
01443     If x = 0.0, return 0.
01444     If x < 0.0, return (- (rationalize (- x))).
01445     If x > 0.0:
01446       Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
01447       exponent, sign).
01448       If m = 0 or e >= 0: return x = m*2^e.
01449       Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
01450       with smallest possible numerator and denominator.
01451       Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
01452         But in this case the result will be x itself anyway, regardless of
01453         the choice of a. Therefore we can simply ignore this case.
01454       Note 2: At first, we need to consider the closed interval [a,b].
01455         but since a and b have the denominator 2^(|e|+1) whereas x itself
01456         has a denominator <= 2^|e|, we can restrict the search to the open
01457         interval (a,b).
01458       So, for given a and b (0 < a < b) we are searching a rational number
01459       y with a <= y <= b.
01460       Recursive algorithm fraction_between(a,b):
01461         c := (ceiling a)
01462         if c < b
01463           then return c       ; because a <= c < b, c integer
01464           else
01465             ; a is not integer (otherwise we would have had c = a < b)
01466             k := c-1          ; k = floor(a), k < a < b <= k+1
01467             return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
01468                               ; note 1 <= 1/(b-k) < 1/(a-k)
01469 
01470   You can see that we are actually computing a continued fraction expansion.
01471 
01472   Algorithm (iterative):
01473     If x is rational, return x.
01474     Call (integer-decode-float x). It returns a m,e,s (mantissa,
01475       exponent, sign).
01476     If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
01477     Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
01478     (positive and already in lowest terms because the denominator is a
01479     power of two and the numerator is odd).
01480     Start a continued fraction expansion
01481       p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
01482     Loop
01483       c := (ceiling a)
01484       if c >= b
01485         then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
01486              goto Loop
01487     finally partial_quotient(c).
01488     Here partial_quotient(c) denotes the iteration
01489       i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
01490     At the end, return s * (p[i]/q[i]).
01491     This rational number is already in lowest terms because
01492     p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
01493 */
01494 
01495 static void
01496 nurat_rationalize_internal(VALUE a, VALUE b, VALUE *p, VALUE *q)
01497 {
01498     VALUE c, k, t, p0, p1, p2, q0, q1, q2;
01499 
01500     p0 = ZERO;
01501     p1 = ONE;
01502     q0 = ONE;
01503     q1 = ZERO;
01504 
01505     while (1) {
01506         c = f_ceil(a);
01507         if (f_lt_p(c, b))
01508             break;
01509         k = f_sub(c, ONE);
01510         p2 = f_add(f_mul(k, p1), p0);
01511         q2 = f_add(f_mul(k, q1), q0);
01512         t = f_reciprocal(f_sub(b, k));
01513         b = f_reciprocal(f_sub(a, k));
01514         a = t;
01515         p0 = p1;
01516         q0 = q1;
01517         p1 = p2;
01518         q1 = q2;
01519     }
01520     *p = f_add(f_mul(c, p1), p0);
01521     *q = f_add(f_mul(c, q1), q0);
01522 }
01523 
01524 /*
01525  * call-seq:
01526  *    rat.rationalize       ->  self
01527  *    rat.rationalize(eps)  ->  rational
01528  *
01529  * Returns a simpler approximation of the value if the optional
01530  * argument eps is given (rat-|eps| <= result <= rat+|eps|), self
01531  * otherwise.
01532  *
01533  *    r = Rational(5033165, 16777216)
01534  *    r.rationalize                    #=> (5033165/16777216)
01535  *    r.rationalize(Rational('0.01'))  #=> (3/10)
01536  *    r.rationalize(Rational('0.1'))   #=> (1/3)
01537  */
01538 static VALUE
01539 nurat_rationalize(int argc, VALUE *argv, VALUE self)
01540 {
01541     VALUE e, a, b, p, q;
01542 
01543     if (argc == 0)
01544         return self;
01545 
01546     if (f_negative_p(self))
01547         return f_negate(nurat_rationalize(argc, argv, f_abs(self)));
01548 
01549     rb_scan_args(argc, argv, "01", &e);
01550     e = f_abs(e);
01551     a = f_sub(self, e);
01552     b = f_add(self, e);
01553 
01554     if (f_eqeq_p(a, b))
01555         return self;
01556 
01557     nurat_rationalize_internal(a, b, &p, &q);
01558     return f_rational_new2(CLASS_OF(self), p, q);
01559 }
01560 
01561 /* :nodoc: */
01562 static VALUE
01563 nurat_hash(VALUE self)
01564 {
01565     st_index_t v, h[2];
01566     VALUE n;
01567 
01568     get_dat1(self);
01569     n = rb_hash(dat->num);
01570     h[0] = NUM2LONG(n);
01571     n = rb_hash(dat->den);
01572     h[1] = NUM2LONG(n);
01573     v = rb_memhash(h, sizeof(h));
01574     return LONG2FIX(v);
01575 }
01576 
01577 static VALUE
01578 f_format(VALUE self, VALUE (*func)(VALUE))
01579 {
01580     VALUE s;
01581     get_dat1(self);
01582 
01583     s = (*func)(dat->num);
01584     rb_str_cat2(s, "/");
01585     rb_str_concat(s, (*func)(dat->den));
01586 
01587     return s;
01588 }
01589 
01590 /*
01591  * call-seq:
01592  *    rat.to_s  ->  string
01593  *
01594  * Returns the value as a string.
01595  *
01596  *    Rational(2).to_s      #=> "2/1"
01597  *    Rational(-8, 6).to_s  #=> "-4/3"
01598  *    Rational('1/2').to_s  #=> "1/2"
01599  */
01600 static VALUE
01601 nurat_to_s(VALUE self)
01602 {
01603     return f_format(self, f_to_s);
01604 }
01605 
01606 /*
01607  * call-seq:
01608  *    rat.inspect  ->  string
01609  *
01610  * Returns the value as a string for inspection.
01611  *
01612  *    Rational(2).inspect      #=> "(2/1)"
01613  *    Rational(-8, 6).inspect  #=> "(-4/3)"
01614  *    Rational('1/2').inspect  #=> "(1/2)"
01615  */
01616 static VALUE
01617 nurat_inspect(VALUE self)
01618 {
01619     VALUE s;
01620 
01621     s = rb_usascii_str_new2("(");
01622     rb_str_concat(s, f_format(self, f_inspect));
01623     rb_str_cat2(s, ")");
01624 
01625     return s;
01626 }
01627 
01628 /* :nodoc: */
01629 static VALUE
01630 nurat_dumper(VALUE self)
01631 {
01632     return self;
01633 }
01634 
01635 /* :nodoc: */
01636 static VALUE
01637 nurat_loader(VALUE self, VALUE a)
01638 {
01639     get_dat1(self);
01640 
01641     dat->num = rb_ivar_get(a, id_i_num);
01642     dat->den = rb_ivar_get(a, id_i_den);
01643 
01644     return self;
01645 }
01646 
01647 /* :nodoc: */
01648 static VALUE
01649 nurat_marshal_dump(VALUE self)
01650 {
01651     VALUE a;
01652     get_dat1(self);
01653 
01654     a = rb_assoc_new(dat->num, dat->den);
01655     rb_copy_generic_ivar(a, self);
01656     return a;
01657 }
01658 
01659 /* :nodoc: */
01660 static VALUE
01661 nurat_marshal_load(VALUE self, VALUE a)
01662 {
01663     rb_check_frozen(self);
01664     rb_check_trusted(self);
01665 
01666     Check_Type(a, T_ARRAY);
01667     if (RARRAY_LEN(a) != 2)
01668         rb_raise(rb_eArgError, "marshaled rational must have an array whose length is 2 but %ld", RARRAY_LEN(a));
01669     if (f_zero_p(RARRAY_PTR(a)[1]))
01670         rb_raise_zerodiv();
01671 
01672     rb_ivar_set(self, id_i_num, RARRAY_PTR(a)[0]);
01673     rb_ivar_set(self, id_i_den, RARRAY_PTR(a)[1]);
01674 
01675     return self;
01676 }
01677 
01678 /* --- */
01679 
01680 VALUE
01681 rb_rational_reciprocal(VALUE x)
01682 {
01683     get_dat1(x);
01684     return f_rational_new_no_reduce2(CLASS_OF(x), dat->den, dat->num);
01685 }
01686 
01687 /*
01688  * call-seq:
01689  *    int.gcd(int2)  ->  integer
01690  *
01691  * Returns the greatest common divisor (always positive).  0.gcd(x)
01692  * and x.gcd(0) return abs(x).
01693  *
01694  *    2.gcd(2)                    #=> 2
01695  *    3.gcd(-7)                   #=> 1
01696  *    ((1<<31)-1).gcd((1<<61)-1)  #=> 1
01697  */
01698 VALUE
01699 rb_gcd(VALUE self, VALUE other)
01700 {
01701     other = nurat_int_value(other);
01702     return f_gcd(self, other);
01703 }
01704 
01705 /*
01706  * call-seq:
01707  *    int.lcm(int2)  ->  integer
01708  *
01709  * Returns the least common multiple (always positive).  0.lcm(x) and
01710  * x.lcm(0) return zero.
01711  *
01712  *    2.lcm(2)                    #=> 2
01713  *    3.lcm(-7)                   #=> 21
01714  *    ((1<<31)-1).lcm((1<<61)-1)  #=> 4951760154835678088235319297
01715  */
01716 VALUE
01717 rb_lcm(VALUE self, VALUE other)
01718 {
01719     other = nurat_int_value(other);
01720     return f_lcm(self, other);
01721 }
01722 
01723 /*
01724  * call-seq:
01725  *    int.gcdlcm(int2)  ->  array
01726  *
01727  * Returns an array; [int.gcd(int2), int.lcm(int2)].
01728  *
01729  *    2.gcdlcm(2)                    #=> [2, 2]
01730  *    3.gcdlcm(-7)                   #=> [1, 21]
01731  *    ((1<<31)-1).gcdlcm((1<<61)-1)  #=> [1, 4951760154835678088235319297]
01732  */
01733 VALUE
01734 rb_gcdlcm(VALUE self, VALUE other)
01735 {
01736     other = nurat_int_value(other);
01737     return rb_assoc_new(f_gcd(self, other), f_lcm(self, other));
01738 }
01739 
01740 VALUE
01741 rb_rational_raw(VALUE x, VALUE y)
01742 {
01743     return nurat_s_new_internal(rb_cRational, x, y);
01744 }
01745 
01746 VALUE
01747 rb_rational_new(VALUE x, VALUE y)
01748 {
01749     return nurat_s_canonicalize_internal(rb_cRational, x, y);
01750 }
01751 
01752 static VALUE nurat_s_convert(int argc, VALUE *argv, VALUE klass);
01753 
01754 VALUE
01755 rb_Rational(VALUE x, VALUE y)
01756 {
01757     VALUE a[2];
01758     a[0] = x;
01759     a[1] = y;
01760     return nurat_s_convert(2, a, rb_cRational);
01761 }
01762 
01763 #define id_numerator rb_intern("numerator")
01764 #define f_numerator(x) rb_funcall((x), id_numerator, 0)
01765 
01766 #define id_denominator rb_intern("denominator")
01767 #define f_denominator(x) rb_funcall((x), id_denominator, 0)
01768 
01769 #define id_to_r rb_intern("to_r")
01770 #define f_to_r(x) rb_funcall((x), id_to_r, 0)
01771 
01772 /*
01773  * call-seq:
01774  *    num.numerator  ->  integer
01775  *
01776  * Returns the numerator.
01777  */
01778 static VALUE
01779 numeric_numerator(VALUE self)
01780 {
01781     return f_numerator(f_to_r(self));
01782 }
01783 
01784 /*
01785  * call-seq:
01786  *    num.denominator  ->  integer
01787  *
01788  * Returns the denominator (always positive).
01789  */
01790 static VALUE
01791 numeric_denominator(VALUE self)
01792 {
01793     return f_denominator(f_to_r(self));
01794 }
01795 
01796 /*
01797  * call-seq:
01798  *    int.numerator  ->  self
01799  *
01800  * Returns self.
01801  */
01802 static VALUE
01803 integer_numerator(VALUE self)
01804 {
01805     return self;
01806 }
01807 
01808 /*
01809  * call-seq:
01810  *    int.denominator  ->  1
01811  *
01812  * Returns 1.
01813  */
01814 static VALUE
01815 integer_denominator(VALUE self)
01816 {
01817     return INT2FIX(1);
01818 }
01819 
01820 /*
01821  * call-seq:
01822  *    flo.numerator  ->  integer
01823  *
01824  * Returns the numerator.  The result is machine dependent.
01825  *
01826  *    n = 0.3.numerator    #=> 5404319552844595
01827  *    d = 0.3.denominator  #=> 18014398509481984
01828  *    n.fdiv(d)            #=> 0.3
01829  */
01830 static VALUE
01831 float_numerator(VALUE self)
01832 {
01833     double d = RFLOAT_VALUE(self);
01834     if (isinf(d) || isnan(d))
01835         return self;
01836     return rb_call_super(0, 0);
01837 }
01838 
01839 /*
01840  * call-seq:
01841  *    flo.denominator  ->  integer
01842  *
01843  * Returns the denominator (always positive).  The result is machine
01844  * dependent.
01845  *
01846  * See numerator.
01847  */
01848 static VALUE
01849 float_denominator(VALUE self)
01850 {
01851     double d = RFLOAT_VALUE(self);
01852     if (isinf(d) || isnan(d))
01853         return INT2FIX(1);
01854     return rb_call_super(0, 0);
01855 }
01856 
01857 /*
01858  * call-seq:
01859  *    nil.to_r  ->  (0/1)
01860  *
01861  * Returns zero as a rational.
01862  */
01863 static VALUE
01864 nilclass_to_r(VALUE self)
01865 {
01866     return rb_rational_new1(INT2FIX(0));
01867 }
01868 
01869 /*
01870  * call-seq:
01871  *    nil.rationalize([eps])  ->  (0/1)
01872  *
01873  * Returns zero as a rational.  The optional argument eps is always
01874  * ignored.
01875  */
01876 static VALUE
01877 nilclass_rationalize(int argc, VALUE *argv, VALUE self)
01878 {
01879     rb_scan_args(argc, argv, "01", NULL);
01880     return nilclass_to_r(self);
01881 }
01882 
01883 /*
01884  * call-seq:
01885  *    int.to_r  ->  rational
01886  *
01887  * Returns the value as a rational.
01888  *
01889  *    1.to_r        #=> (1/1)
01890  *    (1<<64).to_r  #=> (18446744073709551616/1)
01891  */
01892 static VALUE
01893 integer_to_r(VALUE self)
01894 {
01895     return rb_rational_new1(self);
01896 }
01897 
01898 /*
01899  * call-seq:
01900  *    int.rationalize([eps])  ->  rational
01901  *
01902  * Returns the value as a rational.  The optional argument eps is
01903  * always ignored.
01904  */
01905 static VALUE
01906 integer_rationalize(int argc, VALUE *argv, VALUE self)
01907 {
01908     rb_scan_args(argc, argv, "01", NULL);
01909     return integer_to_r(self);
01910 }
01911 
01912 static void
01913 float_decode_internal(VALUE self, VALUE *rf, VALUE *rn)
01914 {
01915     double f;
01916     int n;
01917 
01918     f = frexp(RFLOAT_VALUE(self), &n);
01919     f = ldexp(f, DBL_MANT_DIG);
01920     n -= DBL_MANT_DIG;
01921     *rf = rb_dbl2big(f);
01922     *rn = INT2FIX(n);
01923 }
01924 
01925 #if 0
01926 static VALUE
01927 float_decode(VALUE self)
01928 {
01929     VALUE f, n;
01930 
01931     float_decode_internal(self, &f, &n);
01932     return rb_assoc_new(f, n);
01933 }
01934 #endif
01935 
01936 #define id_lshift rb_intern("<<")
01937 #define f_lshift(x,n) rb_funcall((x), id_lshift, 1, (n))
01938 
01939 /*
01940  * call-seq:
01941  *    flt.to_r  ->  rational
01942  *
01943  * Returns the value as a rational.
01944  *
01945  * NOTE: 0.3.to_r isn't the same as '0.3'.to_r.  The latter is
01946  * equivalent to '3/10'.to_r, but the former isn't so.
01947  *
01948  *    2.0.to_r    #=> (2/1)
01949  *    2.5.to_r    #=> (5/2)
01950  *    -0.75.to_r  #=> (-3/4)
01951  *    0.0.to_r    #=> (0/1)
01952  *
01953  * See rationalize.
01954  */
01955 static VALUE
01956 float_to_r(VALUE self)
01957 {
01958     VALUE f, n;
01959 
01960     float_decode_internal(self, &f, &n);
01961 #if FLT_RADIX == 2
01962     {
01963         long ln = FIX2LONG(n);
01964 
01965         if (ln == 0)
01966             return f_to_r(f);
01967         if (ln > 0)
01968             return f_to_r(f_lshift(f, n));
01969         ln = -ln;
01970         return rb_rational_new2(f, f_lshift(ONE, INT2FIX(ln)));
01971     }
01972 #else
01973     return f_to_r(f_mul(f, f_expt(INT2FIX(FLT_RADIX), n)));
01974 #endif
01975 }
01976 
01977 /*
01978  * call-seq:
01979  *    flt.rationalize([eps])  ->  rational
01980  *
01981  * Returns a simpler approximation of the value (flt-|eps| <= result
01982  * <= flt+|eps|).  if the optional eps is not given, it will be chosen
01983  * automatically.
01984  *
01985  *    0.3.rationalize          #=> (3/10)
01986  *    1.333.rationalize        #=> (1333/1000)
01987  *    1.333.rationalize(0.01)  #=> (4/3)
01988  *
01989  * See to_r.
01990  */
01991 static VALUE
01992 float_rationalize(int argc, VALUE *argv, VALUE self)
01993 {
01994     VALUE e, a, b, p, q;
01995 
01996     if (f_negative_p(self))
01997         return f_negate(float_rationalize(argc, argv, f_abs(self)));
01998 
01999     rb_scan_args(argc, argv, "01", &e);
02000 
02001     if (argc != 0) {
02002         e = f_abs(e);
02003         a = f_sub(self, e);
02004         b = f_add(self, e);
02005     }
02006     else {
02007         VALUE f, n;
02008 
02009         float_decode_internal(self, &f, &n);
02010         if (f_zero_p(f) || f_positive_p(n))
02011             return rb_rational_new1(f_lshift(f, n));
02012 
02013 #if FLT_RADIX == 2
02014         {
02015             VALUE two_times_f, den;
02016 
02017             two_times_f = f_mul(TWO, f);
02018             den = f_lshift(ONE, f_sub(ONE, n));
02019 
02020             a = rb_rational_new2(f_sub(two_times_f, ONE), den);
02021             b = rb_rational_new2(f_add(two_times_f, ONE), den);
02022         }
02023 #else
02024         {
02025             VALUE radix_times_f, den;
02026 
02027             radix_times_f = f_mul(INT2FIX(FLT_RADIX), f);
02028             den = f_expt(INT2FIX(FLT_RADIX), f_sub(ONE, n));
02029 
02030             a = rb_rational_new2(f_sub(radix_times_f, INT2FIX(FLT_RADIX - 1)), den);
02031             b = rb_rational_new2(f_add(radix_times_f, INT2FIX(FLT_RADIX - 1)), den);
02032         }
02033 #endif
02034     }
02035 
02036     if (f_eqeq_p(a, b))
02037         return f_to_r(self);
02038 
02039     nurat_rationalize_internal(a, b, &p, &q);
02040     return rb_rational_new2(p, q);
02041 }
02042 
02043 #include <ctype.h>
02044 
02045 inline static int
02046 issign(int c)
02047 {
02048     return (c == '-' || c == '+');
02049 }
02050 
02051 static int
02052 read_sign(const char **s)
02053 {
02054     int sign = '?';
02055 
02056     if (issign(**s)) {
02057         sign = **s;
02058         (*s)++;
02059     }
02060     return sign;
02061 }
02062 
02063 inline static int
02064 isdecimal(int c)
02065 {
02066     return isdigit((unsigned char)c);
02067 }
02068 
02069 static int
02070 read_digits(const char **s, int strict,
02071             VALUE *num, int *count)
02072 {
02073     char *b, *bb;
02074     int us = 1, ret = 1;
02075     VALUE tmp;
02076 
02077     if (!isdecimal(**s)) {
02078         *num = ZERO;
02079         return 0;
02080     }
02081 
02082     bb = b = ALLOCV_N(char, tmp, strlen(*s) + 1);
02083 
02084     while (isdecimal(**s) || **s == '_') {
02085         if (**s == '_') {
02086             if (strict) {
02087                 if (us) {
02088                     ret = 0;
02089                     goto conv;
02090                 }
02091             }
02092             us = 1;
02093         }
02094         else {
02095             if (count)
02096                 (*count)++;
02097             *b++ = **s;
02098             us = 0;
02099         }
02100         (*s)++;
02101     }
02102     if (us)
02103         do {
02104             (*s)--;
02105         } while (**s == '_');
02106   conv:
02107     *b = '\0';
02108     *num = rb_cstr_to_inum(bb, 10, 0);
02109     ALLOCV_END(tmp);
02110     return ret;
02111 }
02112 
02113 inline static int
02114 islettere(int c)
02115 {
02116     return (c == 'e' || c == 'E');
02117 }
02118 
02119 static int
02120 read_num(const char **s, int numsign, int strict,
02121          VALUE *num)
02122 {
02123     VALUE ip, fp, exp;
02124 
02125     *num = rb_rational_new2(ZERO, ONE);
02126     exp = Qnil;
02127 
02128     if (**s != '.') {
02129         if (!read_digits(s, strict, &ip, NULL))
02130             return 0;
02131         *num = rb_rational_new2(ip, ONE);
02132     }
02133 
02134     if (**s == '.') {
02135         int count = 0;
02136 
02137         (*s)++;
02138         if (!read_digits(s, strict, &fp, &count))
02139             return 0;
02140         {
02141             VALUE l = f_expt10(INT2NUM(count));
02142             *num = f_mul(*num, l);
02143             *num = f_add(*num, fp);
02144             *num = f_div(*num, l);
02145         }
02146     }
02147 
02148     if (islettere(**s)) {
02149         int expsign;
02150 
02151         (*s)++;
02152         expsign = read_sign(s);
02153         if (!read_digits(s, strict, &exp, NULL))
02154             return 0;
02155         if (expsign == '-')
02156             exp = f_negate(exp);
02157     }
02158 
02159     if (numsign == '-')
02160         *num = f_negate(*num);
02161     if (!NIL_P(exp)) {
02162         VALUE l = f_expt10(exp);
02163         *num = f_mul(*num, l);
02164     }
02165     return 1;
02166 }
02167 
02168 inline static int
02169 read_den(const char **s, int strict,
02170          VALUE *num)
02171 {
02172     if (!read_digits(s, strict, num, NULL))
02173         return 0;
02174     return 1;
02175 }
02176 
02177 static int
02178 read_rat_nos(const char **s, int sign, int strict,
02179              VALUE *num)
02180 {
02181     VALUE den;
02182 
02183     if (!read_num(s, sign, strict, num))
02184         return 0;
02185     if (**s == '/') {
02186         (*s)++;
02187         if (!read_den(s, strict, &den))
02188             return 0;
02189         if (!(FIXNUM_P(den) && FIX2LONG(den) == 1))
02190             *num = f_div(*num, den);
02191     }
02192     return 1;
02193 }
02194 
02195 static int
02196 read_rat(const char **s, int strict,
02197          VALUE *num)
02198 {
02199     int sign;
02200 
02201     sign = read_sign(s);
02202     if (!read_rat_nos(s, sign, strict, num))
02203         return 0;
02204     return 1;
02205 }
02206 
02207 inline static void
02208 skip_ws(const char **s)
02209 {
02210     while (isspace((unsigned char)**s))
02211         (*s)++;
02212 }
02213 
02214 static int
02215 parse_rat(const char *s, int strict,
02216           VALUE *num)
02217 {
02218     skip_ws(&s);
02219     if (!read_rat(&s, strict, num))
02220         return 0;
02221     skip_ws(&s);
02222 
02223     if (strict)
02224         if (*s != '\0')
02225             return 0;
02226     return 1;
02227 }
02228 
02229 static VALUE
02230 string_to_r_strict(VALUE self)
02231 {
02232     char *s;
02233     VALUE num;
02234 
02235     rb_must_asciicompat(self);
02236 
02237     s = RSTRING_PTR(self);
02238 
02239     if (!s || memchr(s, '\0', RSTRING_LEN(self)))
02240         rb_raise(rb_eArgError, "string contains null byte");
02241 
02242     if (s && s[RSTRING_LEN(self)]) {
02243         rb_str_modify(self);
02244         s = RSTRING_PTR(self);
02245         s[RSTRING_LEN(self)] = '\0';
02246     }
02247 
02248     if (!s)
02249         s = (char *)"";
02250 
02251     if (!parse_rat(s, 1, &num)) {
02252         VALUE ins = f_inspect(self);
02253         rb_raise(rb_eArgError, "invalid value for convert(): %s",
02254                  StringValuePtr(ins));
02255     }
02256 
02257     if (RB_TYPE_P(num, T_FLOAT))
02258         rb_raise(rb_eFloatDomainError, "Infinity");
02259     return num;
02260 }
02261 
02262 /*
02263  * call-seq:
02264  *    str.to_r  ->  rational
02265  *
02266  * Returns a rational which denotes the string form.  The parser
02267  * ignores leading whitespaces and trailing garbage.  Any digit
02268  * sequences can be separated by an underscore.  Returns zero for null
02269  * or garbage string.
02270  *
02271  * NOTE: '0.3'.to_r isn't the same as 0.3.to_r.  The former is
02272  * equivalent to '3/10'.to_r, but the latter isn't so.
02273  *
02274  *    '  2  '.to_r       #=> (2/1)
02275  *    '300/2'.to_r       #=> (150/1)
02276  *    '-9.2'.to_r        #=> (-46/5)
02277  *    '-9.2e2'.to_r      #=> (-920/1)
02278  *    '1_234_567'.to_r   #=> (1234567/1)
02279  *    '21 june 09'.to_r  #=> (21/1)
02280  *    '21/06/09'.to_r    #=> (7/2)
02281  *    'bwv 1079'.to_r    #=> (0/1)
02282  *
02283  * See Kernel.Rational.
02284  */
02285 static VALUE
02286 string_to_r(VALUE self)
02287 {
02288     char *s;
02289     VALUE num;
02290 
02291     rb_must_asciicompat(self);
02292 
02293     s = RSTRING_PTR(self);
02294 
02295     if (s && s[RSTRING_LEN(self)]) {
02296         rb_str_modify(self);
02297         s = RSTRING_PTR(self);
02298         s[RSTRING_LEN(self)] = '\0';
02299     }
02300 
02301     if (!s)
02302         s = (char *)"";
02303 
02304     (void)parse_rat(s, 0, &num);
02305 
02306     if (RB_TYPE_P(num, T_FLOAT))
02307         rb_raise(rb_eFloatDomainError, "Infinity");
02308     return num;
02309 }
02310 
02311 VALUE
02312 rb_cstr_to_rat(const char *s, int strict) /* for complex's internal */
02313 {
02314     VALUE num;
02315 
02316     (void)parse_rat(s, strict, &num);
02317 
02318     if (RB_TYPE_P(num, T_FLOAT))
02319         rb_raise(rb_eFloatDomainError, "Infinity");
02320     return num;
02321 }
02322 
02323 static VALUE
02324 nurat_s_convert(int argc, VALUE *argv, VALUE klass)
02325 {
02326     VALUE a1, a2, backref;
02327 
02328     rb_scan_args(argc, argv, "11", &a1, &a2);
02329 
02330     if (NIL_P(a1) || (argc == 2 && NIL_P(a2)))
02331         rb_raise(rb_eTypeError, "can't convert nil into Rational");
02332 
02333     switch (TYPE(a1)) {
02334       case T_COMPLEX:
02335         if (k_exact_zero_p(RCOMPLEX(a1)->imag))
02336             a1 = RCOMPLEX(a1)->real;
02337     }
02338 
02339     switch (TYPE(a2)) {
02340       case T_COMPLEX:
02341         if (k_exact_zero_p(RCOMPLEX(a2)->imag))
02342             a2 = RCOMPLEX(a2)->real;
02343     }
02344 
02345     backref = rb_backref_get();
02346     rb_match_busy(backref);
02347 
02348     switch (TYPE(a1)) {
02349       case T_FIXNUM:
02350       case T_BIGNUM:
02351         break;
02352       case T_FLOAT:
02353         a1 = f_to_r(a1);
02354         break;
02355       case T_STRING:
02356         a1 = string_to_r_strict(a1);
02357         break;
02358     }
02359 
02360     switch (TYPE(a2)) {
02361       case T_FIXNUM:
02362       case T_BIGNUM:
02363         break;
02364       case T_FLOAT:
02365         a2 = f_to_r(a2);
02366         break;
02367       case T_STRING:
02368         a2 = string_to_r_strict(a2);
02369         break;
02370     }
02371 
02372     rb_backref_set(backref);
02373 
02374     switch (TYPE(a1)) {
02375       case T_RATIONAL:
02376         if (argc == 1 || (k_exact_one_p(a2)))
02377             return a1;
02378     }
02379 
02380     if (argc == 1) {
02381         if (!(k_numeric_p(a1) && k_integer_p(a1)))
02382             return rb_convert_type(a1, T_RATIONAL, "Rational", "to_r");
02383     }
02384     else {
02385         if ((k_numeric_p(a1) && k_numeric_p(a2)) &&
02386             (!f_integer_p(a1) || !f_integer_p(a2)))
02387             return f_div(a1, a2);
02388     }
02389 
02390     {
02391         VALUE argv2[2];
02392         argv2[0] = a1;
02393         argv2[1] = a2;
02394         return nurat_s_new(argc, argv2, klass);
02395     }
02396 }
02397 
02398 /*
02399  * A rational number can be represented as a paired integer number;
02400  * a/b (b>0).  Where a is numerator and b is denominator.  Integer a
02401  * equals rational a/1 mathematically.
02402  *
02403  * In ruby, you can create rational object with Rational, to_r or
02404  * rationalize method.  The return values will be irreducible.
02405  *
02406  *    Rational(1)      #=> (1/1)
02407  *    Rational(2, 3)   #=> (2/3)
02408  *    Rational(4, -6)  #=> (-2/3)
02409  *    3.to_r           #=> (3/1)
02410  *
02411  * You can also create rational object from floating-point numbers or
02412  * strings.
02413  *
02414  *    Rational(0.3)    #=> (5404319552844595/18014398509481984)
02415  *    Rational('0.3')  #=> (3/10)
02416  *    Rational('2/3')  #=> (2/3)
02417  *
02418  *    0.3.to_r         #=> (5404319552844595/18014398509481984)
02419  *    '0.3'.to_r       #=> (3/10)
02420  *    '2/3'.to_r       #=> (2/3)
02421  *    0.3.rationalize  #=> (3/10)
02422  *
02423  * A rational object is an exact number, which helps you to write
02424  * program without any rounding errors.
02425  *
02426  *    10.times.inject(0){|t,| t + 0.1}              #=> 0.9999999999999999
02427  *    10.times.inject(0){|t,| t + Rational('0.1')}  #=> (1/1)
02428  *
02429  * However, when an expression has inexact factor (numerical value or
02430  * operation), will produce an inexact result.
02431  *
02432  *    Rational(10) / 3   #=> (10/3)
02433  *    Rational(10) / 3.0 #=> 3.3333333333333335
02434  *
02435  *    Rational(-8) ** Rational(1, 3)
02436  *                       #=> (1.0000000000000002+1.7320508075688772i)
02437  */
02438 void
02439 Init_Rational(void)
02440 {
02441     VALUE compat;
02442 #undef rb_intern
02443 #define rb_intern(str) rb_intern_const(str)
02444 
02445     assert(fprintf(stderr, "assert() is now active\n"));
02446 
02447     id_abs = rb_intern("abs");
02448     id_cmp = rb_intern("<=>");
02449     id_convert = rb_intern("convert");
02450     id_eqeq_p = rb_intern("==");
02451     id_expt = rb_intern("**");
02452     id_fdiv = rb_intern("fdiv");
02453     id_floor = rb_intern("floor");
02454     id_idiv = rb_intern("div");
02455     id_inspect = rb_intern("inspect");
02456     id_integer_p = rb_intern("integer?");
02457     id_negate = rb_intern("-@");
02458     id_to_f = rb_intern("to_f");
02459     id_to_i = rb_intern("to_i");
02460     id_to_s = rb_intern("to_s");
02461     id_truncate = rb_intern("truncate");
02462     id_i_num = rb_intern("@numerator");
02463     id_i_den = rb_intern("@denominator");
02464 
02465     rb_cRational = rb_define_class("Rational", rb_cNumeric);
02466 
02467     rb_define_alloc_func(rb_cRational, nurat_s_alloc);
02468     rb_undef_method(CLASS_OF(rb_cRational), "allocate");
02469 
02470 #if 0
02471     rb_define_private_method(CLASS_OF(rb_cRational), "new!", nurat_s_new_bang, -1);
02472     rb_define_private_method(CLASS_OF(rb_cRational), "new", nurat_s_new, -1);
02473 #else
02474     rb_undef_method(CLASS_OF(rb_cRational), "new");
02475 #endif
02476 
02477     rb_define_global_function("Rational", nurat_f_rational, -1);
02478 
02479     rb_define_method(rb_cRational, "numerator", nurat_numerator, 0);
02480     rb_define_method(rb_cRational, "denominator", nurat_denominator, 0);
02481 
02482     rb_define_method(rb_cRational, "+", nurat_add, 1);
02483     rb_define_method(rb_cRational, "-", nurat_sub, 1);
02484     rb_define_method(rb_cRational, "*", nurat_mul, 1);
02485     rb_define_method(rb_cRational, "/", nurat_div, 1);
02486     rb_define_method(rb_cRational, "quo", nurat_div, 1);
02487     rb_define_method(rb_cRational, "fdiv", nurat_fdiv, 1);
02488     rb_define_method(rb_cRational, "**", nurat_expt, 1);
02489 
02490     rb_define_method(rb_cRational, "<=>", nurat_cmp, 1);
02491     rb_define_method(rb_cRational, "==", nurat_eqeq_p, 1);
02492     rb_define_method(rb_cRational, "coerce", nurat_coerce, 1);
02493 
02494 #if 0 /* NUBY */
02495     rb_define_method(rb_cRational, "//", nurat_idiv, 1);
02496 #endif
02497 
02498 #if 0
02499     rb_define_method(rb_cRational, "quot", nurat_quot, 1);
02500     rb_define_method(rb_cRational, "quotrem", nurat_quotrem, 1);
02501 #endif
02502 
02503 #if 0
02504     rb_define_method(rb_cRational, "rational?", nurat_true, 0);
02505     rb_define_method(rb_cRational, "exact?", nurat_true, 0);
02506 #endif
02507 
02508     rb_define_method(rb_cRational, "floor", nurat_floor_n, -1);
02509     rb_define_method(rb_cRational, "ceil", nurat_ceil_n, -1);
02510     rb_define_method(rb_cRational, "truncate", nurat_truncate_n, -1);
02511     rb_define_method(rb_cRational, "round", nurat_round_n, -1);
02512 
02513     rb_define_method(rb_cRational, "to_i", nurat_truncate, 0);
02514     rb_define_method(rb_cRational, "to_f", nurat_to_f, 0);
02515     rb_define_method(rb_cRational, "to_r", nurat_to_r, 0);
02516     rb_define_method(rb_cRational, "rationalize", nurat_rationalize, -1);
02517 
02518     rb_define_method(rb_cRational, "hash", nurat_hash, 0);
02519 
02520     rb_define_method(rb_cRational, "to_s", nurat_to_s, 0);
02521     rb_define_method(rb_cRational, "inspect", nurat_inspect, 0);
02522 
02523     rb_define_private_method(rb_cRational, "marshal_dump", nurat_marshal_dump, 0);
02524     compat = rb_define_class_under(rb_cRational, "compatible", rb_cObject);
02525     rb_define_private_method(compat, "marshal_load", nurat_marshal_load, 1);
02526     rb_marshal_define_compat(rb_cRational, compat, nurat_dumper, nurat_loader);
02527 
02528     /* --- */
02529 
02530     rb_define_method(rb_cInteger, "gcd", rb_gcd, 1);
02531     rb_define_method(rb_cInteger, "lcm", rb_lcm, 1);
02532     rb_define_method(rb_cInteger, "gcdlcm", rb_gcdlcm, 1);
02533 
02534     rb_define_method(rb_cNumeric, "numerator", numeric_numerator, 0);
02535     rb_define_method(rb_cNumeric, "denominator", numeric_denominator, 0);
02536 
02537     rb_define_method(rb_cInteger, "numerator", integer_numerator, 0);
02538     rb_define_method(rb_cInteger, "denominator", integer_denominator, 0);
02539 
02540     rb_define_method(rb_cFloat, "numerator", float_numerator, 0);
02541     rb_define_method(rb_cFloat, "denominator", float_denominator, 0);
02542 
02543     rb_define_method(rb_cNilClass, "to_r", nilclass_to_r, 0);
02544     rb_define_method(rb_cNilClass, "rationalize", nilclass_rationalize, -1);
02545     rb_define_method(rb_cInteger, "to_r", integer_to_r, 0);
02546     rb_define_method(rb_cInteger, "rationalize", integer_rationalize, -1);
02547     rb_define_method(rb_cFloat, "to_r", float_to_r, 0);
02548     rb_define_method(rb_cFloat, "rationalize", float_rationalize, -1);
02549 
02550     rb_define_method(rb_cString, "to_r", string_to_r, 0);
02551 
02552     rb_define_private_method(CLASS_OF(rb_cRational), "convert", nurat_s_convert, -1);
02553 }
02554 
02555 /*
02556 Local variables:
02557 c-file-style: "ruby"
02558 End:
02559 */
02560