MATH

;=========================================================================
;Mathmatical Functions
;=========================================================================
;This set of functions was obtained from the
;"M[MUMPs] Draft Standard Version 18 (Millenium)
;which Ed de Moel, as MDC Proposer and Sponsor kindly supplied to me.
;Ed J.P.M. de Moel
;c/o Jacquard Systems Research
;800 Nelson Street
;Rockville Maryland 20850-2051
;Phone and fax: +1-301-762-8999
;dem...@jacquardsystems.com
;Please note that I include Ed's name here as a matter of
identification and gratitude
;only. He is not providing support for this library.
;
;Since the proposed standard's code already contains many functions,
which are not
;yet part of the standard, I adapted this code to GT.M for LINUX and
put in code that conforms
;to the current syntax . This required changes, so that THIS CODE
DOES NOT CONFROM WITH
;THE PROPOSED STANDARD AND DOES NOT BELONG TO ANY CURRENT STANDARD.
;
;NOT EVERY ROUTINE HAS BEEN TESTED AFTER IT WAS ADAPTED. THE USER
MUST
;VERIFY THE ACCURACY OF THE CODE HIMSELF. IN CASE YOU FIND ANY ERRORS
OR HAVE SUGGESTIONS
;PLEASE MAKE THEM IN THE SAME DISCUSSION THREAD IN WHICH YOU FOUND
THIS MESSAGE.
;All and any liability whatsoever is excluded. Any use is at your own
risk.
;
;All suggestions, or the posting of further functions to be included
into this
;library is welcome. This library may be copied, distributed and
reproduced, provided that
;mentioning is made of the MDC and the sponsor as stated above and
that the sharing of
;suggestions and further functions is encouraged. All new suggestions
and functions should
;also be made available on google's discussion group comp.lang.mumps,
if possible under the
;same discussion thread as this set or routines.
;To find comp.lang.mumps, go to www.google.com and select the groups
bottom, then search form MUMPS.
;
;I would also like to thank Sanchez for making their GT.M available
for LINUX on Sourceforge.com
;It is an excellent program, simple to handle and well documented:
www.sanchez.com

;
;============ CONTENTS
===================================================================
;
;Routine: Key Word
;======== ========
;--- Mathmatical Functions ---
;ABS Absolute Value
;ARCCOSS Trigonometric arccosine in radians, optimized for
speed
;ARCOS Trigonometric arccosine in radians, optimized for
precision
;ARCCOSH Hyperbolic Arccosine in radians
;ARCCOT Trigonomic arccotangent in radians
;ARCCOTH Hyperbolic arccotangent in radians
;ARCCSC Trigonometric arccosecant in radians
;ARCSEC Trigonometric arcsecant in radians
;ARCSINS Trigonometric arcsine in radians optimized for speed
;ARCSIN Trigonometric arcsine in radians optimized for
precision
;ARCSINH Hyperbolic arcsine in radians
;ARCTAN Trigonometric arctangent in radians
;ARCTANH Hyperbolic artangent in radians
;CABS Absolute value of a complex number
;CADD Sum of two complex numbers
;CCOS Trigonometric consine in radians
;CDIV Division of complex numbers
;CEXP Raise e to the power of a complex number
;CLOG Return the logarithm of a complex number
;CMUL Multiply two complex numbers
;COMPLEX Returns the complex representation of the number
specified
;CONJUG Conjugate of a complex number
;COSS Trigonometric cosine in radians, optimized for speed
;COS Trigonometric cosine in radians, optimized for
precision
;COSH Hyperbolic cosine in radians
;COT Trigonometric cotangent in radians
;COTH Hyperbolic cotangent in radians
;CPOWER Raise one complex number to the power of another
complex number
;CSC Trigonometric cosecant of a complex number in
radians
;CSCH Hyperbolic cosecant of a complex number in radians
;CSIN Trigonometric sine of a complex number in radians
;CSUB Subtract one complex number from the other
;DECDMS Convert an angle specified in degrees into the
degree:minutes:seconds notation
;DEGRAD Convert an angle specified in degrees into raidans.
;DEMSDEC Convert an angle specified in degree:minutes:seconds
into degrees
;E Retrun the value of e
;EXP Raise the value of e to the power of a specified
number
;LOG Naperian logarithm ("natural" logarithm)
;LOG10 Briggsian logarithm (logarithm with the base 10)
;MTXADD Add one matrix to the other
;MTXCOF Cofactor
;MTXCOPY Copy a matrix
;MTXDET Determinant of a matrix
;MTXEQU Matrix equation
;MTXINV Invert a matrix
;MTXMUL Multiply one matrix with the other
;MYXSCA Multiply a scalar value with a matrix
;MTXSUB Subtract one matrix from the other
;MTXTRP Transposes one matric into the other
;MTXUNIT Create matrix as a unit matrix
;PI Returns the value of pi
;RADDEG Converts an angle specified in radians into degrees
;SEC Trigonometric secant
;SECH Hyperbolic secant
;SIGN Returns -1 if argument smaller than 0, 0 if argument
0 and +1 if argument greater than 0
;SINS Trigonometric sine, optimized for speed
;SIN Trigonometric sine, optimized for precision
;SINH Hyperbolic sine
;SQRT Root
;TAN Trigonometric tangent
;TANH Hyperbolic tangent
;
;============ Mathmatical Functions
==========================================

