Back to index

plt-scheme  4.2.1
ratfloat.inc
Go to the documentation of this file.
00001 
00002 /* Optimization sometimes causes a problem: n or d is represented in an
00003    extended format instead of a `double'. We don't want to turn off
00004    floatng-point optimizations in the rest of the program, so we use a
00005    little function to defeat the optimization: */
00006 
00007 FP_TYPE DO_FLOAT_DIV(FP_TYPE n, FP_TYPE d)
00008 {
00009   return n / d;
00010 }
00011 
00012 
00013 FP_TYPE SCHEME_RATIONAL_TO_FLOAT(const Scheme_Object *o)
00014 {
00015   Scheme_Rational *r = (Scheme_Rational *)o;
00016   FP_TYPE n, d;
00017   int ns, ds;
00018 
00019   if (SCHEME_INTP(r->num)) {
00020     n = (FP_TYPE)SCHEME_INT_VAL(r->num);
00021     ns = 0;
00022   } else
00023     n = SCHEME_BIGNUM_TO_FLOAT_INF_INFO(r->num, 0, &ns);
00024 
00025   if (SCHEME_INTP(r->denom)) {
00026     d = (FP_TYPE)SCHEME_INT_VAL(r->denom);
00027     ds = 0;
00028   } else
00029     d = SCHEME_BIGNUM_TO_FLOAT_INF_INFO(r->denom, 0, &ds);
00030 
00031   if (ns || ds) {
00032     /* Quick path doesn't necessarily work. The more general
00033        way is adpated from Gambit-C 4.1. */
00034     long nl, dl, p, shift;
00035     Scheme_Object *a[2], *n, *d;
00036     FP_TYPE res;
00037 
00038     a[0] = r->num;
00039     n = scheme_abs(1, a);
00040 
00041     d = r->denom;
00042 
00043     nl = scheme_integer_length(n);
00044     dl = scheme_integer_length(d);
00045 
00046     p = nl - dl;
00047     if (p < 0) {
00048       a[0] = n;
00049       a[1] = scheme_make_integer(-p);
00050       
00051       n = scheme_bitwise_shift(2, a);
00052     } else {
00053       a[0] = d;
00054       a[1] = scheme_make_integer(p);
00055       
00056       d = scheme_bitwise_shift(2, a);
00057     }
00058 
00059     if (scheme_bin_lt(n, d)) {
00060       a[0] = n;
00061       a[1] = scheme_make_integer(1);
00062       
00063       n = scheme_bitwise_shift(2, a);
00064       --p;
00065     }
00066 
00067     shift = p - FLOAT_E_MIN;
00068     if (shift > FLOAT_M_BITS) {
00069       shift = FLOAT_M_BITS;
00070     }
00071     
00072     a[0] = n;
00073     a[1] = scheme_make_integer(shift);
00074     n = scheme_bitwise_shift(2, a);
00075 
00076     n = scheme_bin_div(n, d);
00077     if (SCHEME_RATIONALP(n))
00078       n = scheme_rational_round(n);
00079     
00080     if (SCHEME_INTP(n))
00081       res = (FP_TYPE)SCHEME_INT_VAL(n);
00082     else
00083       res = SCHEME_BIGNUM_TO_FLOAT_INF_INFO(n, 0, NULL);
00084 
00085     res = res * pow(2, p - shift);
00086 
00087     if (SCHEME_INTP(r->num)) {
00088       if (SCHEME_INT_VAL(r->num) < 0)
00089         res = -res;
00090     } else if (!SCHEME_BIGPOS(r->num)) {
00091       res = -res;
00092     }
00093 
00094     return res;
00095   }
00096 
00097   return DO_FLOAT_DIV(n, d);
00098 }
00099 
00100 Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d)
00101 {
00102   double frac, i;
00103   int count, exponent, is_neg;
00104   Scheme_Object *int_part, *frac_part, *frac_num, *frac_denom, *two, *result;
00105 #ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS
00106   int negate;
00107   if (d < 0) {
00108     d = -d;
00109     negate = 1;
00110   } else
00111     negate = 0;
00112 #endif
00113 
00114   SCHEME_CHECK_FLOAT("inexact->exact", d, "exact");
00115 
00116   is_neg = (d < 0);
00117 
00118   frac = modf((double)d, &i);
00119   (void)frexp(d, &exponent);
00120 
00121   int_part = SCHEME_BIGNUM_FROM_FLOAT(i);
00122 
00123   if (!frac) {
00124 #ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS
00125    if (negate)
00126      return scheme_bin_minus(scheme_make_integer(0), int_part);
00127 #endif  
00128     return int_part;
00129   }
00130 
00131   frac_num = scheme_make_integer(0);
00132   frac_denom = one;
00133   two = scheme_make_integer(2);
00134 
00135   count = 0;
00136   while (frac) {
00137     count++;
00138     frac_num = scheme_bin_mult(frac_num, two);
00139     frac_denom = scheme_bin_mult(frac_denom, two);    
00140     frac = modf(ldexp(frac, 1), &i);
00141     if (i) {
00142       if (is_neg)    
00143        frac_num = scheme_bin_minus(frac_num, one);
00144       else
00145        frac_num = scheme_bin_plus(frac_num, one);
00146     }
00147   }
00148 
00149   frac_part = scheme_bin_div(frac_num, frac_denom);
00150 
00151   result = scheme_bin_plus(int_part, frac_part);
00152 
00153 #ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS
00154   if (negate)
00155     return scheme_bin_minus(scheme_make_integer(0), result);
00156 #endif
00157 
00158   return result;
00159 }