From mboxrd@z Thu Jan 1 00:00:00 1970 Path: main.gmane.org!not-for-mail From: Kevin Ryde Newsgroups: gmane.lisp.guile.devel Subject: Re: ratio implementation Date: Wed, 15 Oct 2003 08:56:05 +1000 Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Message-ID: <87r81fzawq.fsf@zip.com.au> References: <3F250809.9030108@ccrma.stanford.edu> <87smmyibk7.fsf@zagadka.ping.de> <3F6637EC.7010004@dirk-herrmanns-seiten.de> <3F66F68B.3070100@ccrma.stanford.edu> <3F6A1F1A.8000507@dirk-herrmanns-seiten.de> <87pth9cbmt.fsf@zagadka.ping.de> <3F8A853D.1020708@ccrma> NNTP-Posting-Host: deer.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: sea.gmane.org 1066173244 16232 80.91.224.253 (14 Oct 2003 23:14:04 GMT) X-Complaints-To: usenet@sea.gmane.org NNTP-Posting-Date: Tue, 14 Oct 2003 23:14:04 +0000 (UTC) Cc: Bill Schottstaedt , guile-devel@gnu.org Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Wed Oct 15 01:14:02 2003 Return-path: Original-Received: from monty-python.gnu.org ([199.232.76.173]) by deer.gmane.org with esmtp (Exim 3.35 #1 (Debian)) id 1A9YMb-0007gx-00 for ; Wed, 15 Oct 2003 01:14:01 +0200 Original-Received: from localhost ([127.0.0.1] helo=monty-python.gnu.org) by monty-python.gnu.org with esmtp (Exim 4.24) id 1A9YLt-0000CH-Ic for guile-devel@m.gmane.org; Tue, 14 Oct 2003 19:13:17 -0400 Original-Received: from list by monty-python.gnu.org with tmda-scanned (Exim 4.24) id 1A9YJR-0007Rm-SY for guile-devel@gnu.org; Tue, 14 Oct 2003 19:10:45 -0400 Original-Received: from mail by monty-python.gnu.org with spam-scanned (Exim 4.24) id 1A9YE7-0003en-QS for guile-devel@gnu.org; Tue, 14 Oct 2003 19:05:47 -0400 Original-Received: from [199.232.41.8] (helo=mx20.gnu.org) by monty-python.gnu.org with esmtp (TLSv1:DES-CBC3-SHA:168) (Exim 4.24) id 1A9Y5e-0000sp-CM for guile-devel@gnu.org; Tue, 14 Oct 2003 18:56:30 -0400 Original-Received: from [61.8.0.36] (helo=snoopy.pacific.net.au) by mx20.gnu.org with esmtp (Exim 4.24) id 1A9Y5X-0007Vc-WA for guile-devel@gnu.org; Tue, 14 Oct 2003 18:56:24 -0400 Original-Received: from mongrel.pacific.net.au (mongrel.pacific.net.au [61.8.0.107]) by snoopy.pacific.net.au (8.12.3/8.12.3/Debian-6.6) with ESMTP id h9EMuFYH021901; Wed, 15 Oct 2003 08:56:15 +1000 Original-Received: from localhost (ppp113.dyn228.pacific.net.au [203.143.228.113]) by mongrel.pacific.net.au (8.12.3/8.12.3/Debian-6.6) with ESMTP id h9EMqKst005947; Wed, 15 Oct 2003 08:52:22 +1000 Original-Received: from gg by localhost with local (Exim 3.35 #1 (Debian)) id 1A9Y5F-0000gK-00; Wed, 15 Oct 2003 08:56:05 +1000 Original-To: Marius Vollmer Mail-Copies-To: never In-Reply-To: (Marius Vollmer's message of "Tue, 14 Oct 2003 14:39:51 +0200") User-Agent: Gnus/5.1003 (Gnus v5.10.3) Emacs/21.3 (gnu/linux) X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.2 Precedence: list List-Id: Developers list for Guile, the GNU extensibility library List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: main.gmane.org gmane.lisp.guile.devel:2874 X-Report-Spam: http://spam.gmane.org/gmane.lisp.guile.devel:2874 --=-=-= Marius Vollmer 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). --=-=-= Content-Type: text/x-csrc Content-Disposition: attachment; filename=sqrt.c /* sqrt and mpz_sqrt+mpz_get_d don't round the same */ #include #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 --=-=-= Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ Guile-devel mailing list Guile-devel@gnu.org http://mail.gnu.org/mailman/listinfo/guile-devel --=-=-=--