ABS(X) Quit $Translate(+X,"-")


ARCCOSS(X) ;
; This version of the function is
; optimized for speed, not for precision.
; The 'precision' parameter is not supported,
; and the precision is at best 2 in 10**-8.
;
New A,N,R,SIGN,XX
If X<-1 Set $Ecode=",M28,"
If X>1 Set $Ecode=",M28,"
Set SIGN=1 Set:X<0 X=-X,SIGN=-1
Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
Set A(6)=0.0066700901,A(7)=-0.0012624911
Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
Set R=$$SQRT^MATH(1-X,11)*R
Quit R*SIGN

ARCCOS(X,PREC) ;
; This version of the function is
; optimized for precision.
;
New L,LIM,K,SIG,SIGS,VALUE
;
If X<-1 Set $Ecode=",M28,"
If X>1 Set $Ecode=",M28,"
Set PREC=$Get(PREC,11)
If $Translate(X,"-")=1 Quit 0
;
Set SIG=$Select(X<0:-1,1:1),VALUE=1-(X*X)
Set X=$$SQRT^MATH(VALUE,PREC)
If $Translate(X,"-")=1 Do Quit VALUE
. Set VALUE=$$PI^MATH()/2*X
. Quit
;
If X>0.9 Do Quit VALUE
. Set SIGS=$Select(X<0:-1,1:1)
. Set VALUE=1/(1/X/X-1)
. Set X=$$SQRT^MATH(VALUE,PREC)
. Set VALUE=$$ARCTAN^MATH(X,PREC)*SIGS
. Quit
Set (VALUE,L)=X
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM)
. Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
. Quit
Quit $Select(SIG<0:$$PI^MATH()-VALUE,1:VALUE)



ARCCOSH(X,PREC) ;
If X<1 Set $Ecode=",M28,"
New SQ
Set PREC=$Get(PREC,11)
Set SQ=$$SQRT^MATH(X*X-1,PREC)
Quit $$LOG^MATH(X+SQ,PREC)



ARCCOT(X,PREC) ;
Set PREC=$Get(PREC,11)
Set X=1/X
Quit $$ARCTAN^MATH(X,PREC)



ARCCOTH(X,PREC) ;
New L1,L2
Set PREC=$Get(PREC,11)
Set L1=$$LOG^MATH(X+1,PREC)
Set L2=$$LOG^MATH(X-1,PREC)
Quit L1-L2/2



ARCCSC(X,PREC) ;
Set PREC=$Get(PREC,11)
Set X=1/X
Quit $$ARCSIN^MATH(X,PREC)



