Re: ratio implementation

[Kevin Ryde <user42@zip.com.au>]

[-- Attachment #1: Type: text/plain, Size: 819 bytes --]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>
>   (sqrt 4) => 2 (exact)
>
> I'm not sure whether exact square roots are imortant: they will only
> be useful when both the numerator and denominator are squares of
> integers and those pairs will be rare, I'd say.

I'd wondered if some sort of isqrt or isqrt+remainder would be of more
value to those wanting integer roots.  (Though I'm aware r5rs says
exact results for exact perfect squares is desirable.)

> We could use GMP for
> computing bignum square roots...

Rob mentioned this to me at one stage, I actually started on an
attempt at it.  Code below, which might sort of work, maybe.

One concern though is to ensure the root calculated is the same for a
"double" input (which uses sqrt()) as compared to some bignum root
(using gmp + conversions).


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: sqrt.c --]
[-- Type: text/x-csrc, Size: 3815 bytes --]


/* sqrt and mpz_sqrt+mpz_get_d don't round the same

   */


#include <float.h>


#ifndef DBL_MANT_DIG
#define DBL_MANT_DIG  53  /* by default assume IEEE double */
#endif


SCM scm_sqrt (SCM x);
SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
            (SCM x),
	    "Return the square root of @var{x}.\n"
	    "\n"
	    "If @var{x} is exact and a perfect square then the return is\n"
	    "exact, otherwise the return is inexact (either real or\n"
	    "complex).")
#define FUNC_NAME s_scm_sqrt
{
  double root_d;
  SCM    root;
  int    neg;

  if (SCM_INUMP (x))
    {
      long x_l = SCM_INUM (x);
      neg = (x_l < 0);
      x_l = (neg ? -x_l : x_l);

      /* inums which fit the mantissa of a double are handled with sqrt().
         On 32-bit systems inums always fit, but on 64-bit systems we need
         to check.  */
#if SCM_I_FIXNUM_BIT > DBL_MANT_DIG
      if ((unsigned long) x_l < (1UL << DBL_MANT_DIG))
#endif
        {
          root_d = sqrt ((double) x_l);
          if (! neg && floor (root_d) == root_d)
            {
              /* perfect square, return as inum */
              long root_l = (long) root_d;
              return SCM_MAKINUM (neg ? -root_l : root_l);
            }
          goto inexact;
        }

      /* inums which don't fit the mantissa of a double are handled with the
         bignum code, to ensure full accuracy in the result  */
      root = scm_i_mkbig ();
      mpz_set_si (SCM_I_BIG_MPZ (root), x_l);
      x = root;
      goto big;
    }
  else if (SCM_BIGP (x))
    {
      int  scale;

      root = scm_i_mkbig ();
    big:
      if (mpz_perfect_square_p (SCM_I_BIG_MPZ (x)))
        {
          /* Perfect square, return as exact.
             Enhancement: mpz_perfect_square_p already calculates the root
             when proving squareness, it'd be nice to have some sort of
             "square root if perfect square" function.  */
          mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
          scm_remember_upto_here_1 (x);
          return scm_i_normbig (root);
        }

      neg = (mpz_sgn (SCM_I_BIG_MPZ (x)) < 0);
      mpz_abs (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
      scm_remember_upto_here_1 (x);

      /* We need DBL_MANT_DIG bits of root, which means twice that of
         operand.  If we've got less, then pad with zeros, effectively
         generating fractional bits for the root.  If we've got more, then
         truncate so as to save work in mpz_sqrt.  */
      scale = ((long) mpz_sizeinbase (SCM_I_BIG_MPZ (root), 2)
               - 1 - 2*DBL_MANT_DIG) / 2;
      if (scale > 0)
        mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                         2 * scale);
      else
        mpz_mul_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                      - 2 * scale);

      mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root));
      root_d = ldexp (mpz_get_d (SCM_I_BIG_MPZ (root)), scale);
      scm_remember_upto_here_1 (root);
      goto inexact;
    }
  else if (SCM_REALP (x))
    {
      double x_d = SCM_REAL_VALUE (x);
      neg = (x_d < 0.0);
      root_d = sqrt (fabs (x_d));

    inexact:
      if (neg)
        return scm_make_complex (0.0, root_d);
      else
        return scm_make_real (root_d);
    }
  else if (SCM_COMPLEXP (x))
    {
      /* sqrt(x) = (make-polar (sqrt (magnitude x)) (/ (angle x) 2))
         ENHANCE-ME: Use csqrt() when available.  */
      double re = SCM_COMPLEX_REAL (x);
      double im = SCM_COMPLEX_IMAG (x);
      double mag = sqrt (hypot (re, im));
      double ang = atan2 (re, im) / 2.0;
      double s, c;
#if HAVE_SINCOS
      sincos (ang, &s, &c);
#else
      s = sin (ang);
      c = cos (ang);
#endif
      return scm_make_complex (mag * c, mag * s);
    }
  else
    SCM_WRONG_TYPE_ARG (SCM_ARG1, x);
}
#undef FUNC_NAME

[-- Attachment #3: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel