%record %format cxf (%real     real,imag)
%externalroutinespec dhex(%longrealname x)
%record %format dcxf(%longreal real,imag)

%externalrealfnspec sin   %alias "f_sin"   ( %realname x )
%externalrealfnspec cos   %alias "f_cos"   ( %realname x )
%externalrealfnspec sinh  %alias "f_sinh"  ( %realname x )
%externalrealfnspec cosh  %alias "f_cosh"  ( %realname x )
%externalrealfnspec exp   %alias "f_exp"   ( %realname x)
%externalrealfnspec sqrt  %alias "f_sqrt"  ( %realname x )
%externalrealfnspec log   %alias "f_log"   ( %realname x )
%externalrealfnspec atan2 %alias "f_atan2" ( %realname x, y)

%externallongrealfnspec dsin   %alias "f_dsin"   ( %longrealname x )
%externallongrealfnspec dcos   %alias "f_dcos"   ( %longrealname x )
%externallongrealfnspec dsinh  %alias "f_dsinh"  ( %longrealname x )
%externallongrealfnspec dcosh  %alias "f_dcosh"  ( %longrealname x )
%externallongrealfnspec dexp   %alias "f_dexp"   ( %longrealname x )
%externallongrealfnspec dlog   %alias "f_dlog"   ( %longrealname x )
%externallongrealfnspec datan2 %alias "f_datan2" ( %longrealname x, y)
%externallongrealfnspec dsqrt  %alias "f_dsqrt"  ( %longrealname x )

%constantreal zero = 0.0,
              half = 0.5,
              SQRTHALF = R'40B504F3',
              one  = 1.0,
              two  = 2.0

%constantlongreal lzero = 0.0,
                  lhalf = 0.5,
                  DSQRTHALF = R'40B504F333F9DE65',
                  lone  = 1.0,
                  ltwo  = 2.0

!SQRTHALF= 7.0710678118654753@ -1

{------------------------------------------------------------------------}
%externalroutine MATHS ERROR ROUTINE %alias "f_mle" (%integer error number)
    !
    !
%externalroutinespec fstop %alias "f_stop" (%integer i,j)
%externalroutinespec Ndiag %alias "s#ndiag"(%integer i,j,k,l)
%integer R10