ARCSEC(X,PREC) ;
Set PREC=$Get(PREC,11)
Set X=1/X
Quit $$ARCCOS^MATH(X,PREC)



ARCSINS(X) ;
; This version of the function is
; optimized for speed, not for precision.
; The 'precision' parameter is not supported,
; and the precision is at best 2 in 10**-8.
;
New A,N,R,SIGN,XX
If X<-1 Set $Ecode=",M28,"
If X>1 Set $Ecode=",M28,"
Set SIGN=1 Set:X<0 X=-X,SIGN=-1
Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
Set A(6)=0.0066700901,A(7)=-0.0012624911
Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
Set R=$$SQRT^MATH(1-X,11)*R
Set R=$$PI^MATH()/2-R
Quit R*SIGN



ARCSIN(X,PREC) ;
;Option 2, optimized for precision, not speed.
New L,LIM,K,SIGS,VALUE
Set PREC=$Get(PREC,11)
If $Translate(X,"-")=1 Do Quit VALUE
. Set VALUE=$$PI^MATH()/2*X
. Quit
If X>0.99999 Do Quit VALUE
. Set SIGS=$Select(X<0:-1,1:1)
. Set VALUE=1/(1/X/X-1)
. Set X=$$SQRT^MATH(VALUE,PREC)
. Set VALUE=$$ARCTAN^MATH(X,PREC)*SIGS
. Quit
Set (VALUE,L)=X
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM)
. Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
. Quit
Quit VALUE



ARCSINH(X,PREC) ;
If X<1 Set $Ecode=",M28,"
New SQ
Set PREC=$Get(PREC,11)
Set SQ=$$SQRT^MATH(X*X+1,PREC)
Quit $$LOG^MATH(X+SQ,PREC)



ARCTAN(X,PREC) ;
New FOLD,HI,L,LIM,LO,K,SIGN,SIGS,SIGT,VALUE
Set PREC=$Get(PREC,11)
Set LO=0.0000000001,HI=9999999999
Set SIGT=$Select(X<0:-1,1:1),X=$Translate(X,"-")
Set X=$Select(X<LO:LO,X>HI:HI,1:X)
Set FOLD=$Select(X'<1:0,1:1)
Set X=$Select(FOLD:1/X,1:X)
Set L=X,VALUE=$$PI^MATH()/2-(1/X),SIGN=1
If X<1.3 Do Quit VALUE
. Set X=$Select(FOLD:1/X,1:X),VALUE=1/((1/X/X)+1)
. Set X=$$SQRT^MATH(VALUE,PREC)
. If $Translate(X,"-")=1 Do Quit
. . Set VALUE=$$PI^MATH()/2*X
. . Quit
. If X>0.9 Do Quit
. . Set SIGS=$Select(X<0:-1,1:1)
. . Set VALUE=1/(1/X/X-1)
. . Set X=$$SQRT^MATH(VALUE)
. . Set VALUE=$$ARCTAN(X,10)
. . Set VALUE=VALUE*SIGS
. . Quit
. Set (VALUE,L)=X
. Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
. For K=3:2 Do Quit:($Translate(L,"-")<LIM)
. . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
. . Quit
. Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
. Quit
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(1/L,"-")<LIM)
. Set L=L*X*X,VALUE=VALUE+(1/(K*L)*SIGN)
. Set SIGN=SIGN*-1
. Quit
Set VALUE=$Select(FOLD:$$PI^MATH()/2-VALUE,1:VALUE)
Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
Quit VALUE

ARCTANH(X,PREC) ;
If X<-1 Set $Ecode=",M28,"
If X>1 Set $Ecode=",M28,"
Set PREC=$Get(PREC,11)
Quit $$LOG^MATH(1+X/(1-X),PREC)/2

CABS(Z) ;
New ZRE,ZIM
Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
Quit $$SQRT^MATH(ZRE*ZRE+(ZIM*ZIM))



