Le Jeudi 29 Mai 2003 00:29, Frank Heckenbach a écrit :
Boris Herman wrote:
If there's interest, I could add ArcSin and ArcCos to the runtime library.
This would be great!
I've found out that I get more precision by using this self-written function:
function ArcSin(x:double):double; const c1=1e-6; begin if abs(x-1)<c1 then ArcSin:=pi/2 else if abs(1+x)<c1 then ArcSin:=-pi/2 else ArcSin:=ArcTan(x/sqrt(1-sqr(x))); end;
I hope Emil will comment on the numeric stability of this function etc. I'm a bit worried about the c1 definition. It might be suitable for the Double type on your platform, but I think a general function should try to work with MinReal, EpsReal etc.
This code gives only something like simple precision on my machine (AMD athlon, that is IEEE floating point arithmetic)
It would probably be much safer to rely upon the libm library, which is (hopefully) optimized at assembly level for most hardwares. Assembly is the only way to be both accurate and efficient.
If one wants to reinvent the wheel using only standard instructions in Pascal, here is a solutions accurate to 1E-15 if the hardware has true IEEE double precision. Near 1, i use the Taylor expansion of arcsin(1-v^2) Please note that this is not optimized, I invent a wheel much more like that of a merovingian ox-cart than that of a Ferrari. **************************************************************************** function arcsin(x : double):double; var u,v,w : double; begin if (abs(x) < 0.99995) then arcsin:=arctan(x/sqrt(1-sqr(x))) else begin u:=1-abs(x); v:=sqrt(u); w:=1.5707963267948966192 -v*((u*0.26516504294495532165e-1 +0.11785113019775792073)*u +1.4142135623730950488); if (x>0) then arcsin:=w else arcsin := -w; end; end; { arcsin } ***************************************************************************** One should trigger an exception if abs(x) > 1 This is not included here.