%conststring(67) %array maths errors(1:47)=                %c
{dlogsmall    = 1}  "DLOG arg negative or zero",
{dsqrtneg     = 2}  "DSQRT arg negative",
{dexplarge    = 3}  "DEXP arg greater than 174.6731",
{dexpsmall    = 4}  "DEXP arg less than -180.2182, not used",
{dsinlarge    = 5}  "DSIN arg magnitude greater than 7.205D16",
{dasinlarge   = 6}  "DASIN arg magnitude greater than 1.0",
{dcoslarge    = 7}  "DCOS arg magnitude greater than 7.205D16   ",
{dacoslarge   = 8}  "DACOS arg magnitude greater than 1.0",
{dtanlarge    = 9}  "DTAN/DCOTAN arg magnitude greater than 3.521D15",
{dcotansmall  = 10} "DCOTAN arg small in magnitude",
{darcsinlarge = 11} "Not used",
{darccoslarge = 12} "Not used",
{darctan2zero = 13} "DATAN2 args zero",
{dsinhlarge   = 14} "DSINH arg magnitude greater than 175.3662 ",
{dcoshlarge   = 15} "DCOSH arg magnitude greater than 175.3662 ",
{dpowerneg    = 16} "Negative D.P. value raised to a non-integer power",
{dpowerzero   = 17} "D.P. zero raised to non-positive power   ",
{ sqrtneg     = 18} "SQRT arg negative",
{ explarge    = 19} "EXP arg greater than 174.6731",
{ expsmall    = 20} "EXP arg less than -180.2182, not used",
{ alogsmall   = 21} "ALOG arg negative or zero",
{ powerzero   = 22} "REAL zero raised to non-positive power   ",
{ powerneg    = 23} "Negative REAL value raised to a non-integer power",
{ powerbig    = 24} "REAL value raised to too large a REAL power ",
{ powersmall  = 25} "REAL value to REAL power underflows. not used",
{ sinlarge    = 26} "SIN arg magnitude greater than 5.274E7",
{ coslarge    = 27} "COS arg magnitude greater than 5.274E7",
{ tanlarge    = 28} "TAN/COTAN arg magnitude greater than 2.633E7",
{ asinlarge   = 29} "ASIN arg magnitude greater than 1.0",
{ acoslarge   = 30} "ACOS arg magnitude greater than 1.0",
{ arctan2zero = 31} "ATAN2 args zero",
{ sinhlarge   = 32} "SINH arg magnitude greater than 175.3662 ",
{ coshlarge   = 33} "COSH arg magnitude greater than 175.3662 ",
{dgamlarge    = 34} "DGAMMA arg magnitude greater than 57.0",
{dgamnegint   = 35} "DGAMMA arg near zero or negative integer",
{dlgamneg     = 36} "DLGAMA arg is negative",
{dlgamlarge   = 37} "DLGAMA arg is greater than 4.29D73",
{ipowexpneg   = 38} "Integer raised to negative integer power",
{ipowerzero   = 39} "Integer zero raised to non-positive power",
{ipowlarge    = 40} "Integer to integer power overflows",
{ cotansmall  = 41} "COTAN arg small in magnitude",
{ gamlarge    = 42} "GAMMA arg magnitude greater than 57.0",
{ gamnegint   = 43} "GAMMA arg near zero or negative integer",
{ lgamneg     = 44} "ALGAMA arg is negative",
{ lgamlarge   = 45} "ALGAMA arg is greater than 4.29E73",
{ dpowerlarge = 46} "D.P. value raised to too large a D.P. power",
{ cxpowzero   = 47} "COMPLEX zero raised to non-positive INTEGER power"


