🚀 go-pugleaf

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.