Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!swrinde!elroy.jpl.nasa.gov!ucla-cs!ucla-ma!euphemia!pmontgom From: pmontgom@euphemia.math.ucla.edu (Peter Montgomery) Newsgroups: comp.arch Subject: Re: float/float = integer, remainder Message-ID: <1991May17.090214.7830@math.ucla.edu> Date: 17 May 91 09:02:14 GMT References: <9105150230.AA26299@ucbvax.Berkeley.EDU> Sender: news@math.ucla.edu Organization: UCLA Mathematics Dept. Lines: 81 In article <9105150230.AA26299@ucbvax.Berkeley.EDU> jbs@WATSON.IBM.COM writes: > Peter Montgomery gives an example (computing sin(1234.) in >4 digit decimal) which purports to show how to use a float/float -> >integer, remainder instruction for argument reduction. There are some >problems with his method (as I will discuss later). In any case it is >besides the point. I am not claiming you can not use this instruction >for argument reduction. I am claiming it will not enable you to do >argument reduction any faster. Then I misunderstood you. I believe that the remainder is needed seldom enough that it may be reasonable to require up to (say) five floating point instructions to get it. But it is essential that one be able to request this operation in a portable way which compilers can translate into efficient object code producing an accurate result, such as Q = round(Y/X); (take quotient, round to nearest integer, keep in floating format) R = subtract_product(Y, Q, X); (Y - Q*X without intermediate rounding) Writing R = Y - X*round(Y/X) (say) is unacceptable if IEEE rounding is done after the multiply, as it will not produce the remainder which the application wants. The point of my previous posting was that the remainder Herman Rubin requests can be useful. Now the question is how one can obtain that remainder on current and proposed systems. I asked "How can Herman Rubin get this remainder accurately with the hardware and software commonly available?" The IBM RS/6000 has such an instruction, but how does one write code that works there as well as the SPARC series, the MIPS R3000 series, and elsewhere? Shearer goes on to mention how using ln(2) = z1 + z2 where z1 has only 42 significant bits we can compute exponentials since the quotient will be at most 11 bits. It appears sound. > For trigonometric functions we let n be the closest integer to >x/(pi/2) and y=x-n*(pi/2). The computation of n and y is trickier for >reasons: > 1) n may be very large. > 2) y must be computed with small relative error not just small >absolute error (since sin has a zero at zero). > Peter Montgomery's method does not deal with these problems. >Since the quotient is returned as an integer it will overflow whenever >n>2**31. His example assumes his integers have as many digits as his >floats. The intermediate quotient can be indeed stored as a float to guard against this very problem, using a separate instruction to convert it to integer if the application knows it will be in range. >His example also assumes 8 digits of pi/2 suffice for accurate argument >argument reduction of 4 digit numbers. This is not true as can be seen >by trying sin(355.). Yes, if one wants four significant figures of sin(355) ~ -.00003014. But for many applications, the absolute error in trigonometric functions is more important than the relative error; for example, it matters little whether one approximates the complex number 10000*(cos(355), sin(355)) as (-10000, 0) or as (-10000, -0.3014). On a SUN 4, if I compute sin(2549491779) using double precision arithmetic, I get 4.507E-10 whereas Maple gives 4.474E-10. Despite the huge quotient, the result is accurate to 11 decimal places, although to only 2 significant digits; however it is the decimal places which is important for trig functions. If this error is unacceptable, one needs a more accurate approximation to pi/2, as Shearer suggests. And since the subsequent corrections won't have the direct form of a remainder from divide, putting the remainder directly in the floating point divide instruction is probably not the best solution. But the functionality is needed somehow. -- Peter L. Montgomery pmontgom@MATH.UCLA.EDU Department of Mathematics, UCLA, Los Angeles, CA 90024-1555 If I spent as much time on my dissertation as I do reading news, I'd graduate.