(* ONE of JRS, Noe, PZa, or TjM must be defined for this to compile; 
   they select the method.  PZa and TjM require {$N+} to compile,
   and an adequate FPU to run. *)


program INVCOS { jrs@merlyn.demon.co.uk 96/05/18 } ;

(*
  <Quote taken from Jon Shemitz's "c.l.p.b FAQ 0.5" of 1995>

  Borland Pascal does not have ArcSin or ArcCos functions.
  It does have an ArcTan function, and the online help for
  that function gives the following conversion formulae:

  ArcSin(x) = ArcTan (x/sqrt (1-sqr (x)))
  ArcCos(x) = ArcTan (sqrt (1-sqr (x)) /x)

  Dr. John Stockton, jrs@merlin.dclf.npl.co.uk, points out
  that ArcTan will always return a value between -Pi / 2
  and Pi / 2.  Also, there are two angles in the range from
  0 to 2 Pi for any given sine or cosine value, even though
  the formulae above will only give you one of them. For this
  reason, the formula   ArcCos(x) = (Pi / 2) - ArcSin(X)
  from the CRC Standard Mathematical Tables handbook is
  generally better.

  <end quote>


  This still does not cover one of the two points; Borland's
  ArcCos is actually ***WRONG*** for negative arguments, in
  that cos(arccos(x))<>x, as the following code shows. ArcSin
  itself as given is inoperable for Abs(X)=1; for X such that
  Abs(X)>sqrt(0.5) it may be better to get ArcSin from
  ArcTan(reciprocal).


  DRAFT REPLACEMENT by JRS@dclf.npl.co.uk :

  Borland Pascal does not have ArcSin or ArcCos functions.
  It does have an ArcTan function, and the online help for
  that function gives the following conversion formulae:

  o        ArcSin(x) = ArcTan (x/sqrt (1-sqr (x)))
  o        ArcCos(x) = ArcTan (sqrt (1-sqr (x)) /x)

  These are unsatisfactory for program statements, although
  trigonometrically correct -
  1: ArcSin will fail, divide-by-zero, if Abs(x)=1.0; likewise
  for ArcCos if x=0.0; these correspond to angles +/- 0.5*Pi.
  2: ArcTan will always return a value _between_ -0.5*Pi and
  +0.5*Pi exclusive, which is nearly OK for ArcSin; but for
  ArcCos the range must be 0 to Pi inclusive (in fact one
  need only add Pi when x<0.0, getting a different angle
  with the same tangent).  With a corrected ArcSin, the better
  formula ArcCos(x) = 0.5*Pi - ArcSin(x) [from the CRC Standard
  Mathematical Tables handbook] always gives the Principal Value.

  NOTE also, there are two angles in the range from 0 to 2*Pi for
  any given sine or cosine (or tangent) value, either of which
  might be appropriate; but the functions ArcSin2, ArcCos2 below
  will only give the Principal Value.

  But ATAN2(Y, X), as provided by a 486's FPU's FPATAN, is more
  useful than ArcTan(X) itself [& ArcTan(x)=ATAN2(1, x)]; it
  gives the best way of doing ArcSin, ArcCos.

  END DRAFT

  Prof. Seppo Pynnonen, sjp@uwasa.fi, knows about these things;
  also Norbert Juffa?.  Someone clever could produce a nice
  little MOREMATH.PAS unit requiring {$N+}{$E-} and 486, using
  the better FPU op-codes and asm...end stuff to do all trig.

  Later : Now alternative ArcCos from Tony Noe ...

  Later : Paul Zank of Lockheed Sanders has a method using ATAN2
  which is in all PC FPUs - it looks good, except in {$N-} ...

  Later : Taken Terje Mathisen's idea from news:*.asm.x86 ...

  *)


{$IFDEF WINDOWS} uses WinCrt ; {$ENDIF}

type float = {$IFOPT N+} extended {$ELSE} real {$ENDIF} ;

{as BP7 help}
function ArcSin(x : float) : float { FAILS if Abs(x)=1.0 } ;
begin ArcSin := ArcTan(x/Sqrt(1-Sqr(x))) end ;

{as BP7 help}
function ArcCos(x : float) : float { FAILS if x=0.0; WRONG if x<0.0 } ;
begin
  ArcCos := ArcTan(Sqrt(1-Sqr(x))/x) { +Pi*Ord(x<0) will correct it }
end ;


const ID = { Define one of these }
  {$IFDEF JRS} 'JRS' {$ENDIF}
  {$IFDEF Noe} 'Noe' {$ENDIF}
  {$IFDEF PZa} 'PZa' {$ENDIF}
  {$IFDEF TjM} 'TjM' {$ENDIF}
  ;

{$IFDEF TjM}
(* Adapted
  From: Terje Mathisen <Terje.Mathisen@hda.hydro.com>
  Date: Tue, 30 Apr 1996 15:33:26 +0200
  o        tan(y) = sin(y) / cos(y)
  o        tan(y) = sqrt(1-cos(y)^2) / cos(y)
  o  so,   arccos(x) = FPATAN(sqrt(1-x^2),x)
  *)