!Print the Error Message:
   !
   !
 print string ("
Maths Library Error"); write(error number,3); print string (": ");

%if error number>=  1 %and                                                %c
    error number<= 47 %then print string (maths errors (error number))    %c
                      %else print string ("no text")

newline

!Now Stop With Diagnostics:
          !
!          fstop (-255, {skip} 2 {stack frames})
      *ST_10,R10
      Ndiag(0,integer(R10+40),0,0)
      %stop

%end; !of MATHS ERROR ROUTINE



%externalroutine cmpcos %alias "f_ccos" (%record(cxf) %name r, z)
{------------------------------------------------------------------------}

  r_real = cos(z_real) * cosh(z_imag)
  r_imag = - sin(z_real) * sinh(z_imag)
%end


{------------------------------------------------------------------------}
%externalroutine cmpsin %alias "f_csin" ( %record(cxf) %name r, z)
{------------------------------------------------------------------------}

  r_real = sin(z_real) * cosh(z_imag)
  r_imag = cos(z_real) * sinh(z_imag)

%end



{------------------------------------------------------------------------}
%externalroutine cmpexp %alias "f_cexp" ( %record(cxf) %name r, z)
{------------------------------------------------------------------------}

%real expx

  expx = exp(z_real)
  r_real = expx * cos(z_imag)
  r_imag = expx * sin(z_imag)

%end

{------------------------------------------------------------------------}
%realfn cabs( %real realpt, imagpt )
{------------------------------------------------------------------------}
  { Note: this is params by value - not to be seen by user }
  {  See: f_cabs }

  %real temp, t

  %if realpt < zero %then realpt = -realpt
  %if imagpt < zero %then imagpt = -imagpt
  %if imagpt > realpt %start
    temp = realpt
    realpt = imagpt
    imagpt = temp
  %finish
  %if (realpt+imagpt) = realpt %then %result = realpt
  temp = imagpt/realpt
  t = one + temp*temp  
  %result = realpt*sqrt( t )

%end


{------------------------------------------------------------------------}
%externalrealfunction cmpabs %alias "f_cabs" ( %record(cxf) %name z)
{------------------------------------------------------------------------}

  %result = cabs(z_real,z_imag)

%end




{------------------------------------------------------------------------}
%externalroutine cmplog %alias "f_clog" ( %record(cxf) %name r, z)
{------------------------------------------------------------------------}

%real temp

  r_imag = atan2(z_imag, z_real)
  temp   = cabs(z_real, z_imag)
  r_real = log( temp )

%end


{------------------------------------------------------------------------}
%externalroutine cmpsqrt %alias "f_csqrt" ( %record(cxf) %name r, z)
{------------------------------------------------------------------------}

%real XR,XI,H

XR= MOD(Z_REAL)
XI= Z_IMAG
H= CABS(XR,XI)
%IF XR>ONE %THEN H= HALF*XR %AND H= SQRT(H+CABS(H,XI*HALF))  %C
           %ELSE H= SQRT(XR + CABS(XR,XI))*SQRTHALF
%IF XI#ZERO %THEN XI= XI/(H+H)
%IF Z_REAL>ZERO %START
  R_REAL= H
  R_IMAG= XI
%FINISHELSE %IF XI>=ZERO %START
  R_REAL= XI
  R_IMAG= H
%FINISHELSESTART
  R_REAL= -XI
  R_IMAG= -H
%FINISH
%END





{----------------------------------------------------------------------}
%externalroutine cdcos %alias "f_cdcos" (%record(dcxf) %name r, z)
{----------------------------------------------------------------------}

  r_real =   dcos(z_real) * dcosh(z_imag)
  r_imag = - dsin(z_real) * dsinh(z_imag)
%end

{----------------------------------------------------------------------}
%externalroutine cdsin %alias "f_cdsin" ( %record(dcxf) %name r, z)
{----------------------------------------------------------------------}

  r_real = dsin(z_real) * dcosh(z_imag)
  r_imag = dcos(z_real) * dsinh(z_imag)

%end

{----------------------------------------------------------------------}
%externalroutine cdexp %alias "f_cdexp" ( %record(dcxf) %name r, z)
{----------------------------------------------------------------------}
  %longreal dexpx

  dexpx   = dexp(z_real)
  r_real = dexpx * dcos(z_imag)
  r_imag = dexpx * dsin(z_imag)

%end


{------------------------------------------------------------------}
%longrealfn zabs %alias "f_zabs"( %longreal realpt, imagpt )
{------------------------------------------------------------------}
   { Note: this is params by value - not to be seen by users }
   {  See: f_cdabs }

%longreal temp, t

  %if realpt < lzero %then realpt = -realpt
  %if imagpt < lzero %then imagpt = -imagpt
  %if imagpt > realpt %start
    temp   = realpt
    realpt = imagpt
    imagpt = temp
  %finish
  %if (realpt+imagpt) = realpt %then %result = realpt

  temp    = imagpt/realpt
  t       = lone + temp*temp { overflow! }
  %result = realpt*Dsqrt( t ) 

%end


{------------------------------------------------------------------------}
%externallongrealfunction cdabs %alias "f_cdabs" ( %record(dcxf) %name z)
{------------------------------------------------------------------------}
%LONGREAL T

T= ZABS(Z_REAL,Z_IMAG)
DHEX(T)
   %result = zabs(z_real,z_imag)

%end


{----------------------------------------------------------------------}
%externalroutine cdlog %alias "f_cdlog" ( %record(dcxf) %name r, z)
{----------------------------------------------------------------------}
  %longreal t

  r_imag = datan2(z_imag, z_real)
  t      = zabs(z_real, z_imag)
  r_real = dlog( t )

%end


{------------------------------------------------------------------}
%externalroutine cdsqrt %alias "f_cdsqrt" ( %record(dcxf) %name r, z)
{------------------------------------------------------------------}
%LONGREAL XR,XI,H

XR= MOD(Z_REAL)
XI= Z_IMAG
%IF XR>LONE %THEN H= LHALF*XR %AND H= DSQRT(H+ZABS(H,XI*LHALF))  %C
           %ELSE H= DSQRT(XR + ZABS(XR,XI))*DSQRTHALF
%IF XI#LZERO %THEN XI= XI/(H+H)
%IF Z_REAL>LZERO %START
  R_REAL= H
  R_IMAG= XI
%FINISHELSE %IF XI>=LZERO %START
  R_REAL= XI
  R_IMAG= H
%FINISHELSESTART
  R_REAL= -XI
  R_IMAG= -H
%FINISH
%END


%externalroutinespec error %alias "f_mle" (%integer type)
!*

!---------------------------------------------------------------------------}
%externalroutine powcc %alias "f_powcc" (%record(cxf)%name r,a,b)
!---------------------------------------------------------------------------}

  %real logr,logi, x, y
   x=cmpabs(a)
  logr = log(x)
  logi = atan2(a_real,a_imag)

  x    = (logr*b_real) - (logi*b_imag)
  x    = exp(x)
  y    = (logr*b_imag) + (logi*b_real)

  r_real   = x*cos(y)
  r_imag   = x*sin(y)

