Back to index

plt-scheme  4.2.1
rational.c
Go to the documentation of this file.
00001 /*
00002   MzScheme
00003   Copyright (c) 2004-2009 PLT Scheme Inc.
00004   Copyright (c) 1995-2001 Matthew Flatt
00005  
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Library General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Library General Public License for more details.
00015 
00016     You should have received a copy of the GNU Library General Public
00017     License along with this library; if not, write to the Free
00018     Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
00019     Boston, MA 02110-1301 USA.
00020 
00021   libscheme
00022   Copyright (c) 1994 Brent Benson
00023   All rights reserved.
00024 */
00025 
00026 #include "schpriv.h"
00027 #include <ctype.h>
00028 #include <math.h>
00029 
00030 static Scheme_Object *one = scheme_make_integer(1);
00031 
00032 static Scheme_Object *make_rational(const Scheme_Object *n, const Scheme_Object *d,
00033                                 int normalize)
00034 {
00035   Scheme_Rational *r;
00036 
00037   r = (Scheme_Rational *)scheme_malloc_small_dirty_tagged(sizeof(Scheme_Rational));
00038   r->so.type = scheme_rational_type;
00039   CLEAR_KEY_FIELD(&r->so);
00040   r->num = (Scheme_Object *)n;
00041   r->denom = (Scheme_Object *)d;
00042   
00043   return (normalize 
00044          ? scheme_rational_normalize((Scheme_Object *)r) 
00045          : (Scheme_Object *)r);
00046 }
00047 
00048 Scheme_Object *scheme_make_rational(const Scheme_Object *n, const Scheme_Object *d)
00049 {
00050   return make_rational(scheme_bignum_normalize(n), 
00051                      scheme_bignum_normalize(d), 1);
00052 }
00053 
00054 Scheme_Object *scheme_integer_to_rational(const Scheme_Object *n)
00055 {
00056   return make_rational(n, one, 0);
00057 }
00058 
00059 #ifdef MZ_XFORM
00060 START_XFORM_SKIP;
00061 #endif
00062 
00063 Scheme_Object *scheme_make_small_rational(long n, Small_Rational *s)
00064 {
00065   s->so.type = scheme_rational_type;
00066   s->num = scheme_make_integer(n);
00067   s->denom = one;
00068 
00069   return (Scheme_Object *)s;
00070 }
00071 
00072 Scheme_Object *scheme_make_small_bn_rational(Scheme_Object *n, Small_Rational *s)
00073 {
00074   s->so.type = scheme_rational_type;
00075   s->num = n;
00076   s->denom = one;
00077 
00078   return (Scheme_Object *)s;
00079 }
00080 
00081 #ifdef MZ_XFORM
00082 END_XFORM_SKIP;
00083 #endif
00084 
00085 int scheme_is_rational_positive(const Scheme_Object *o)
00086 {
00087   Scheme_Rational *r = (Scheme_Rational *)o;
00088 
00089   if (SCHEME_INTP(r->num))
00090     return (SCHEME_INT_VAL(r->num) > 0);
00091   else 
00092     return SCHEME_BIGPOS(r->num);
00093 }
00094 
00095 Scheme_Object *scheme_rational_normalize(const Scheme_Object *o)
00096 {
00097   Scheme_Rational *r = (Scheme_Rational *)o;
00098   Scheme_Object *gcd, *tmpn;
00099   int negate = 0;
00100 
00101   if (r->num == scheme_exact_zero)
00102     return scheme_make_integer(0);
00103 
00104   if (SCHEME_INTP(r->denom)) {
00105     if (SCHEME_INT_VAL(r->denom) < 0) {
00106       tmpn = scheme_make_integer_value(-SCHEME_INT_VAL(r->denom));
00107       r->denom = tmpn;
00108       negate = 1;
00109     }
00110   } else if (!SCHEME_BIGPOS(r->denom)) {
00111     tmpn = scheme_bignum_negate(r->denom);
00112     r->denom = tmpn;
00113     negate = 1;
00114   }
00115 
00116   if (negate) {
00117     if (SCHEME_INTP(r->num)) {
00118       tmpn = scheme_make_integer_value(-SCHEME_INT_VAL(r->num));
00119       r->num = tmpn;
00120     } else {
00121       tmpn = scheme_bignum_negate(r->num);
00122       r->num = tmpn;
00123     }
00124   }
00125   
00126   if (r->denom == one)
00127     return r->num;
00128 
00129   gcd = scheme_bin_gcd(r->num, r->denom);
00130 
00131   if (gcd == one)
00132     return (Scheme_Object *)o;
00133 
00134   tmpn = scheme_bin_quotient(r->num, gcd);
00135   r->num = tmpn;
00136   tmpn = scheme_bin_quotient(r->denom, gcd);
00137   r->denom = tmpn;
00138 
00139   if (r->denom == one)
00140     return r->num;
00141 
00142   return (Scheme_Object *)r;
00143 }
00144 
00145 Scheme_Object *scheme_rational_numerator(const Scheme_Object *n)
00146 {
00147   return ((Scheme_Rational *)n)->num;
00148 }
00149 
00150 Scheme_Object *scheme_rational_denominator(const Scheme_Object *n)
00151 {
00152   return ((Scheme_Rational *)n)->denom;
00153 }
00154 
00155 Scheme_Object *scheme_make_fixnum_rational(long n, long d)
00156 {
00157   /* This function is called to implement division on small integers,
00158      so don't allocate unless necessary. */
00159   Small_Rational s;
00160   Scheme_Object *o;
00161   
00162   s.so.type = scheme_rational_type;
00163   s.num = scheme_make_integer(n);
00164   s.denom = scheme_make_integer(d);
00165 
00166   o = scheme_rational_normalize((Scheme_Object *)&s);
00167   if (o == (Scheme_Object *)&s)
00168     return make_rational(s.num, s.denom, 0);
00169   else
00170     return o;
00171 }
00172 
00173 int scheme_rational_eq(const Scheme_Object *a, const Scheme_Object *b)
00174 {
00175   Scheme_Rational *ra = (Scheme_Rational *)a;
00176   Scheme_Rational *rb = (Scheme_Rational *)b;
00177 
00178   if (SCHEME_INTP(ra->num) && SCHEME_INTP(rb->num)) {
00179     if (ra->num != rb->num)
00180       return 0;
00181   } else if (SCHEME_BIGNUMP(ra->num) && SCHEME_BIGNUMP(rb->num)) {
00182     if (!scheme_bignum_eq(ra->num, rb->num))
00183       return 0;
00184   } else
00185     return 0;
00186 
00187   if (SCHEME_INTP(ra->denom) && SCHEME_INTP(rb->denom)) {
00188     if (ra->denom != rb->denom)
00189       return 0;
00190   } else if (SCHEME_BIGNUMP(ra->denom) && SCHEME_BIGNUMP(rb->denom)) {
00191     if (!scheme_bignum_eq(ra->denom, rb->denom))
00192       return 0;
00193   } else
00194     return 0;
00195 
00196   return 1;
00197 }
00198 
00199 static int rational_lt(const Scheme_Object *a, const Scheme_Object *b, int or_eq)
00200 {
00201   Scheme_Rational *ra = (Scheme_Rational *)a;
00202   Scheme_Rational *rb = (Scheme_Rational *)b;
00203   Scheme_Object *ma, *mb;
00204 
00205   ma = scheme_bin_mult(ra->num, rb->denom);
00206   mb = scheme_bin_mult(rb->num, ra->denom);
00207 
00208   if (SCHEME_INTP(ma) && SCHEME_INTP(mb)) {
00209     if (or_eq)
00210       return (SCHEME_INT_VAL(ma) <= SCHEME_INT_VAL(mb));
00211     else
00212       return (SCHEME_INT_VAL(ma) < SCHEME_INT_VAL(mb));
00213   } else if (SCHEME_BIGNUMP(ma) && SCHEME_BIGNUMP(mb)) {
00214     if (or_eq)
00215       return scheme_bignum_le(ma, mb);
00216     else
00217       return scheme_bignum_lt(ma, mb);
00218   } else if (SCHEME_BIGNUMP(mb)) {
00219     return SCHEME_BIGPOS(mb);
00220   } else
00221     return !SCHEME_BIGPOS(ma);
00222 }
00223 
00224 int scheme_rational_lt(const Scheme_Object *a, const Scheme_Object *b)
00225 {
00226   return rational_lt(a, b, 0);
00227 }
00228 
00229 int scheme_rational_gt(const Scheme_Object *a, const Scheme_Object *b)
00230 {
00231   return !rational_lt(a, b, 1);
00232 }
00233 
00234 int scheme_rational_le(const Scheme_Object *a, const Scheme_Object *b)
00235 {
00236   return rational_lt(a, b, 1);
00237 }
00238 
00239 int scheme_rational_ge(const Scheme_Object *a, const Scheme_Object *b)
00240 {
00241   return !rational_lt(a, b, 0);
00242 }
00243 
00244 Scheme_Object *scheme_rational_negate(const Scheme_Object *o)
00245 {
00246   Scheme_Rational *r = (Scheme_Rational *)o;
00247 
00248   return make_rational(scheme_bin_minus(scheme_make_integer(0),
00249                                    r->num), 
00250                      r->denom, 0);
00251 }
00252 
00253 Scheme_Object *scheme_rational_add(const Scheme_Object *a, const Scheme_Object *b)
00254 {
00255   Scheme_Rational *ra = (Scheme_Rational *)a;
00256   Scheme_Rational *rb = (Scheme_Rational *)b;
00257   Scheme_Object *ac, *bd, *sum, *cd;
00258   int no_normalize = 0;
00259 
00260   if (SCHEME_INTP(ra->denom) && (SCHEME_INT_VAL(ra->denom) == 1)) {
00261     /* Swap, to take advantage of the next optimization */
00262     Scheme_Rational *rx = ra;
00263     ra = rb;
00264     rb = rx;
00265   }
00266   if (SCHEME_INTP(rb->denom) && (SCHEME_INT_VAL(rb->denom) == 1)) {
00267     /* From Brad Lucier: */
00268     /*    (+ p/q n) = (make-rational (+ p (* n q)) q), no normalize */
00269     ac = ra->num;
00270     cd = ra->denom;
00271     no_normalize = 1;
00272   } else {
00273     ac = scheme_bin_mult(ra->num, rb->denom);
00274     cd = scheme_bin_mult(ra->denom, rb->denom);
00275   }
00276 
00277   bd = scheme_bin_mult(ra->denom, rb->num);
00278   sum = scheme_bin_plus(ac, bd);
00279 
00280   if (no_normalize)
00281     return make_rational(sum, cd, 0);
00282   else
00283     return scheme_make_rational(sum, cd);
00284 }
00285 
00286 Scheme_Object *scheme_rational_subtract(const Scheme_Object *a, const Scheme_Object *b)
00287 {
00288   return scheme_rational_add(a, scheme_rational_negate(b));
00289 }
00290 
00291 Scheme_Object *scheme_rational_add1(const Scheme_Object *n)
00292 {
00293   Small_Rational s;
00294 
00295   return scheme_rational_add(scheme_make_small_rational(1, &s), n);
00296 }
00297 
00298 Scheme_Object *scheme_rational_sub1(const Scheme_Object *n)
00299 {
00300   Small_Rational s;
00301 
00302   return scheme_rational_add(n, scheme_make_small_rational(-1, &s));
00303 }
00304 
00305 Scheme_Object *scheme_rational_multiply(const Scheme_Object *a, const Scheme_Object *b)
00306 {
00307   Scheme_Rational *ra = (Scheme_Rational *)a;
00308   Scheme_Rational *rb = (Scheme_Rational *)b;
00309   Scheme_Object *gcd_ps, *gcd_rq, *p_, *r_, *q_, *s_;
00310 
00311   /* From Brad Lucier: */
00312   /* (* p/q r/s) => (make-rational (* (quotient p (gcd p s))
00313                                       (quotient r (gcd r q)))
00314                                    (* (quotient q (gcd r q))
00315                                       (quotient s (gcd p s)))) */
00316   
00317   gcd_ps = scheme_bin_gcd(ra->num, rb->denom);
00318   gcd_rq = scheme_bin_gcd(rb->num, ra->denom);
00319 
00320   p_ = scheme_bin_quotient(ra->num, gcd_ps);
00321   r_ = scheme_bin_quotient(rb->num, gcd_rq);
00322 
00323   q_ = scheme_bin_quotient(ra->denom, gcd_rq);
00324   s_ = scheme_bin_quotient(rb->denom, gcd_ps);
00325 
00326   p_ = scheme_bin_mult(p_, r_);
00327   q_ = scheme_bin_mult(q_, s_);
00328 
00329   return scheme_make_rational(p_, q_);
00330 }
00331 
00332 Scheme_Object *scheme_rational_max(const Scheme_Object *a, const Scheme_Object *b)
00333 {
00334   int lt;
00335   lt = scheme_rational_lt(a, b);
00336   return scheme_rational_normalize(lt ? b : a);
00337 }
00338 
00339 Scheme_Object *scheme_rational_min(const Scheme_Object *a, const Scheme_Object *b)
00340 {
00341   int lt;
00342   lt = scheme_rational_lt(a, b);
00343   return scheme_rational_normalize(lt ? a : b);
00344 }
00345 
00346 static Scheme_Object *negate_simple(Scheme_Object *v)
00347 {
00348   if (SCHEME_INTP(v))
00349     return scheme_make_integer_value(-SCHEME_INT_VAL(v));
00350   else
00351     return scheme_bignum_negate(v);
00352 }
00353 
00354 Scheme_Object *scheme_rational_divide(const Scheme_Object *n, const Scheme_Object *d)
00355 { 
00356   Scheme_Rational *rd = (Scheme_Rational *)d, *rn = (Scheme_Rational *)n;
00357   Scheme_Rational d_inv;
00358 
00359   /* Check for [negative] inverse, which is easy */
00360   if ((SCHEME_INTP(rn->num) && ((SCHEME_INT_VAL(rn->num) == 1)
00361                             || (SCHEME_INT_VAL(rn->num) == -1)))
00362       && (SCHEME_INTP(rn->denom) && SCHEME_INT_VAL(rn->denom) == 1)) {
00363     int negate = (SCHEME_INT_VAL(rn->num) == -1);
00364     if (SCHEME_INTP(rd->num)) {
00365       if ((SCHEME_INT_VAL(rd->num) == 1)) {
00366        if (negate)
00367          return negate_simple(rd->denom);
00368        else
00369          return rd->denom;
00370       }
00371       if (SCHEME_INT_VAL(rd->num) == -1) {
00372        if (negate)
00373          return rd->denom;
00374        else
00375          return negate_simple(rd->denom);
00376       }
00377     }
00378     if (((SCHEME_INTP(rd->num))
00379         && (SCHEME_INT_VAL(rd->num) < 0))
00380        || (!SCHEME_INTP(rd->num)
00381            && !SCHEME_BIGPOS(rd->num))) {
00382       Scheme_Object *v;
00383       v = negate ? rd->denom : negate_simple(rd->denom);
00384       return make_rational(v, negate_simple(rd->num), 0);
00385     } else {
00386       Scheme_Object *v;
00387       v = negate ? negate_simple(rd->denom) : rd->denom;
00388       return make_rational(v, rd->num, 0);
00389     }
00390   }
00391   
00392   d_inv.so.type = scheme_rational_type;
00393   d_inv.denom = rd->num;
00394   d_inv.num = rd->denom;
00395 
00396   return scheme_rational_multiply(n, (Scheme_Object *)&d_inv);
00397 }
00398 
00399 Scheme_Object *scheme_rational_power(const Scheme_Object *o, const Scheme_Object *p)
00400 {
00401   double b, e, v;
00402 
00403   if (((Scheme_Rational *)p)->denom == one) {
00404     Scheme_Object *a[2], *n;
00405     a[0] = ((Scheme_Rational *)o)->num;
00406     a[1] = ((Scheme_Rational *)p)->num;
00407     n = scheme_expt(2, a);
00408     a[0] = ((Scheme_Rational *)o)->denom;
00409     return make_rational(n, scheme_expt(2, a), 0);
00410   }
00411 
00412   if (scheme_is_rational_positive(o)) {
00413     b = scheme_rational_to_double(o);
00414     e = scheme_rational_to_double(p);
00415 
00416     v = pow(b, e);
00417 
00418 #ifdef USE_SINGLE_FLOATS_AS_DEFAULT
00419     return scheme_make_float(v);
00420 #else
00421     return scheme_make_double(v);
00422 #endif
00423   } else {
00424     return scheme_complex_power(scheme_real_to_complex(o),
00425                             scheme_real_to_complex(p));
00426   }
00427 }
00428 
00429 Scheme_Object *scheme_rational_truncate(const Scheme_Object *o)
00430 {
00431   Scheme_Rational *r = (Scheme_Rational *)o;
00432 
00433   return scheme_bin_quotient(r->num, r->denom);
00434 }
00435 
00436 Scheme_Object *scheme_rational_floor(const Scheme_Object *o)
00437 {
00438   if (scheme_is_rational_positive(o))
00439     return scheme_rational_truncate(o);
00440   else {
00441     Scheme_Object *r;
00442     r = scheme_rational_truncate(o);
00443     return scheme_sub1(1, &r);
00444   }
00445 }
00446 
00447 Scheme_Object *scheme_rational_ceiling(const Scheme_Object *o)
00448 {
00449   if (!scheme_is_rational_positive(o))
00450     return scheme_rational_truncate(o);
00451   else {
00452     Scheme_Object *r;
00453     r = scheme_rational_truncate(o);
00454     return scheme_add1(1, &r);
00455   }  
00456 }
00457 
00458 Scheme_Object *scheme_rational_round(const Scheme_Object *o)
00459 {
00460   Scheme_Rational *r = (Scheme_Rational *)o;
00461   Scheme_Object *q, *qd, *delta, *half;
00462   int more = 0, can_eq_half, negative;
00463 
00464   negative = !scheme_is_rational_positive(o);
00465   
00466   q = scheme_bin_quotient(r->num, r->denom);
00467 
00468   /* Get remainder absolute value: */
00469   qd = scheme_bin_mult(q, r->denom);
00470   if (negative)
00471     delta = scheme_bin_minus(qd, r->num);
00472   else
00473     delta = scheme_bin_minus(r->num, qd);
00474 
00475   half = scheme_bin_quotient(r->denom, scheme_make_integer(2));
00476   can_eq_half = SCHEME_FALSEP(scheme_odd_p(1, &r->denom));
00477 
00478   if (SCHEME_INTP(half) && SCHEME_INTP(delta)) {
00479     if (can_eq_half && (SCHEME_INT_VAL(delta) == SCHEME_INT_VAL(half)))
00480       more = SCHEME_TRUEP(scheme_odd_p(1, &q));
00481     else
00482       more = (SCHEME_INT_VAL(delta) > SCHEME_INT_VAL(half));
00483   } else if (SCHEME_BIGNUMP(delta) && SCHEME_BIGNUMP(half)) {
00484     if (can_eq_half && (scheme_bignum_eq(delta, half)))
00485       more = SCHEME_TRUEP(scheme_odd_p(1, &q));      
00486     else
00487       more = !scheme_bignum_lt(delta, half);
00488   } else
00489     more = SCHEME_BIGNUMP(delta);
00490 
00491   if (more) {
00492     if (negative)
00493       q = scheme_sub1(1, &q);
00494     else
00495       q = scheme_add1(1, &q);      
00496   }
00497 
00498   return q;
00499 }
00500 
00501 
00502 Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o)
00503 {
00504   Scheme_Rational *r = (Scheme_Rational *)o;
00505   Scheme_Object *n, *d;
00506   double v;
00507 
00508   n = scheme_integer_sqrt(r->num);
00509   if (!SCHEME_DBLP(n)) {
00510     d = scheme_integer_sqrt(r->denom);
00511     if (!SCHEME_DBLP(d))
00512       return make_rational(n, d, 0);
00513   }
00514 
00515   v = sqrt(scheme_rational_to_double(o));
00516 
00517 #ifdef USE_SINGLE_FLOATS_AS_DEFAULT
00518   return scheme_make_float(v);
00519 #else
00520   return scheme_make_double(v);
00521 #endif
00522 }
00523 
00524 #define FP_TYPE double
00525 #define SCHEME_RATIONAL_TO_FLOAT scheme_rational_to_double
00526 #define SCHEME_RATIONAL_FROM_FLOAT scheme_rational_from_double
00527 #define SCHEME_BIGNUM_TO_FLOAT_INF_INFO scheme_bignum_to_double_inf_info
00528 #define SCHEME_CHECK_FLOAT scheme_check_double
00529 #define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_double
00530 #define DO_FLOAT_DIV scheme__do_double_div
00531 #define FLOAT_E_MIN -1074
00532 #define FLOAT_M_BITS 52
00533 #include "ratfloat.inc"
00534 
00535 #ifdef MZ_USE_SINGLE_FLOATS
00536 # undef FP_TYPE
00537 # undef SCHEME_RATIONAL_TO_FLOAT
00538 # undef SCHEME_RATIONAL_FROM_FLOAT
00539 # undef SCHEME_BIGNUM_TO_FLOAT_INF_INFO
00540 # undef SCHEME_BIGNUM_FROM_FLOAT
00541 # undef SCHEME_CHECK_FLOAT
00542 # undef DO_FLOAT_DIV 
00543 # undef FLOAT_E_MIN
00544 # undef FLOAT_M_BITS
00545 
00546 #define FP_TYPE float
00547 #define SCHEME_RATIONAL_TO_FLOAT scheme_rational_to_float
00548 #define SCHEME_RATIONAL_FROM_FLOAT scheme_rational_from_float
00549 #define SCHEME_BIGNUM_TO_FLOAT_INF_INFO scheme_bignum_to_float_inf_info
00550 #define SCHEME_CHECK_FLOAT scheme_check_float
00551 #define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_float
00552 #define DO_FLOAT_DIV scheme__do_float_div
00553 #define FLOAT_E_MIN -127
00554 #define FLOAT_M_BITS 23
00555 #include "ratfloat.inc"
00556 #endif
00557