2002-03-05 John W. Eaton * numbers.c (scm_divide): Adapt code from libstdc++/f2c to void potential overflow problems. Index: libguile/numbers.c =================================================================== RCS file: /cvs/guile/guile-core/libguile/numbers.c,v retrieving revision 1.157 diff -u -r1.157 numbers.c --- libguile/numbers.c 1 Mar 2002 00:19:20 -0000 1.157 +++ libguile/numbers.c 6 Mar 2002 03:46:53 -0000 @@ -3691,6 +3691,33 @@ #undef FUNC_NAME +/* The code below for complex division is adapted from the GNU + libstdc++, which adapted it from f2c's libF77, and is subject to + this copyright: */ + +/**************************************************************** +Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore. + +Permission to use, copy, modify, and distribute this software +and its documentation for any purpose and without fee is hereby +granted, provided that the above copyright notice appear in all +copies and that both that the copyright notice and this +permission notice and warranty disclaimer appear in supporting +documentation, and that the names of AT&T Bell Laboratories or +Bellcore or any of their entities not be used in advertising or +publicity pertaining to distribution of the software without +specific, written prior permission. + +AT&T and Bellcore disclaim all warranties with regard to this +software, including all implied warranties of merchantability +and fitness. In no event shall AT&T or Bellcore be liable for +any special, indirect or consequential damages or any damages +whatsoever resulting from loss of use, data or profits, whether +in an action of contract, negligence or other tortious action, +arising out of or in connection with the use or performance of +this software. +****************************************************************/ + SCM_GPROC1 (s_divide, "/", scm_tc7_asubr, scm_divide, g_divide); /* Divide the first argument by the product of the remaining arguments. If called with one argument @var{z1}, 1/@var{z1} is @@ -3724,8 +3751,15 @@ } else if (SCM_COMPLEXP (x)) { double r = SCM_COMPLEX_REAL (x); double i = SCM_COMPLEX_IMAG (x); - double d = r * r + i * i; - return scm_make_complex (r / d, -i / d); + if (r <= i) { + double t = r / i; + double d = i * (1.0 + t * t); + return scm_make_complex (t / d, -1.0 / d); + } else { + double t = i / r; + double d = r * (1.0 + t * t); + return scm_make_complex (1.0 / d, -t / d); + } } else { SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide); } @@ -3765,8 +3799,15 @@ { double r = SCM_COMPLEX_REAL (y); double i = SCM_COMPLEX_IMAG (y); - double d = r * r + i * i; - return scm_make_complex ((a * r) / d, (-a * i) / d); + if (r <= i) { + double t = r / i; + double d = i * (1.0 + t * t); + return scm_make_complex ((a * t) / d, -a / d); + } else { + double t = i / r; + double d = r * (1.0 + t * t); + return scm_make_complex (a / d, -(a * t) / d); + } } } else { SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); @@ -3870,9 +3911,15 @@ } else if (SCM_COMPLEXP (y)) { double ry = SCM_COMPLEX_REAL (y); double iy = SCM_COMPLEX_IMAG (y); - double d = ry * ry + iy * iy; - return scm_make_complex ((rx * ry + ix * iy) / d, - (ix * ry - rx * iy) / d); + if (ry <= iy) { + double t = ry / iy; + double d = iy * (1.0 + t * t); + return scm_make_complex ((rx * t + ix) / d, (ix * t - rx) / d); + } else { + double t = iy / ry; + double d = ry * (1.0 + t * t); + return scm_make_complex ((rx + ix * t) / d, (ix - rx * t) / d); + } } else { SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); }