function ArcCos2(x : extended) : extended ;
assembler ; asm
    fld [x]           { X }
    fld st(0)         { X   X }
    fmul st,st(0)     { X^2 X }
    fld1              { 1 X^2 X }
    fsubrp st(1),st   { 1-X^2 X }
    FSQRT             { 1-X^2 X }
    fxch st(1)        { X 1-X^2 }
    fpatan            { Uses both arguments }
    end {ArcCos2} ;

function ArcSin2(x : extended) : extended ;
assembler ; asm
    fld [x]           { X }
    fld st(0)         { X   X }
    fmul st,st(0)     { X^2 X }
    fld1              { 1 X^2 X }
    fsubrp st(1),st   { 1-X^2 X }
    FSQRT             { 1-X^2 X }
    fpatan            { Uses both arguments }
    end {ArcSin2} ;
(* Terje.Mathisen@hda.hydro.com *)
{$ENDIF TjM}


{$IFDEF PZa} { From: Paul Zank <PZANK@mailgw.sanders.lockheed.com> }
function ArcTan2(y : extended ; x : extended) : extended ;
{ The '87, '287, '387, i486, and i586 have an op code to do ATAN2!
  It is advertised to be well behaved and execute in 290 cycles on
  a i486.  This would be about 3 uSec on a 100MHz i486 machine!  It
  is also supposed to be very accurate.  More information can be
  found in the "i486 Microprocessor - Programmers Reference Manual" }
assembler ;
asm  fld [y] ; fld [x] ; fpatan  end ;

function ArcSin2(x : extended) : extended;
begin
  ArcSin2 :=  ArcTan2(x, Sqrt(1-Sqr(x))) ;
  { ArcTan2 is used instead of ArcTan 
    to avoid exception when x = 1 or -1 }
  end ;

function ArcCos2(x : extended) : extended ;
begin
  ArcCos2 := ArcTan2(Sqrt(1-Sqr (x)), x) ;
  { ArcTan2 is used instead of ArcTan to avoid exception when x = 0}
  end ;
{$ENDIF PZa}


{$IFDEF Noe} { Idea from: Tony Noe <noe@sspectra.com> }
function ArcCos2(x : float) : float ;
begin if X=-1.0 then ArcCos2 := Pi {jrs} else
    ArcCos2 := 2.0 * ArcTan(Sqrt((1.0 - x) / (1.0 + x))) {?} end ;

function ArcSin2(x : float) : float ;
begin ArcSin2 := 0.5*Pi-ArcCos2(x) end ;
{$ENDIF}


{$IFDEF JRS}
function ArcSin2(x : float) : float ;
var T : float ;
begin T := Sqrt(1-Sqr(x)) ;
  if T>=0.7071 {about 1/Root2} then ArcSin2 := ArcTan(x/T) else
    begin T := ArcTan(T/x) ;
    if x>=0.0 then ArcSin2 := 0.5*Pi-T else ArcSin2 := -0.5*Pi-T end ;
  end {OK} ;

function ArcCos2(x : float) : float ;
begin ArcCos2 := 0.5*Pi-ArcSin2(x) {OK} end ;
{$ENDIF}


var J, K : integer ; X, Y, S, T : float ;
const C : array [boolean] of char =
  {$IFDEF WINDOWS} 'x ' {$ELSE} ' '#251 {$ENDIF} ;


BEGIN ;
{$IFDEF WINDOWS}
with WindowSize do begin X := 660 ; Y := 480 end ;
{$ENDIF}
Writeln('INVCOS - see program for ', ID, ' algorithm.'^M^J, 'X':6,
  'Cos(ArcCos(X))':17,
  'ArcSin2(X)':14,
  'ArcCos2(X)':14,
  'Cos(ArcCos2(X))':17, '  '#127'*10^',
  {$IFOPT N+} '21' {$ELSE} '13' {$ENDIF} ) ;
S := 0.0 ; K := 0 ;
for J := -10 to 10 do begin X := 0.1*J ; Write(X:6:3) ;
  if J=0 then Write('Error 200 !':17) {ArcCos(0) divides by 0}
    else Write(Cos(ArcCos(X)):17:6) ;
  Write(ArcSin2(X):14:6) ;
  Write(ArcCos2(X):14:6) ;
  Y := Cos(ArcCos2(X)) ; Write(Y:17:6, C[X=Y]:3) ;
  T := (X-Y) * {$IFOPT N+} 1E21 {$ELSE} 1E13 {$ENDIF} ;
  Write(T:6:0) ; S := S + Sqr(T) ; Inc(K) ;
  Writeln end ;
Writeln('':60, ID, ' : RMS =', Sqrt(S/K):6:0) ;
{$IFNDEF WINDOWS} Write('<cr>') ; Readln ; {$ENDIF}
END. {
    JRS@dclf.npl.co.uk     JRS@merlyn.demon.co.uk
}