CADD(X,Y) ;
New XRE,XIM,YRE,YIM
Set XRE=+X,XIM=+$Piece(X,"%",2)
Set YRE=+Y,YIM=+$Piece(Y,"%",2)
Quit XRE+YRE_"%"_(XIM+YIM)



CCOS(Z,PREC) ;
New E1,E2,IA
Set PREC=$Get(PREC,11)
Set IA=$$CMUL^MATH(Z,"0%1")
Set E1=$$CEXP^MATH(IA,PREC)
Set IA=-IA_"%"_(-$Piece(IA,"%",2))
Set E2=$$CEXP^MATH(IA,PREC)
Set IA=$$CADD^MATH(E1,E2)
Quit $$CMUL^MATH(IA,"0.5%0")



CDIV(X,Y) ;
New D,IM,RE,XIM,XRE,YIM,YRE
Set XRE=+X,XIM=+$Piece(X,"%",2)
Set YRE=+Y,YIM=+$Piece(Y,"%",2)
Set D=YRE*YRE+(YIM*YIM)
Set RE=XRE*YRE+(XIM*YIM)/D
Set IM=XIM*YRE-(XRE*YIM)/D
Quit RE_"%"_IM



CEXP(Z,PREC) ;
New R,ZIM,ZRE
Set PREC=$Get(PREC,11)
Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
Set R=$$EXP^MATH(ZRE,PREC)
Quit R*$$COS^MATH(ZIM,PREC)_"%"_(R*$$SIN^MATH(ZIM,PREC))



CLOG(Z,PREC) ;
New ABS,ARG,ZIM,ZRE
Set PREC=$Get(PREC,11)
Set ABS=$$CABS^MATH(Z)
Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
Set ARG=$$ARCTAN^MATH(ZIM/ZRE,PREC)
Quit $$LOG^MATH(ABS,PREC)_"%"_ARG



CMUL(X,Y) ;
New XIM,XRE,YIM,YRE
Set XRE=+X,XIM=+$Piece(X,"%",2)
Set YRE=+Y,YIM=+$Piece(Y,"%",2)
Quit XRE*YRE-(XIM*YIM)_"%"_(XRE*YIM+(XIM*YRE))



COMPLEX(X) Quit +X_"%0"



CONJUG(Z) ;
New ZIM,ZRE
Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
Quit ZRE_"%"_(-ZIM)



COSS(X) ;
; This version of the function is
; optimized for speed, not for precision.
; The 'precision' parameter is not supported,
; and the precision is at best 1 in 10**-9.
; Note that this function does not accept its
; parameter in degrees, minutes and seconds.
;
New A,N,PI,R,SIGN,XX
;
; This approximation only works for 0 <= x <= pi/2
; so reduce angle to correct quadrant.
;
Set PI=$$PI^MATH(),X=X#(PI*2),SIGN=1
Set:X>PI X=2*PI-X
Set:X*2>PI X=PI-X,SIGN=-1
;
Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
Quit R*SIGN



COS(X,PREC) ;
New L,LIM,K,SIGN,VALUE
; Option 2, optimized for precision, not speed.
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=1,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
. Quit
Quit VALUE



COSH(X,PREC) ;
New E,F,I,P,R,T,XX
Set PREC=$Get(PREC,11)+1
Set @("E=1E-"_PREC)
Set XX=X*X,F=1,(P,R,T)=1,I=1
For Set T=T*XX,F=I+1*I*F,R=T/F+R,P=P-R/R,I=I+2 If -E<P,P<E Quit
Quit R

COT(X,PREC) ;
New C,L,LIM,K,SIGN,VALUE
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
Set (VALUE,L)=1,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
. Quit
Set C=VALUE
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=X,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
. Quit
If 'VALUE Quit "INFINITE"
Quit VALUE=C/VALUE

COTH(X,PREC) ;
New SINH
If 'X Quit "INFINITE"
;
Set PREC=$Get(PREC,11)
Set SINH=$$SINH^MATH(X,PREC)
If 'SINH Quit "INFINITE"
Quit $$COSH^MATH(X,PREC)/SINH