%end


{--------------------------------------------------------------------}
%routine zdiv (%longrealname cr,ci,%longreal ar,ai,br,bi)
{--------------------------------------------------------------------}

  %longreal ratio, den, abr, abi

  abr = br
  %if abr  < lzero %then abr = - abr
  abi = bi
  %if abi < lzero %then abi = - abi
  %if abr <= abi %start
    %if abi = lzero %then error(47) { fatal("complex division by zero"); } 
    ratio = br/bi
    den   = bi * (1 + ratio*ratio)
    cr    = (ar*ratio + ai) / den
    ci    = (ai*ratio - ar) / den
  %finishelsestart
    ratio = bi / br
    den   = br * (1 + ratio*ratio)
    cr    = (ar + ai*ratio) / den
    ci    = (ai - ar*ratio) / den
  %finish

%end

{--------------------------------------------------------------------}
%externalroutine powzi %alias "f_powzi" %c
                       (%record(dcxf) %name p, a, %integer n) 
{--------------------------------------------------------------------}
  %record(dcxf) x
  %longreal t;

  p_real = 1; p_imag = 0
 
  %return %if n=0

  %if n>0 %then x = a %else n = -n %andc
                            zdiv(x_real,x_imag,p_real,p_imag,a_real,a_imag) 

  %cycle
    %if n&1#0 %start
      t      = p_real*x_real - p_imag*x_imag
      p_imag = p_real*x_imag + p_imag*x_real
      p_real = t
    %finish
    n = n>>1
    %if n#0 %start
      t = x_real*x_real - x_imag*x_imag
      x_imag = 2*x_real*x_imag
      x_real=t
    %finishelse %exit
  %repeat
%end


{----------------------------------------------------------------------------}
%externalroutine powci %alias "f_powci"  %c
                       (%record(cxf) %name p,a ,%integer b)
{----------------------------------------------------------------------------}
   	                            { p = a**b  }
 %record(dcxf) dp,da
  
  da_real = a_real  { coerce to 64 bit real }
  da_imag = a_imag

  powzi(dp, da, b)
  p_real = dp_real  { coerce to 32 bit real}
  p_imag = dp_imag

%end






!---------------------------------------------------------------------------}
%externalroutine powzz %alias "f_powzz" (%record(dcxf)%name r,a,b)
!---------------------------------------------------------------------------}
   %longreal logr,logi, x, y
   logr = cdabs(a)
   logr = dlog(logr)
   logi = datan2(a_real,a_imag)
 
   x    = (logr*b_real) - (logi*b_imag)
   x    = dexp(x)
   y    = (logr*b_imag) + (logi*b_real)
 
   r_real   = x*dcos(y)
   r_imag   = x*dsin(y)
 
%end
!
%endoffile