doubles and cos(pi / 2)

[Matthew Fluet]

I think it has to do with the way in which SML/NJ actually implements the
trigonometric functions.  This is from the 110.0.6 version of the
compiler, /boot/math64.sml

    val pi = 3.14159265358979323846
    val negone = ~1.0
    val zero = 0.0
    val half = 0.5
    val one = 1.0
    val two = 2.0
    val four = 4.0

    val maxint = 4503599627370496.0
    fun realround x = if x>=0.0 then x+maxint-maxint else x-maxint+maxint

    local
      val S0 = ~1.6666666666666463126E~1
      val S1 =  8.3333333332992771264E~3
      val S2 = ~1.9841269816180999116E~4
      val S3 =  2.7557309793219876880E~6
      val S4 = ~2.5050225177523807003E~8
      val S5 =  1.5868926979889205164E~10
    in
    fun sin__S z = (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
    end

    local
      val C0 =  4.1666666666666504759E~2
      val C1 = ~1.3888888888865301516E~3
      val C2 =  2.4801587269650015769E~5
      val C3 = ~2.7557304623183959811E~7
      val C4 =  2.0873958177697780076E~9
      val C5 = ~1.1250289076471311557E~11
    in
    fun cos__C z = (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
    end

val PIo4   =  7.853981633974483096E~1
val PIo2   =  1.5707963267948966192E0
val PI3o4  =  2.3561944901923449288E0
val PI     =  pi
val PI2    =  6.2831853071795864769E0
val oneOver2Pi = 0.1591549430918953357688837633725143620345

local
    val thresh =  2.6117239648121182150E~1
in  fun S y = y + y * sin__S(y*y)
    fun C y =
	let val yy = y*y
	    val c = cos__C yy
	    val Y = yy/two
	in  if Y < thresh then one - (Y - c)
	    else half - (Y - half - c)
	end
end
    fun sin x =  (* This function propagages Inf's and Nan's correctly. *)
	let (* x may be finite, inf, or nan at this point. *)
	    val xover2pi = x * oneOver2Pi
	    val x = PI2*(xover2pi - realround(xover2pi))
	       (* now, probably,  ~pi <= x <= pi, except on vaxes.
		  x may be a nan, but cannot now be an inf. *)
	    fun lessPIo2 x = if x>PIo4 then C(PIo2-x) else S x
	    fun lessPI x = if x>PIo2 then lessPIo2(PI-x) else lessPIo2 x
	    fun positive x = if x>PI then sin(x-PI2) (* exceedingly rare *)
			             else lessPI x
	 in if x>=0.0 
		then positive x
	        else ~(positive(~x))
	end

    fun cos x = sin(PIo2-x)

    fun tan x = sin x / cos x


So, your basic expansion version of sin, with trig identities to cos and
tan.  Plug all this into MLton and compute cos(Math.pi / 2.0) and you do
indeed get 0.0.

They do a similar expansion for exponential and logarithmic functions.  Of
course, they now can't make use of any special processor instructions that
might compute these values.

This of course doesn't answer the question about whether or not it is IEEE
correct.  And if you make use of the identity  cos x = sin(pi/2 - x)  then
MLton does return 0.0 for cos(pi/2).  

On Wed, 4 Oct 2000, Stephen Weeks wrote:

> 
> > I don't understand
> > how the OCAML entries got this test right, 
> ...
> > OCAML		6.12303176911e-17
> 
> I should have mentioned that this was OCAML 2.02.  I just downloaded 3.0 and
> indeed, it prints the result as zero.  So now OCAML and SML/NJ agree the result
> is zero, while everyone else thinks the result is 6.123E-17.
>