CPOWER(Z,N,PREC) ;
New AR,NIM,NRE,PHI,PI,R,RHO,TH,ZIM,ZRE
;
Set PREC=$Get(PREC,11)
Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
Set NRE=+N,NIM=+$Piece(N,"%",2)
If 'ZRE,'ZIM,'NRE,'NIM Set $Ecode=",M28,"
;
If 'ZRE,'ZIM Quit "0%0"
;
Set PI=$$PI^MATH()
Set R=$$SQRT^MATH(ZRE*ZRE+(ZIM*ZIM),PREC)
;
If ZRE Set TH=$$ARCTAN^MATH(ZIM/ZRE,PREC)
Else Set TH=$SELECT(ZIM>0:PI/2,1:-PI/2)
;
Set RHO=$$LOG^MATH(R,PREC)
Set AR=$$EXP^MATH(RHO*NRE-(TH*NIM),PREC)
Set PHI=RHO*NIM+(NRE*TH)
Quit AR*$$COS^MATH(PHI,PREC)_"%"_(AR*$$SIN^MATH(PHI,PREC))

CSC(X,PREC) ;
New L,LIM,K,SIGN,VALUE
;
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
;
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=X,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
. Quit
If 'VALUE Quit "INFINITE"
Quit 1/VALUE

CSCH(X,PREC) Quit 1/$$SINH^MATH(X,$Get(PREC,11))

CSIN(Z,PREC) ;
New IA,E1,E2
;
Set PREC=$Get(PREC,11)
;
Set IA=$$CMUL^MATH(Z,"0%1")
Set E1=$$CEXP^MATH(IA,PREC)
Set IA=-IA_"%"_(-$Piece(IA,"%",2))
Set E2=$$CEXP^MATH(IA,PREC)
Set IA=$$CSUB^MATH(E1,E2)
Set IA=$$CMUL^MATH(IA,"0.5%0")
Quit $$CMUL^MATH("0%-1",IA)



CSUB(X,Y) ;
New XIM,XRE,YIM,YRE
Set XRE=+X,XIM=+$Piece(X,"%",2)
Set YRE=+Y,YIM=+$Piece(Y,"%",2)
Quit XRE-YRE_"%"_(XIM-YIM)



