(* ONE of JRS, Noe, PZa, TjM, JRF, TS1, or TS3 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 { www.merlyn.demon.co.uk 96/12/14 } ;

(*
  <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, www.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 www.merlyn.demon.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 JRF} 'JRF' {$ENDIF}
  {$IFDEF TS1} 'TS1'  {$ENDIF}
  {$IFDEF TS3} 'TS3'  {$ENDIF}
  ;

(** N.B. What if, due to error, the argument is greater in magnitude
  than 1.0 - (a) infinitesimally, (b) significantly ? **)

{$IFDEF TS1}
{ Adapted from presumed draft TSFAQP #109 of 96/11/23 }
function ArcSin2(x : float) : float ;
const tiny = 1.0E-37 ;
begin
  if abs(x-1.0)<tiny then ArcSin2 := Pi/2.0
    else if abs(x+1.0) < tiny then ArcSin2 := -Pi/2.0
    else ArcSin2 := ArcTan(x/Sqrt(1.0-x*x)) ;
  end {ArcSin2} ;

function ArcCos2(x : float) : float ;
const tiny = 1.0E-37 ;
begin
  if abs(x)<tiny then ArcCos2 := Pi/2
    else ArcCos2 := ArcTan(Sqrt(1.0-x*x)/x) + Pi*Ord(x < 0.0) ;
  end {ArcCos2} ;
{$ENDIF TS1}

{$IFDEF TS3}
{ Adapted from presumed draft TSFAQP #109 of 96/11/26 }
function ArcSin2(x : float) : float ;
const { small = 1.0E-37 ; }
HalfPi = Pi/2.0 ;
begin
  if (x < -1.0) or (x > +1.0) then begin
    Writeln ('ArcSin2 argument ', x, ' out of range [-1,+1]') ;
    HALT end ;
  if x = +1.0 {- small} then ArcSin2 := +HalfPi
    else if x = -1.0 {+ small} then ArcSin2 := -HalfPi
    else ArcSin2 := ArcTan(x/Sqrt(1.0-Sqr(x))) ;
  end {ArcSin2} ;

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

{$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             { sqrt(1-X^2) X }
    fxch st(1)        { X sqrt(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             { sqrt(1-X^2) X }
    fpatan            { Uses both arguments }
    end {ArcSin2} ;
(*
  It is possible that the FSUB instruction should
  have been FSUBP instead of FSUBRP, and you might
  need to interpose an FXCH ST(1) between the last
  two instructions, but otherwise, it should work.

  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 JRF}
{ Derived from part of JRFPAS by J.R. Ferguson, Amsterdam, The Netherlands
  e-mail: j.r.ferguson@iname.com   web: www.xs4all.nl/~ferguson }
const HalfPi = 0.5*Pi ;

function ArcSin2(x : float) : float ;
begin
  if x=1.0 then ArcSin2 := HalfPi
    else if x=-1.0 then ArcSin2 := -HalfPi
    else ArcSin2 := ArcTan(x/sqrt(1.0-sqr(x))) ;
  end ;

function ArcCos2(x : float) : float ;
begin
  if x=1.0 then ArcCos2 := 0.0
    else if x=-1.0 then ArcCos2 := Pi
    else ArcCos2 := HalfPi-ArcTan(x/sqrt(1.0-sqr(x))) ;
  end ;
{$ENDIF JRF}


{$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 Noe}


{$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 JRS}


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


const
{$IFOPT N+}
Q1 = 10 ; Q2 = 24 ;
{$ELSE}
Q1 =  8 ; Q2 = 18 ;
{$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('':52, ID,
  ' $N', {$IFOPT N+} '+' {$ELSE} '-' {$ENDIF},
  ' $E', {$IFOPT E+} '+' {$ELSE} '-' {$ENDIF},
  ' : RMS =', Sqrt(S/K):6:0) ;
Write('<cr>') ; Readln ;
X := 1.0 ;
Writeln('X':Q1, 'T=ArcSin2(1.0-X)':Q2, '(2.0*T/Pi)-1.0':Q2) ;
while X>1.0E-31 do begin X := 0.01*X ; Y := 1.0-X ; T := ArcSin2(Y) ;
  Writeln(X:0, #32, T, #32, (2.0*T/Pi)-1.0) end ;
{$IFNDEF WINDOWS} Write('<cr>') ; Readln ; {$ENDIF}
END.     www.merlyn.demon.co.uk