DECDMS(X,PREC) New T
Set PREC=$Get(PREC,5)
Set X=X#360*3600
Set T=PREC-$Length(X\1)
Set X=+$Justify(X,0,$Select(T'<0:T,1:0))
Quit X\3600_":"_(X\60#60)_":"_(X#60)



DEGRAD(X) Quit X*3.14159265358979/180



DMSDEC(X) ;
Quit $Piece(X,":")+($Piece(X,":",2)/60)+($Piece(X,":",3)/3600)



E() Quit 2.71828182845905



EXP(X,PREC) ;
New L,LIM,K,VALUE
Set PREC=$Get(PREC,11)
Set L=X,VALUE=X+1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=2:1 Set L=L*X/K,VALUE=VALUE+L Quit:($Translate(L,"-")<LIM)
Quit VALUE



LOG(X,PREC) ;
New L,LIM,M,N,K,VALUE
If X'>0 Set $Ecode=",M28,"
Set PREC=$Get(PREC,11)
Set M=1
;
For N=0:1 Quit:(X/M)<10 Set M=M*10
;
If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
Set X=X/M
Set X=(X-1)/(X+1),(VALUE,L)=X
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
Set VALUE=VALUE*2+(N*2.30258509298749)
Quit VALUE



LOG10(X,PREC) ;
New L,LIM,M,N,K,VALUE
If X'>0 Set $Ecode=",M28,"
Set PREC=$Get(PREC,11)
Set M=1
For N=0:1 Quit:(X/M)<10 Set M=M*10
If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
Set X=X/M
Set X=(X-1)/(X+1),(VALUE,L)=X
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
Set VALUE=VALUE*2+(N*2.30258509298749)
Quit VALUE/2.30258509298749



MTXADD(A,B,R,ROWS,COLS) ;
; Add A[ROWS,COLS] to B[ROWS,COLS],
; result goes to R[ROWS,COLS]
If $Data(A)<10 Quit 0
If $Data(B)<10 Quit 0
If $Get(ROWS)<1 Quit 0
If $Get(COLS)<1 Quit 0
;
New ROW,COL,ANY
For ROW=1:1:ROWS For COL=1:1:COLS Do
. Kill R(ROW,COL) Set ANY=0
. Set:$Data(A(ROW,COL))#2 ANY=1
. Set:$Data(B(ROW,COL))#2 ANY=1
. Set:ANY R(ROW,COL)=$Get(A(ROW,COL))+$Get(B(ROW,COL))
. Quit
Quit 1



MTXCOF(A,I,K,N) ;
; Compute cofactor for element [i,k]
; in matrix A[N,N]
New T,R,C,RR,CC
Set CC=0 For C=1:1:N Do:C'=K
. Set CC=CC+1,RR=0
. For R=1:1:N Set:R'=I RR=RR+1,T(RR,CC)=$Get(A(R,C))
. Quit
Quit $$MTXDET^MATH(.T,N-1)



MTXCOPY(A,R,ROWS,COLS) ;
; Copy A[ROWS,COLS] to R[ROWS,COLS]
If $Data(A)<10 Quit 0
If $Get(ROWS)<1 Quit 0
If $Get(COLS)<1 Quit 0
;
New ROW,COL
For ROW=1:1:ROWS For COL=1:1:COLS Do
. Kill R(ROW,COL)
. Set:$Data(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)
. Quit
Quit 1



MTXDET(A,N) ;
; Compute determinant of matrix A[N,N]
If $Data(A)<10 Quit ""
If $Get(N)<1 Quit ""
;
; First the simple cases
;
If N=1 Quit $Get(A(1,1))
If N=2 Quit $Get(A(1,1))*$Get(A(2,2))-($Get(A(1,2))*$Get(A(2,1)))
;
New DET,I,SIGN
;
; Det A = sum (k=1:n) element (i,k) * cofactor [i,k]
;
Set DET=0,SIGN=1
For I=1:1:N Do
. Set DET=$Get(A(1,I))*$$MTXCOF^MATH(.A,1,I,N)*SIGN+DET
. Set SIGN=-SIGN
. Quit
Quit DET



MTXEQU(A,B,R,N,M) ;
; Solve matrix equation A [M,M] * R [M,N] = B [M,N]
If $Get(M)<1 Quit ""
If $Get(N)<1 Quit ""
If '$$MTXDET^MATH(.A,M) Quit 0
;
New I,I1,J,J1,J2,K,L,T,T1,T2,TEMP,X
;
Set X=$$MTXCOPY^MATH(.A,.T,N,N)
Set X=$$MTXCOPY^MATH(.B,.R,N,M)
;
; Reduction of matrix A
; Steps of reduction are counted by index K
;
For K=1:1:N-1 Do
. ;
. ; Search for largest coefficient of T
. ; (denoted by TEMP)
. ; in first column of reduced system
. ;
. Set TEMP=0,J2=K
. For J1=K:1:N Do
. . Quit:$TRanslate($Get(T(J1,K)),"-")>$TRanslate(TEMP,"-")
. . Set TEMP=T(J1,K),J2=J1
. . Quit
. ;
. ; Exchange row number K with row number J2,
. ; if necessary
. ;
. Do:J2'=K
. . ;
. . For J=K:1:N Do
. . . Set T1=$Get(T(K,J)),T2=$Get(T(J2,J))
. . . Kill T(K,J),T(J2,J)
. . . If T1'="" Set T(J2,J)=T1
. . . If T2'="" Set T(K,J)=T2
. . . Quit
. . For J=1:1:M Do
. . . Set T1=$Get(R(K,J)),T2=$Get(R(J2,J))
. . . Kill R(K,J),R(J2,J)
. . . If T1'="" Set R(J2,J)=T1
. . . If T2'="" Set R(K,J)=T2
. . . Quit
. . Quit
. ;
. ; Actual reduction
. ;
. For I=K+1:1:N Do
. . For J=K+1:1:N Do
. . . Quit:'$Get(T(K,K))
. . . Set T(I,J)=-$Get(T(K,J))*$Get(T(I,K))/T(K,K)+$Get(T(I,J))
. . . Quit
. . For J=1:1:M Do
. . . Quit:'$Get(T(K,K))
. . . Set R(I,J)=-$Get(R(K,J))*$Get(T(I,K))/T(K,K)+$Get(R(I,J))
. . . Quit
. . Quit
. Quit
;
; Backsubstitution
;
For J=1:1:M Do
. If $Get(T(N,N)) Set R(N,J)=$Get(R(N,J))/T(N,N)
. If N-1>0 For I1=1:1:N-1 Do
. . Set I=N-I1
. . For L=I+1:1:N Do
. . . Set R(I,J)=-$Get(T(I,L))*$Get(R(L,J))+$Get(R(I,J))
. . . Quit
. . If $Get(T(I,I)) Set R(I,J)=$Get(R(I,J))/$Get(T(I,I))
. . Quit
. Quit
Quit $Select(M=N:$$MTXDET^MATH(.R,M),1:1)



MTXINV(A,R,N) ;
; Invert A[N,N], result goes to R[N,N]
If $Data(A)<10 Quit 0
If $Get(N)<1 Quit 0
;
New T,X
Set X=$$MTXUNIT^MATH(.T,N)
Quit $$MTXEQU^MATH(.A,.T,.R,N,N)



MTXMUL(A,B,R,M,L,N) ;
; Multiply A[M,L] by B[L,N], result goes to R[M,N]
If $Data(A)<10 Quit 0
If $Data(B)<10 Quit 0
If $Get(L)<1 Quit 0
If $Get(M)<1 Quit 0
If $Get(N)<1 Quit 0
;
New I,J,K,SUM,ANY
For I=1:1:M For J=1:1:N Do
. Set (SUM,ANY)=0
. Kill R(I,J)
. For K=1:1:L Do
. . Set:$Data(A(I,K))#2 ANY=1
. . Set:$Data(B(K,J))#2 ANY=1
. . Set SUM=$Get(A(I,K))*$Get(B(K,J))+SUM
. . Quit
. Set:ANY R(I,J)=SUM
. Quit
Quit 1



MTXSCA(A,R,ROWS,COLS,S) ;
; Multiply A[ROWS,COLS] with the scalar S,
; result goes to R[ROWS,COLS]
If $Data(A)<10 Quit 0
If $Get(ROWS)<1 Quit 0
If $Get(COLS)<1 Quit 0
If '($Data(S)#2) Quit 0
;
New ROW,COL
For ROW=1:1:ROWS For COL=1:1:COLS Do
. Kill R(ROW,COL)
. Set:$Data(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)*S
. Quit
Quit 1



MTXSUB(A,B,R,ROWS,COLS) ;
; Subtract B[ROWS,COLS] from A[ROWS,COLS],
; result goes to R[ROWS,COLS]
If $Data(A)<10 Quit 0
If $Data(B)<10 Quit 0
If $Get(ROWS)<1 Quit 0
If $Get(COLS)<1 Quit 0
;
New ROW,COL,ANY
For ROW=1:1:ROWS For COL=1:1:COLS Do
. Kill R(ROW,COL) Set ANY=0
. Set:$Data(A(ROW,COL))#2 ANY=1
. Set:$Data(B(ROW,COL))#2 ANY=1
. Set:ANY R(ROW,COL)=$Get(A(ROW,COL))-$Get(B(ROW,COL))
. Quit
Quit 1



MTXTRP(A,R,M,N) ;
; Transpose A[M,N], result goes to R[N,M]
If $Data(A)<10 Quit 0
If $Get(M)<1 Quit 0
If $Get(N)<1 Quit 0
;
New I,J,K,D1,V1,D2,V2
For I=1:1:M+N-1 For J=1:1:I+1\2 Do
. Set K=I-J+1
. If K=J Do Quit
. . Set V1=$Get(A(J,J)),D1=$Data(A(J,J))#2
. . If J'>N,J'>M Kill R(J,J) Set:D1 R(J,J)=V1
. . Quit
. ;
. Set V1=$Get(A(K,J)),D1=$Data(A(K,J))#2
. Set V2=$Get(A(J,K)),D2=$Data(A(J,K))#2
. If K'>M,J'>N Kill R(K,J) Set:D2 R(K,J)=V2
. If J'>M,K'>N Kill R(J,K) Set:D1 R(J,K)=V1
. Quit
Quit 1



MTXUNIT(R,N,SPARSE) ;
; Create a unit matrix R[N,N]
If $Get(N)<1 Quit 0
;
New ROW,COL
For ROW=1:1:N For COL=1:1:N Do
. Kill R(ROW,COL)
. If $Get(SPARSE) Quit:ROW'=COL
. Set R(ROW,COL)=$Select(ROW=COL:1,1:0)
. Quit
Quit 1



PI() Quit 3.14159265358979



RADDEG(X) Quit X*180/3.14159265358979

SEC(X,PREC) ;
New L,LIM,K,SIGN,VALUE
;
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=1,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
. Quit
If 'VALUE Quit "INFINITE"
Quit 1/VALUE



SECH(X,PREC) Quit 1/$$COSH^MATH(X,$Get(PREC,11))



SIGN(X) Quit $Select(X<0:-1,X>0:1,1:0)



SINS(X) ;
;
; This version of the function is
; optimized for speed, not for precision.
; The 'precision' parameter is not supported,
; and the precision is at best 1 in 10**-9.
; Note that this function does not accept its
; parameter in degrees, minutes and seconds.
;
New A,N,PI,R,SIGN,XX
;
; This approximation only works for 0 <= x <= pi/2
; so reduce angle to correct quadrant.
;
Set PI=$$PI^MATH(),X=X#(PI*2),SIGN=1
Set:X>PI X=2*PI-X,SIGN=-1
Set:X*2<PI X=PI-X
Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
Quit R*SIGN



SIN(X,PREC) ;
New L,LIM,K,SIGN,VALUE
;
; Option 2, optimized for precision, not speed
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=X,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
. Quit
Quit VALUE

SINH(X,PREC) ;
New E,F,I,P,R,T,XX
;
Set PREC=$Get(PREC,11)+1
Set @("E=1E-"_PREC)
Set XX=X*X,F=1,I=2,(P,R,T)=X
For Set T=T*XX,F=I+1*I*F,R=T/F+R,P=P-R/R,I=I+2 If -E<P,P<E Quit
Quit R

SQRT(X,PREC) ;
If X<0 Set $Ecode=",M28,"
If X=0 Quit 0
;
Set PREC=$Get(PREC,11)
If X<1 Quit 1/$$SQRT^MATH(1/X,PREC)
;
New P,R,E
Set PREC=$Get(PREC,11)+1
Set @("E=1E-"_PREC)
;
Set R=X
For Set P=R,R=X/R+R/2,P=P-R/R If -E<P,P<E Quit
Quit R

TAN(X,PREC) ;
New L,LIM,K,S,SIGN,VALUE
;
; The official description does not mention that
; the function may also be called with the first
; parameter in degrees, minutes and seconds.
Set:X[":" X=$$DMSDEC^MATH(X)
;
Set PREC=$Get(PREC,11)
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=X,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
. Quit
Set S=VALUE
Set X=X#(2*$$PI^MATH())
Set (VALUE,L)=1,SIGN=-1
Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
. Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
. Quit
If 'VALUE Quit "INFINITE"
Quit S/VALUE



TANH(X,PREC) ;
Set PREC=$Get(PREC,11)
;
Quit $$SINH^MATH(X,PREC)/$$COSH^MATH(X,PREC)