1 | ;"16-Feb-1999, 16:54:35
|
---|
2 | ;"Routine Save for all M[UMPS] Library Functions
|
---|
3 | ;
|
---|
4 | ;" Unless otherwise noted, the code below
|
---|
5 | ;" was approved in document X11/95-11
|
---|
6 | ;
|
---|
7 | ;" If corrections have been applied,
|
---|
8 | ;" first the original line appears,
|
---|
9 | ;" with three semicolons at the beginning of the line.
|
---|
10 | ;
|
---|
11 | ;" Then the source of the correction is acknowledged,
|
---|
12 | ;" then the corrected line appears, followed by a
|
---|
13 | ;" line containing three semicolons.
|
---|
14 | ;
|
---|
15 | ;"Downloaded from http://www.jacquardsystems.com/Examples/lib/mlibfunc.rs
|
---|
16 | ;"on 5/21/07
|
---|
17 | ABS(X) Quit $Translate(+X,"-")
|
---|
18 | ;===
|
---|
19 | ;
|
---|
20 | ;
|
---|
21 | ARCCOS(X) ;
|
---|
22 | ;;;" ;" Number ~~
|
---|
23 | ;" Winfried Gerum (8 June 1995)
|
---|
24 | ;" Comment: This version of the function is
|
---|
25 | ;" optimized for speed, not for precision.
|
---|
26 | ;" The 'precision' parameter is not supported,
|
---|
27 | ;" and the precision is at best 2 in 10**-8.
|
---|
28 | ;;;
|
---|
29 | ;
|
---|
30 | New A,N,R,SIGN,XX
|
---|
31 | If X<-1 Set $Ecode=",M28,"
|
---|
32 | If X>1 Set $Ecode=",M28,"
|
---|
33 | Set SIGN=1 Set:X<0 X=-X,SIGN=-1
|
---|
34 | Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
|
---|
35 | Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
|
---|
36 | Set A(6)=0.0066700901,A(7)=-0.0012624911
|
---|
37 | Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
|
---|
38 | ;
|
---|
39 | ;;;" Set R=$%SQRT^MATH(1-X)*R ;" Number ~~
|
---|
40 | ;" Winfried Gerum (8 June 1995)
|
---|
41 | Set R=$%SQRT^MATH(1-X,11)*R
|
---|
42 | ;;;
|
---|
43 | ;
|
---|
44 | Quit R*SIGN
|
---|
45 | ;===
|
---|
46 | ;
|
---|
47 | ;
|
---|
48 | ARCCOS(X,PREC) ;
|
---|
49 | ;
|
---|
50 | ;;;" New L,LIM,K,SIG,SIGS ;" Number ~~
|
---|
51 | ;" Winfried Gerum (8 June 1995)
|
---|
52 | New L,LIM,K,SIG,SIGS,VALUE
|
---|
53 | ;;;
|
---|
54 | ;
|
---|
55 | If X<-1 Set $Ecode=",M28,"
|
---|
56 | If X>1 Set $Ecode=",M28,"
|
---|
57 | Set PREC=$Get(PREC,11)
|
---|
58 | ;
|
---|
59 | ;;;" If $Translate(X,"-")=1 Set VALUE=0 Quit ;" Number ~~
|
---|
60 | ;" Winfried Gerum (8 June 1995)
|
---|
61 | ;" Eli Reidler (28 June 1996)
|
---|
62 | If $Translate(X,"-")=1 Quit 0
|
---|
63 | ;;;
|
---|
64 | ;
|
---|
65 | Set SIG=$Select(X<0:-1,1:1),VALUE=1-(X*X)
|
---|
66 | ;
|
---|
67 | ;;;" Set X=$%SQRT^MATH(VALUE) ;" Number ~~
|
---|
68 | ;" Winfried Gerum (8 June 1995)
|
---|
69 | Set X=$%SQRT^MATH(VALUE,PREC)
|
---|
70 | ;;;
|
---|
71 | ;
|
---|
72 | ;;;" If $Translate(X,"-")=1 Do Quit ;" Number ~~
|
---|
73 | ;" Winfried Gerum (8 June 1995)
|
---|
74 | ;" Eli Reidler (28 June 1996)
|
---|
75 | If $Translate(X,"-")=1 Do Quit VALUE
|
---|
76 | . ;;;
|
---|
77 | . ;
|
---|
78 | . Set VALUE=$%PI^MATH()/2*X
|
---|
79 | . Quit
|
---|
80 | ;
|
---|
81 | ;;;" If X>0.9 Do Quit ;" Number ~~
|
---|
82 | ;" Winfried Gerum (8 June 1995)
|
---|
83 | ;" Eli Reidler (28 June 1996)
|
---|
84 | If X>0.9 Do Quit VALUE
|
---|
85 | . ;;;
|
---|
86 | . ;
|
---|
87 | . Set SIGS=$Select(X<0:-1,1:1)
|
---|
88 | . Set VALUE=1/(1/X/X-1)
|
---|
89 | . ;
|
---|
90 | . ;;;" Set X=$%SQRT^MATH(VALUE) ;" Number ~~
|
---|
91 | . ;" Winfried Gerum (8 June 1995)
|
---|
92 | . Set X=$%SQRT^MATH(VALUE,PREC)
|
---|
93 | . ;;;
|
---|
94 | . ;
|
---|
95 | . ;
|
---|
96 | . ;;;" Set VALUE=$%ARCTAN^MATH(X,10)*SIGS ;" Number ~~
|
---|
97 | . ;" Winfried Gerum (8 June 1995)
|
---|
98 | . Set VALUE=$%ARCTAN^MATH(X,PREC)*SIGS
|
---|
99 | . ;;;
|
---|
100 | ;
|
---|
101 | . Quit
|
---|
102 | Set (VALUE,L)=X
|
---|
103 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
104 | For K=3:2 Do Quit:($Translate(L,"-")<LIM)
|
---|
105 | . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
|
---|
106 | . Quit
|
---|
107 | Quit $Select(SIG<0:$%PI^MATH()-VALUE,1:VALUE)
|
---|
108 | ;===
|
---|
109 | ;
|
---|
110 | ;
|
---|
111 | ARCCOSH(X,PREC) ;
|
---|
112 | If X<1 Set $Ecode=",M28,"
|
---|
113 | New SQ
|
---|
114 | ;
|
---|
115 | ;;;" ;" Number ~~
|
---|
116 | ;" Winfried Gerum (8 June 1995)
|
---|
117 | ;" Alan Frank (October 1995)
|
---|
118 | Set PREC=$Get(PREC,11)
|
---|
119 | ;;;
|
---|
120 | ;
|
---|
121 | Set SQ=$%SQRT^MATH(X*X-1,PREC)
|
---|
122 | Quit $%LOG^MATH(X+SQ,PREC)
|
---|
123 | ;===
|
---|
124 | ;
|
---|
125 | ;
|
---|
126 | ARCCOT(X,PREC) ;
|
---|
127 | Set PREC=$Get(PREC,11)
|
---|
128 | Set X=1/X
|
---|
129 | Quit $%ARCTAN^MATH(X,PREC)
|
---|
130 | ;===
|
---|
131 | ;
|
---|
132 | ;
|
---|
133 | ARCCOTH(X,PREC) ;
|
---|
134 | New L1,L2
|
---|
135 | ;
|
---|
136 | ;;;" ;" Number ~~
|
---|
137 | ;" Winfried Gerum (8 June 1995)
|
---|
138 | ;" Alan Frank (October 1995)
|
---|
139 | Set PREC=$Get(PREC,11)
|
---|
140 | ;;;
|
---|
141 | ;
|
---|
142 | Set L1=$%LOG^MATH(X+1,PREC)
|
---|
143 | Set L2=$%LOG^MATH(X-1,PREC)
|
---|
144 | Quit L1-L2/2
|
---|
145 | ;===
|
---|
146 | ;
|
---|
147 | ;
|
---|
148 | ARCCSC(X,PREC) ;
|
---|
149 | Set PREC=$Get(PREC,11)
|
---|
150 | Set X=1/X
|
---|
151 | Quit $%ARCSIN^MATH(X,PREC)
|
---|
152 | ;===
|
---|
153 | ;
|
---|
154 | ;
|
---|
155 | ARCSEC(X,PREC) ;
|
---|
156 | Set PREC=$Get(PREC,11)
|
---|
157 | Set X=1/X
|
---|
158 | Quit $%ARCCOS^MATH(X,PREC)
|
---|
159 | ;===
|
---|
160 | ;
|
---|
161 | ;
|
---|
162 | ARCSIN(X) ;
|
---|
163 | ;;;" ;" Number ~~
|
---|
164 | ;" Winfried Gerum (8 June 1995)
|
---|
165 | ;" Comment: This version of the function is
|
---|
166 | ;" optimized for speed, not for precision.
|
---|
167 | ;" The 'precision' parameter is not supported,
|
---|
168 | ;" and the precision is at best 2 in 10**-8.
|
---|
169 | ;;;
|
---|
170 | ;
|
---|
171 | New A,N,R,SIGN,XX
|
---|
172 | If X<-1 Set $Ecode=",M28,"
|
---|
173 | If X>1 Set $Ecode=",M28,"
|
---|
174 | Set SIGN=1 Set:X<0 X=-X,SIGN=-1
|
---|
175 | Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
|
---|
176 | Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
|
---|
177 | Set A(6)=0.0066700901,A(7)=-0.0012624911
|
---|
178 | Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
|
---|
179 | ;
|
---|
180 | ;;;" Set R=$%SQRT^MATH(1-X)*R ;" Number ~~
|
---|
181 | ;" Winfried Gerum (8 June 1995)
|
---|
182 | Set R=$%SQRT^MATH(1-X,11)*R
|
---|
183 | ;;;
|
---|
184 | ;
|
---|
185 | Set R=$%PI^MATH()/2-R
|
---|
186 | Quit R*SIGN
|
---|
187 | ;===
|
---|
188 | ;
|
---|
189 | ;
|
---|
190 | ARCSIN(X,PREC) ;
|
---|
191 | New L,LIM,K,SIGS,VALUE
|
---|
192 | Set PREC=$Get(PREC,11)
|
---|
193 | ;
|
---|
194 | ;;;" If $Translate(X,"-")=1 Do Quit ;" Number ~~
|
---|
195 | ;" Winfried Gerum (8 June 1995)
|
---|
196 | ;" Eli Reidler (28 June 1996)
|
---|
197 | If $Translate(X,"-")=1 Do Quit VALUE
|
---|
198 | . ;;;
|
---|
199 | . ;
|
---|
200 | . Set VALUE=$%PI^MATH()/2*X
|
---|
201 | . Quit
|
---|
202 | ;
|
---|
203 | ;;;" If X>0.99999 Do Quit ;" Number ~~
|
---|
204 | ;" Winfried Gerum (8 June 1995)
|
---|
205 | ;" Eli Reidler (28 June 1996)
|
---|
206 | If X>0.99999 Do Quit VALUE
|
---|
207 | . ;;;
|
---|
208 | . ;
|
---|
209 | . Set SIGS=$Select(X<0:-1,1:1)
|
---|
210 | . Set VALUE=1/(1/X/X-1)
|
---|
211 | . ;
|
---|
212 | . ;;;" Set X=$%SQRT^MATH(VALUE) ;" Number ~~
|
---|
213 | . ;" Winfried Gerum (8 June 1995)
|
---|
214 | . Set X=$%SQRT^MATH(VALUE,PREC)
|
---|
215 | . ;;;
|
---|
216 | . ;
|
---|
217 | . ;;;" Set VALUE=$%ARCTAN^MATH(X,10)*SIGS ;" Number ~~
|
---|
218 | . ;" Winfried Gerum (8 June 1995)
|
---|
219 | . Set VALUE=$%ARCTAN^MATH(X,PREC)*SIGS
|
---|
220 | . ;;;
|
---|
221 | . ;
|
---|
222 | . Quit
|
---|
223 | Set (VALUE,L)=X
|
---|
224 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
225 | For K=3:2 Do Quit:($Translate(L,"-")<LIM)
|
---|
226 | . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
|
---|
227 | . Quit
|
---|
228 | Quit VALUE
|
---|
229 | ;===
|
---|
230 | ;
|
---|
231 | ;
|
---|
232 | ARCSINH(X,PREC) ;
|
---|
233 | If X<1 Set $Ecode=",M28,"
|
---|
234 | New SQ
|
---|
235 | ;
|
---|
236 | ;;;" ;" Number ~~
|
---|
237 | ;" Winfried Gerum (8 June 1995)
|
---|
238 | ;" Alan Frank (October 1995)
|
---|
239 | Set PREC=$Get(PREC,11)
|
---|
240 | ;;;
|
---|
241 | ;
|
---|
242 | Set SQ=$%SQRT^MATH(X*X+1,PREC)
|
---|
243 | Quit $%LOG^MATH(X+SQ,PREC)
|
---|
244 | ;===
|
---|
245 | ;
|
---|
246 | ;
|
---|
247 | ARCTAN(X,PREC) ;
|
---|
248 | New FOLD,HI,L,LIM,LO,K,SIGN,SIGS,SIGT,VALUE
|
---|
249 | Set PREC=$Get(PREC,11)
|
---|
250 | Set LO=0.0000000001,HI=9999999999
|
---|
251 | Set SIGT=$Select(X<0:-1,1:1),X=$Translate(X,"-")
|
---|
252 | Set X=$Select(X<LO:LO,X>HI:HI,1:X)
|
---|
253 | ;
|
---|
254 | ;;;" Set FOLD=$Select(X'<1:0,1:1), ;" Number ~~
|
---|
255 | ;" Eli Reidler (28 June 1996)
|
---|
256 | Set FOLD=$Select(X'<1:0,1:1)
|
---|
257 | ;;;
|
---|
258 | ;
|
---|
259 | Set X=$Select(FOLD:1/X,1:X)
|
---|
260 | Set L=X,VALUE=$%PI^MATH()/2-(1/X),SIGN=1
|
---|
261 | ;
|
---|
262 | ;;;" If X<1.3 Do Quit ;" Number ~~
|
---|
263 | ;" Winfried Gerum (8 June 1995)
|
---|
264 | ;" Eli Reidler (28 June 1996)
|
---|
265 | If X<1.3 Do Quit VALUE
|
---|
266 | . ;;;
|
---|
267 | . ;
|
---|
268 | . Set X=$Select(FOLD:1/X,1:X),VALUE=1/((1/X/X)+1)
|
---|
269 | . ;
|
---|
270 | . ;;;" Set $%SQRT^MATH(VALUE) ;" Number ~~
|
---|
271 | . ;" Winfried Gerum (8 June 1995)
|
---|
272 | . ;" Eli Reidler (28 June 1996)
|
---|
273 | . Set X=$%SQRT^MATH(VALUE,PREC)
|
---|
274 | . ;;;
|
---|
275 | . ;
|
---|
276 | . If $Translate(X,"-")=1 Do Quit
|
---|
277 | . . Set VALUE=$%PI^MATH()/2*X
|
---|
278 | . . Quit
|
---|
279 | . If X>0.9 Do Quit
|
---|
280 | . . Set SIGS=$Select(X<0:-1,1:1)
|
---|
281 | . . Set VALUE=1/(1/X/X-1)
|
---|
282 | . . Set X=$%SQRT^MATH(VALUE)
|
---|
283 | . . Set VALUE=$$ARCTAN(X,10)
|
---|
284 | . . Set VALUE=VALUE*SIGS
|
---|
285 | . . Quit
|
---|
286 | . Set (VALUE,L)=X
|
---|
287 | . Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
288 | . For K=3:2 Do Quit:($Translate(L,"-")<LIM)
|
---|
289 | . . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
|
---|
290 | . . Quit
|
---|
291 | . Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
|
---|
292 | . Quit
|
---|
293 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
294 | For K=3:2 Do Quit:($Translate(1/L,"-")<LIM)
|
---|
295 | . ;
|
---|
296 | . ;;;" Set L=L*X*X,VALUE=VALUE+(1/(K*L)*SIGN), ;" Number ~~
|
---|
297 | . ;" Eli Reidler (28 June 1996)
|
---|
298 | . Set L=L*X*X,VALUE=VALUE+(1/(K*L)*SIGN)
|
---|
299 | . ;;;
|
---|
300 | . ;
|
---|
301 | . Set SIGN=SIGN*-1
|
---|
302 | . Quit
|
---|
303 | Set VALUE=$Select(FOLD:$%PI^MATH()/2-VALUE,1:VALUE)
|
---|
304 | Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
|
---|
305 | Quit VALUE
|
---|
306 | ;===
|
---|
307 | ;
|
---|
308 | ;
|
---|
309 | ARCTANH(X,PREC) ;
|
---|
310 | If X<-1 Set $Ecode=",M28,"
|
---|
311 | If X>1 Set $Ecode=",M28,"
|
---|
312 | ;
|
---|
313 | ;;;" ;" Number ~~
|
---|
314 | ;" Winfried Gerum (8 June 1995)
|
---|
315 | ;" Alan Frank (October 1995)
|
---|
316 | Set PREC=$Get(PREC,11)
|
---|
317 | ;;;
|
---|
318 | ;
|
---|
319 | Quit $%LOG^MATH(1+X/(1-X),PREC)/2
|
---|
320 | ;===
|
---|
321 | ;
|
---|
322 | ;
|
---|
323 | CABS(Z) ;
|
---|
324 | New ZRE,ZIM
|
---|
325 | Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
|
---|
326 | Quit $%SQRT^MATH(ZRE*ZRE+(ZIM*ZIM))
|
---|
327 | ;===
|
---|
328 | ;
|
---|
329 | ;
|
---|
330 | CADD(X,Y) ;
|
---|
331 | New XRE,XIM,YRE,YIM
|
---|
332 | Set XRE=+X,XIM=+$Piece(X,"%",2)
|
---|
333 | Set YRE=+Y,YIM=+$Piece(Y,"%",2)
|
---|
334 | Quit XRE+YRE_"%"_(XIM+YIM)
|
---|
335 | ;===
|
---|
336 | ;
|
---|
337 | ;
|
---|
338 | CCOS(Z,PREC) ;
|
---|
339 | New E1,E2,IA
|
---|
340 | ;
|
---|
341 | ;;;" ;" Number ~~
|
---|
342 | ;" Alan Frank (October 1995)
|
---|
343 | Set PREC=$Get(PREC,11)
|
---|
344 | ;;;
|
---|
345 | ;
|
---|
346 | Set IA=$%CMUL^MATH(Z,"0%1")
|
---|
347 | Set E1=$%CEXP^MATH(IA,PREC)
|
---|
348 | Set IA=-IA_"%"_(-$Piece(IA,"%",2))
|
---|
349 | Set E2=$%CEXP^MATH(IA,PREC)
|
---|
350 | Set IA=$%CADD^MATH(E1,E2)
|
---|
351 | Quit $%CMUL^MATH(IA,"0.5%0")
|
---|
352 | ;===
|
---|
353 | ;
|
---|
354 | ;
|
---|
355 | CDIV(X,Y) ;
|
---|
356 | New D,IM,RE,XIM,XRE,YIM,YRE
|
---|
357 | Set XRE=+X,XIM=+$Piece(X,"%",2)
|
---|
358 | Set YRE=+Y,YIM=+$Piece(Y,"%",2)
|
---|
359 | Set D=YRE*YRE+(YIM*YIM)
|
---|
360 | Set RE=XRE*YRE+(XIM*YIM)/D
|
---|
361 | Set IM=XIM*YRE-(XRE*YIM)/D
|
---|
362 | Quit RE_"%"_IM
|
---|
363 | ;===
|
---|
364 | ;
|
---|
365 | ;
|
---|
366 | CEXP(Z,PREC) ;
|
---|
367 | New R,ZIM,ZRE
|
---|
368 | ;
|
---|
369 | ;;;" ;" Number ~~
|
---|
370 | ;" Alan Frank (October 1995)
|
---|
371 | Set PREC=$Get(PREC,11)
|
---|
372 | ;;;
|
---|
373 | ;
|
---|
374 | Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
|
---|
375 | Set R=$%EXP^MATH(ZRE,PREC)
|
---|
376 | Quit R*$%COS^MATH(ZIM,PREC)_"%"_(R*$%SIN^MATH(ZIM,PREC))
|
---|
377 | ;===
|
---|
378 | ;
|
---|
379 | ;
|
---|
380 | CLOG(Z,PREC) ;
|
---|
381 | New ABS,ARG,ZIM,ZRE
|
---|
382 | ;
|
---|
383 | ;;;" ;" Number ~~
|
---|
384 | ;" Alan Frank (October 1995)
|
---|
385 | Set PREC=$Get(PREC,11)
|
---|
386 | ;;;
|
---|
387 | ;
|
---|
388 | Set ABS=$%CABS^MATH(Z)
|
---|
389 | Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
|
---|
390 | ;
|
---|
391 | ;;;" Set ARG=$%ARCTAN^MATH(ZIM,ZRE,PREC) ;" Number ~~
|
---|
392 | ;" Alan Frank (October 1995)
|
---|
393 | Set ARG=$%ARCTAN^MATH(ZIM/ZRE,PREC)
|
---|
394 | ;;;
|
---|
395 | ;
|
---|
396 | Quit $%LOG^MATH(ABS,PREC)_"%"_ARG
|
---|
397 | ;===
|
---|
398 | ;
|
---|
399 | ;
|
---|
400 | CMUL(X,Y) ;
|
---|
401 | New XIM,XRE,YIM,YRE
|
---|
402 | Set XRE=+X,XIM=+$Piece(X,"%",2)
|
---|
403 | Set YRE=+Y,YIM=+$Piece(Y,"%",2)
|
---|
404 | Quit XRE*YRE-(XIM*YIM)_"%"_(XRE*YIM+(XIM*YRE))
|
---|
405 | ;===
|
---|
406 | ;
|
---|
407 | ;
|
---|
408 | COMPLEX(X) Quit +X_"%0"
|
---|
409 | ;===
|
---|
410 | ;
|
---|
411 | ;
|
---|
412 | CONJUG(Z) ;
|
---|
413 | New ZIM,ZRE
|
---|
414 | Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
|
---|
415 | Quit ZRE_"%"_(-ZIM)
|
---|
416 | ;===
|
---|
417 | ;
|
---|
418 | ;
|
---|
419 | COS(X,PREC) ;
|
---|
420 | New L,LIM,K,SIGN,VALUE
|
---|
421 | ;
|
---|
422 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
423 | ;" Winfried Gerum (8 June 1995)
|
---|
424 | ;" Comment: The official description does not mention than
|
---|
425 | ;" the function may also be called with the first
|
---|
426 | ;" parameter in degrees, minutes and seconds.
|
---|
427 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
428 | ;;;
|
---|
429 | ;
|
---|
430 | Set PREC=$Get(PREC,11)
|
---|
431 | Set X=X#(2*$%PI^MATH())
|
---|
432 | Set (VALUE,L)=1,SIGN=-1
|
---|
433 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
434 | For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
435 | . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
|
---|
436 | . Quit
|
---|
437 | Quit VALUE
|
---|
438 | ;===
|
---|
439 | ;
|
---|
440 | ;
|
---|
441 | COS(X) ;
|
---|
442 | ;;;" ;" Number ~~
|
---|
443 | ;" Winfried Gerum (8 June 1995)
|
---|
444 | ;" Comment: This version of the function is
|
---|
445 | ;" optimized for speed, not for precision.
|
---|
446 | ;" The 'precision' parameter is not supported,
|
---|
447 | ;" and the precision is at best 1 in 10**-9.
|
---|
448 | ;" Note that this function does not accept its
|
---|
449 | ;" parameter in degrees, minutes and seconds.
|
---|
450 | ;;;
|
---|
451 | ;
|
---|
452 | New A,N,PI,R,SIGN,XX
|
---|
453 | ;
|
---|
454 | ;" This approximation only works for 0 <= x <= pi/2
|
---|
455 | ;" so reduce angle to correct quadrant.
|
---|
456 | ;
|
---|
457 | Set PI=$%PI^MATH(),X=X#(PI*2),SIGN=1
|
---|
458 | Set:X>PI X=2*PI-X
|
---|
459 | Set:X*2>PI X=PI-X,SIGN=-1
|
---|
460 | ;
|
---|
461 | Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
|
---|
462 | Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
|
---|
463 | Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
|
---|
464 | Quit R*SIGN
|
---|
465 | ;===
|
---|
466 | ;
|
---|
467 | ;
|
---|
468 | COSH(X,PREC) ;
|
---|
469 | ;
|
---|
470 | ;;;" New F,I,P,R,T,XX ;" Number ~~
|
---|
471 | ;" Winfried Gerum (8 June 1995)
|
---|
472 | New E,F,I,P,R,T,XX
|
---|
473 | ;;;
|
---|
474 | ;
|
---|
475 | Set PREC=$Get(PREC,11)+1
|
---|
476 | Set @("E=1E-"_PREC)
|
---|
477 | Set XX=X*X,F=1,(P,R,T)=1,I=1
|
---|
478 | 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
|
---|
479 | Quit R
|
---|
480 | ;===
|
---|
481 | ;
|
---|
482 | ;
|
---|
483 | COT(X,PREC) ;
|
---|
484 | New C,L,LIM,K,SIGN,VALUE
|
---|
485 | ;
|
---|
486 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
487 | ;" Winfried Gerum (8 June 1995)
|
---|
488 | ;" Comment: The official description does not mention than
|
---|
489 | ;" the function may also be called with the first
|
---|
490 | ;" parameter in degrees, minutes and seconds.
|
---|
491 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
492 | ;;;
|
---|
493 | ;
|
---|
494 | Set PREC=$Get(PREC,11)
|
---|
495 | Set (VALUE,L)=1,SIGN=-1
|
---|
496 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
497 | For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
498 | . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
|
---|
499 | . Quit
|
---|
500 | Set C=VALUE
|
---|
501 | Set X=X#(2*$%PI^MATH())
|
---|
502 | Set (VALUE,L)=X,SIGN=-1
|
---|
503 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
504 | For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
505 | . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
|
---|
506 | . Quit
|
---|
507 | If 'VALUE Quit "INFINITE"
|
---|
508 | Quit VALUE=C/VALUE
|
---|
509 | ;===
|
---|
510 | ;
|
---|
511 | ;
|
---|
512 | COTH(X,PREC) ;
|
---|
513 | New SINH
|
---|
514 | If 'X Quit "INFINITE"
|
---|
515 | ;
|
---|
516 | ;;;" ;" Number ~~
|
---|
517 | ;" Winfried Gerum (8 June 1995)
|
---|
518 | ;" Alan Frank (October 1995)
|
---|
519 | Set PREC=$Get(PREC,11)
|
---|
520 | ;;;
|
---|
521 | ;
|
---|
522 | Set SINH=$%SINH^MATH(X,PREC)
|
---|
523 | If 'SINH Quit "INFINITE"
|
---|
524 | Quit $%COSH^MATH(X,PREC)/SINH
|
---|
525 | ;===
|
---|
526 | ;
|
---|
527 | ;
|
---|
528 | CPOWER(Z,N,PREC) ;
|
---|
529 | New AR,NIM,NRE,PHI,PI,R,RHO,TH,ZIM,ZRE
|
---|
530 | ;
|
---|
531 | ;;;" ;" Number ~~
|
---|
532 | ;" Alan Frank (October 1995)
|
---|
533 | Set PREC=$Get(PREC,11)
|
---|
534 | ;;;
|
---|
535 | ;
|
---|
536 | Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
|
---|
537 | Set NRE=+N,NIM=+$Piece(N,"%",2)
|
---|
538 | If 'ZRE,'ZIM,'NRE,'NIM Set $Ecode=",M28,"
|
---|
539 | ;
|
---|
540 | ;;;" If 'ZRE,'ZIM Quit "0%0% ;" Number ~~
|
---|
541 | ;" Eli Reidler (28 June 1996)
|
---|
542 | If 'ZRE,'ZIM Quit "0%0"
|
---|
543 | ;;;
|
---|
544 | ;
|
---|
545 | Set PI=$%PI^MATH()
|
---|
546 | ;
|
---|
547 | ;;;" Set R=$%SQRT^MATH(ZRE*ZRE+(ZIM*ZIM,PREC)) ;" Number ~~
|
---|
548 | ;" Winfried Gerum (8 June 1995)
|
---|
549 | ;" Eli Reidler (28 June 1996)
|
---|
550 | Set R=$%SQRT^MATH(ZRE*ZRE+(ZIM*ZIM),PREC)
|
---|
551 | ;;;
|
---|
552 | ;
|
---|
553 | ;
|
---|
554 | ;;;" If ZRE Set TH=$%ARCTAN^MATH(ZIM,ZRE,PREC) ;" Number ~~
|
---|
555 | ;" Alan Frank (October 1995)
|
---|
556 | If ZRE Set TH=$%ARCTAN^MATH(ZIM/ZRE,PREC)
|
---|
557 | ;;;
|
---|
558 | ;
|
---|
559 | ;;;" Else Set TH=$SELECT(ZRE>0:PI/2,1:-PI/2) ;" Number ~~
|
---|
560 | ;" Winfried Gerum (8 June 1995)
|
---|
561 | Else Set TH=$SELECT(ZIM>0:PI/2,1:-PI/2)
|
---|
562 | ;;;
|
---|
563 | ;
|
---|
564 | Set RHO=$%LOG^MATH(R,PREC)
|
---|
565 | Set AR=$%EXP^MATH(RHO*NRE-(TH*NIM),PREC)
|
---|
566 | Set PHI=RHO*NIM+(NRE*TH)
|
---|
567 | Quit AR*$%COS^MATH(PHI,PREC)_"%"_(AR*$%SIN^MATH(PHI,PREC))
|
---|
568 | ;===
|
---|
569 | ;
|
---|
570 | ;
|
---|
571 | CSC(X,PREC) ;
|
---|
572 | New L,LIM,K,SIGN,VALUE
|
---|
573 | ;
|
---|
574 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
575 | ;" Winfried Gerum (8 June 1995)
|
---|
576 | ;" Comment: The official description does not mention than
|
---|
577 | ;" the function may also be called with the first
|
---|
578 | ;" parameter in degrees, minutes and seconds.
|
---|
579 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
580 | ;;;
|
---|
581 | ;
|
---|
582 | ;;;" Set PREC=$Select($Data(PREC)#2:PREC,1:10) ;" Number ~~
|
---|
583 | ;" Winfried Gerum (8 June 1995)
|
---|
584 | Set PREC=$Get(PREC,11)
|
---|
585 | ;;;
|
---|
586 | ;
|
---|
587 | Set X=X#(2*$%PI^MATH())
|
---|
588 | Set (VALUE,L)=X,SIGN=-1
|
---|
589 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
590 | For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
591 | . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
|
---|
592 | . Quit
|
---|
593 | If 'VALUE Quit "INFINITE"
|
---|
594 | Quit 1/VALUE
|
---|
595 | ;===
|
---|
596 | ;
|
---|
597 | ;
|
---|
598 | ;
|
---|
599 | CSCH(X,PREC) ;;;Quit 1/$%SINH^MATH(X,PREC) ;" Number ~~
|
---|
600 | ;" Winfried Gerum (8 June 1995)
|
---|
601 | ;" Alan Frank (October 1995)
|
---|
602 | Quit 1/$%SINH^MATH(X,$Get(PREC,11))
|
---|
603 | ;;;
|
---|
604 | ;
|
---|
605 | ;===
|
---|
606 | ;
|
---|
607 | ;
|
---|
608 | CSIN(Z,PREC) ;
|
---|
609 | New IA,E1,E2
|
---|
610 | ;
|
---|
611 | ;;;" ;" Number ~~
|
---|
612 | ;" Alan Frank (October 1995)
|
---|
613 | Set PREC=$Get(PREC,11)
|
---|
614 | ;;;
|
---|
615 | ;
|
---|
616 | Set IA=$%CMUL^MATH(Z,"0%1")
|
---|
617 | Set E1=$%CEXP^MATH(IA,PREC)
|
---|
618 | Set IA=-IA_"%"_(-$Piece(IA,"%",2))
|
---|
619 | Set E2=$%CEXP^MATH(IA,PREC)
|
---|
620 | Set IA=$%CSUB^MATH(E1,E2)
|
---|
621 | Set IA=$%CMUL^MATH(IA,"0.5%0")
|
---|
622 | Quit $%CMUL^MATH("0%-1",IA)
|
---|
623 | ;===
|
---|
624 | ;
|
---|
625 | ;
|
---|
626 | CSUB(X,Y) ;
|
---|
627 | New XIM,XRE,YIM,YRE
|
---|
628 | Set XRE=+X,XIM=+$Piece(X,"%",2)
|
---|
629 | Set YRE=+Y,YIM=+$Piece(Y,"%",2)
|
---|
630 | Quit XRE-YRE_"%"_(XIM-YIM)
|
---|
631 | ;===
|
---|
632 | ;
|
---|
633 | ;
|
---|
634 | DECDMS(X,PREC) ;
|
---|
635 | Set PREC=$Get(PREC,5)
|
---|
636 | Set X=X#360*3600
|
---|
637 | Set X=+$Justify(X,0,$Select((PREC-$Length(X\1))'<0:PREC-$Length(X\1),1:0))
|
---|
638 | Quit X\3600_":"_(X\60#60)_":"_(X#60)
|
---|
639 | ;===
|
---|
640 | ;
|
---|
641 | ;
|
---|
642 | DEGRAD(X) Quit X*3.14159265358979/180
|
---|
643 | ;===
|
---|
644 | ;
|
---|
645 | ;
|
---|
646 | DMSDEC(X) ;
|
---|
647 | Quit $Piece(X,":")+($Piece(X,":",2)/60)+($Piece(X,":",3)/3600)
|
---|
648 | ;===
|
---|
649 | ;
|
---|
650 | ;
|
---|
651 | E() Quit 2.71828182845905
|
---|
652 | ;===
|
---|
653 | ;
|
---|
654 | ;
|
---|
655 | EXP(X,PREC) ;
|
---|
656 | New L,LIM,K,VALUE
|
---|
657 | Set PREC=$Get(PREC,11)
|
---|
658 | Set L=X,VALUE=X+1
|
---|
659 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
660 | For K=2:1 Set L=L*X/K,VALUE=VALUE+L Quit:($Translate(L,"-")<LIM)
|
---|
661 | Quit VALUE
|
---|
662 | ;===
|
---|
663 | ;
|
---|
664 | ;
|
---|
665 | LOG(X,PREC) ;
|
---|
666 | New L,LIM,M,N,K,VALUE
|
---|
667 | If X'>0 Set $Ecode=",M28,"
|
---|
668 | Set PREC=$Get(PREC,11)
|
---|
669 | Set M=1
|
---|
670 | ;
|
---|
671 | ;;;" If X>0 For N=0:1 Quit:(X/M)<10 Set M=M*10 ;" Number ~~
|
---|
672 | ;" Winfried Gerum (8 June 1995)
|
---|
673 | For N=0:1 Quit:(X/M)<10 Set M=M*10
|
---|
674 | ;;;
|
---|
675 | ;
|
---|
676 | If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
|
---|
677 | Set X=X/M
|
---|
678 | Set X=(X-1)/(X+1),(VALUE,L)=X
|
---|
679 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
680 | For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
|
---|
681 | Set VALUE=VALUE*2+(N*2.30258509298749)
|
---|
682 | Quit VALUE
|
---|
683 | ;===
|
---|
684 | ;
|
---|
685 | ;
|
---|
686 | LOG10(X,PREC) ;
|
---|
687 | New L,LIM,M,N,K,VALUE
|
---|
688 | If X'>0 Set $Ecode=",M28,"
|
---|
689 | Set PREC=$Get(PREC,11)
|
---|
690 | Set M=1
|
---|
691 | ;
|
---|
692 | ;;;" If X>0 For N=0:1 Quit:(X/M)<10 Set M=M*10 ;" Number ~~
|
---|
693 | ;" Winfried Gerum (8 June 1995)
|
---|
694 | For N=0:1 Quit:(X/M)<10 Set M=M*10
|
---|
695 | ;;;
|
---|
696 | ;
|
---|
697 | If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
|
---|
698 | Set X=X/M
|
---|
699 | Set X=(X-1)/(X+1),(VALUE,L)=X
|
---|
700 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
701 | For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
|
---|
702 | Set VALUE=VALUE*2+(N*2.30258509298749)
|
---|
703 | Quit VALUE/2.30258509298749
|
---|
704 | ;===
|
---|
705 | ;
|
---|
706 | ;
|
---|
707 | MTXADD(A,B,R,ROWS,COLS) ;
|
---|
708 | ;" Add A[ROWS,COLS] to B[ROWS,COLS],
|
---|
709 | ;" result goes to R[ROWS,COLS]
|
---|
710 | IF $DATA(A)<10 QUIT 0
|
---|
711 | IF $DATA(B)<10 QUIT 0
|
---|
712 | IF $GET(ROWS)<1 QUIT 0
|
---|
713 | IF $GET(COLS)<1 QUIT 0
|
---|
714 | ;
|
---|
715 | NEW ROW,COL,ANY
|
---|
716 | FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
|
---|
717 | . KVALUE R(ROW,COL) SET ANY=0
|
---|
718 | . SET:$DATA(A(ROW,COL))#2 ANY=1
|
---|
719 | . SET:$DATA(B(ROW,COL))#2 ANY=1
|
---|
720 | . SET:ANY R(ROW,COL)=$GET(A(ROW,COL))+$GET(B(ROW,COL))
|
---|
721 | . QUIT
|
---|
722 | QUIT 1
|
---|
723 | ;===
|
---|
724 | ;
|
---|
725 | ;
|
---|
726 | MTXCOF(A,I,K,N) ;
|
---|
727 | ;" Compute cofactor for element [i,k]
|
---|
728 | ;" in matrix A[N,N]
|
---|
729 | NEW T,R,C,RR,CC
|
---|
730 | SET CC=0 FOR C=1:1:N DO:C'=K
|
---|
731 | . SET CC=CC+1,RR=0
|
---|
732 | . FOR R=1:1:N SET:R'=I RR=RR+1,T(RR,CC)=$GET(A(R,C))
|
---|
733 | . QUIT
|
---|
734 | QUIT $%MTXDET^MATH(.T,N-1)
|
---|
735 | ;===
|
---|
736 | ;
|
---|
737 | ;
|
---|
738 | MTXCOPY(A,R,ROWS,COLS) ;
|
---|
739 | ;" Copy A[ROWS,COLS] to R[ROWS,COLS]
|
---|
740 | IF $DATA(A)<10 QUIT 0
|
---|
741 | IF $GET(ROWS)<1 QUIT 0
|
---|
742 | IF $GET(COLS)<1 QUIT 0
|
---|
743 | ;
|
---|
744 | NEW ROW,COL
|
---|
745 | FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
|
---|
746 | . KVALUE R(ROW,COL)
|
---|
747 | . SET:$DATA(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)
|
---|
748 | . QUIT
|
---|
749 | QUIT 1
|
---|
750 | ;===
|
---|
751 | ;
|
---|
752 | ;
|
---|
753 | MTXDET(A,N) ;
|
---|
754 | ;" Compute determinant of matrix A[N,N]
|
---|
755 | IF $DATA(A)<10 QUIT ""
|
---|
756 | IF $GET(N)<1 QUIT ""
|
---|
757 | ;
|
---|
758 | ;" First the simple cases
|
---|
759 | ;
|
---|
760 | IF N=1 QUIT $GET(A(1,1))
|
---|
761 | IF N=2 QUIT $GET(A(1,1))*$GET(A(2,2))-($GET(A(1,2))*$GET(A(2,1)))
|
---|
762 | ;
|
---|
763 | NEW DET,I,SIGN
|
---|
764 | ;
|
---|
765 | ;" Det A = sum (k=1:n) element (i,k) * cofactor [i,k]
|
---|
766 | ;
|
---|
767 | SET DET=0,SIGN=1
|
---|
768 | FOR I=1:1:N DO
|
---|
769 | . SET DET=$GET(A(1,I))*$%MTXCOF^MATH(.A,1,I,N)*SIGN+DET
|
---|
770 | . SET SIGN=-SIGN
|
---|
771 | . QUIT
|
---|
772 | QUIT DET
|
---|
773 | ;===
|
---|
774 | ;
|
---|
775 | ;
|
---|
776 | MTXEQU(A,B,R,N,M) ;
|
---|
777 | ;" Solve matrix equation A [M,M] * R [M,N] = B [M,N]
|
---|
778 | IF $GET(M)<1 QUIT ""
|
---|
779 | IF $GET(N)<1 QUIT ""
|
---|
780 | ;;;IF '$%MTXDET^MATH(.A) QUIT 0
|
---|
781 | ;" Ed de Moel, 29 August 1999
|
---|
782 | IF '$%MTXDET^MATH(.A,M) QUIT 0
|
---|
783 | ;;;
|
---|
784 | ;
|
---|
785 | NEW I,I1,J,J1,J2,K,L,T,T1,T2,TEMP,X
|
---|
786 | ;
|
---|
787 | SET X=$%MTXCOPY^MATH(.A,.T,N,N)
|
---|
788 | SET X=$%MTXCOPY^MATH(.B,.R,N,M)
|
---|
789 | ;
|
---|
790 | ;" Reduction of matrix A
|
---|
791 | ;" Steps of reduction are counted by index K
|
---|
792 | ;
|
---|
793 | FOR K=1:1:N-1 DO
|
---|
794 | . ;
|
---|
795 | . ;" Search for largest coefficient of T
|
---|
796 | . ;" (denoted by TEMP)
|
---|
797 | . ;" in first column of reduced system
|
---|
798 | . ;
|
---|
799 | . SET TEMP=0,J2=K
|
---|
800 | . FOR J1=K:1:N DO
|
---|
801 | . . QUIT:$TRANSLATE($GET(T(J1,K)),"-")>$TRANSLATE(TEMP,"-")
|
---|
802 | . . SET TEMP=T(J1,K),J2=J1
|
---|
803 | . . QUIT
|
---|
804 | . ;
|
---|
805 | . ;" Exchange row number K with row number J2,
|
---|
806 | . ;" if necessary
|
---|
807 | . ;
|
---|
808 | . DO:J2'=K
|
---|
809 | . . ;
|
---|
810 | . . FOR J=K:1:N DO
|
---|
811 | . . . SET T1=$GET(T(K,J)),T2=$GET(T(J2,J))
|
---|
812 | . . . KILL T(K,J),T(J2,J)
|
---|
813 | . . . IF T1'="" SET T(J2,J)=T1
|
---|
814 | . . . IF T2'="" SET T(K,J)=T2
|
---|
815 | . . . QUIT
|
---|
816 | . . FOR J=1:1:M DO
|
---|
817 | . . . SET T1=$GET(R(K,J)),T2=$GET(R(J2,J))
|
---|
818 | . . . KILL R(K,J),R(J2,J)
|
---|
819 | . . . IF T1'="" SET R(J2,J)=T1
|
---|
820 | . . . IF T2'="" SET R(K,J)=T2
|
---|
821 | . . . QUIT
|
---|
822 | . . QUIT
|
---|
823 | . ;
|
---|
824 | . ;" Actual reduction
|
---|
825 | . ;
|
---|
826 | . FOR I=K+1:1:N DO
|
---|
827 | . . FOR J=K+1:1:N DO
|
---|
828 | . . . QUIT:'$GET(T(K,K))
|
---|
829 | . . . SET T(I,J)=-$GET(T(K,J))*$GET(T(I,K))/T(K,K)+$GET(T(I,J))
|
---|
830 | . . . QUIT
|
---|
831 | . . FOR J=1:1:M DO
|
---|
832 | . . . QUIT:'$GET(T(K,K))
|
---|
833 | . . . SET R(I,J)=-$GET(R(K,J))*$GET(T(I,K))/T(K,K)+$GET(R(I,J))
|
---|
834 | . . . QUIT
|
---|
835 | . . QUIT
|
---|
836 | . QUIT
|
---|
837 | ;
|
---|
838 | ;" Backsubstitution
|
---|
839 | ;
|
---|
840 | FOR J=1:1:M DO
|
---|
841 | . IF $GET(T(N,N)) SET R(N,J)=$GET(R(N,J))/T(N,N)
|
---|
842 | . IF N-1>0 FOR I1=1:1:N-1 DO
|
---|
843 | . . SET I=N-I1
|
---|
844 | . . FOR L=I+1:1:N DO
|
---|
845 | . . . SET R(I,J)=-$GET(T(I,L))*$GET(R(L,J))+$GET(R(I,J))
|
---|
846 | . . . QUIT
|
---|
847 | . . IF $GET(T(I,I)) SET R(I,J)=$GET(R(I,J))/$GET(T(I,I))
|
---|
848 | . . QUIT
|
---|
849 | . QUIT
|
---|
850 | ;;;QUIT $%MTXDET^MATH(.R)
|
---|
851 | ;" Ed de Moel, 29 Aug 1999
|
---|
852 | QUIT $SELECT(M=N:$%MTXDET^MATH(.R,M),1:1)
|
---|
853 | ;;;
|
---|
854 | ;===
|
---|
855 | ;
|
---|
856 | MTXINV(A,R,N) ;
|
---|
857 | ;" Invert A[N,N], result goes to R[N,N]
|
---|
858 | IF $DATA(A)<10 QUIT 0
|
---|
859 | IF $GET(N)<1 QUIT 0
|
---|
860 | ;
|
---|
861 | NEW T,X
|
---|
862 | SET X=$%MTXUNIT^MATH(.T,N)
|
---|
863 | QUIT $%MTXEQU^MATH(.A,.T,.R,N,N)
|
---|
864 | ;===
|
---|
865 | ;
|
---|
866 | ;
|
---|
867 | MTXMUL(A,B,R,M,L,N) ;
|
---|
868 | ;" Multiply A[M,L] by B[L,N], result goes to R[M,N]
|
---|
869 | IF $DATA(A)<10 QUIT 0
|
---|
870 | IF $DATA(B)<10 QUIT 0
|
---|
871 | IF $GET(L)<1 QUIT 0
|
---|
872 | IF $GET(M)<1 QUIT 0
|
---|
873 | IF $GET(N)<1 QUIT 0
|
---|
874 | ;
|
---|
875 | NEW I,J,K,SUM,ANY
|
---|
876 | FOR I=1:1:M FOR J=1:1:N DO
|
---|
877 | . SET (SUM,ANY)=0
|
---|
878 | . KVALUE R(I,J)
|
---|
879 | . FOR K=1:1:L DO
|
---|
880 | . . SET:$DATA(A(I,K))#2 ANY=1
|
---|
881 | . . SET:$DATA(B(K,J))#2 ANY=1
|
---|
882 | . . SET SUM=$GET(A(I,K))*$GET(B(K,J))+SUM
|
---|
883 | . . QUIT
|
---|
884 | . SET:ANY R(I,J)=SUM
|
---|
885 | . QUIT
|
---|
886 | QUIT 1
|
---|
887 | ;===
|
---|
888 | ;
|
---|
889 | ;
|
---|
890 | MTXSCA(A,R,ROWS,COLS,S) ;
|
---|
891 | ;" Multiply A[ROWS,COLS] with the scalar S,
|
---|
892 | ;" result goes to R[ROWS,COLS]
|
---|
893 | IF $DATA(A)<10 QUIT 0
|
---|
894 | IF $GET(ROWS)<1 QUIT 0
|
---|
895 | IF $GET(COLS)<1 QUIT 0
|
---|
896 | IF '($DATA(S)#2) QUIT 0
|
---|
897 | ;
|
---|
898 | NEW ROW,COL
|
---|
899 | FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
|
---|
900 | . KVALUE R(ROW,COL)
|
---|
901 | . SET:$DATA(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)*S
|
---|
902 | . QUIT
|
---|
903 | QUIT 1
|
---|
904 | ;===
|
---|
905 | ;
|
---|
906 | ;
|
---|
907 | MTXSUB(A,B,R,ROWS,COLS) ;
|
---|
908 | ;" Subtract B[ROWS,COLS] from A[ROWS,COLS],
|
---|
909 | ;" result goes to R[ROWS,COLS]
|
---|
910 | IF $DATA(A)<10 QUIT 0
|
---|
911 | IF $DATA(B)<10 QUIT 0
|
---|
912 | IF $GET(ROWS)<1 QUIT 0
|
---|
913 | IF $GET(COLS)<1 QUIT 0
|
---|
914 | ;
|
---|
915 | NEW ROW,COL,ANY
|
---|
916 | FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
|
---|
917 | . KVALUE R(ROW,COL) SET ANY=0
|
---|
918 | . SET:$DATA(A(ROW,COL))#2 ANY=1
|
---|
919 | . SET:$DATA(B(ROW,COL))#2 ANY=1
|
---|
920 | . ;
|
---|
921 | . ;;;" SET:ANY R(ROW,COL)=$GET(A(ROW,COL)-$GET(B(ROW,COL)) ;" Number ~~
|
---|
922 | . ;" Eli Reidler (28 June 1996)
|
---|
923 | . SET:ANY R(ROW,COL)=$GET(A(ROW,COL))-$GET(B(ROW,COL))
|
---|
924 | . ;;;
|
---|
925 | . ;
|
---|
926 | . QUIT
|
---|
927 | QUIT 1
|
---|
928 | ;===
|
---|
929 | ;
|
---|
930 | ;
|
---|
931 | MTXTRP(A,R,M,N) ;
|
---|
932 | ;" Transpose A[M,N], result goes to R[N,M]
|
---|
933 | IF $DATA(A)<10 QUIT 0
|
---|
934 | IF $GET(M)<1 QUIT 0
|
---|
935 | IF $GET(N)<1 QUIT 0
|
---|
936 | ;
|
---|
937 | NEW I,J,K,D1,V1,D2,V2
|
---|
938 | FOR I=1:1:M+N-1 FOR J=1:1:I+1\2 DO
|
---|
939 | . SET K=I-J+1
|
---|
940 | . IF K=J DO QUIT
|
---|
941 | . . SET V1=$GET(A(J,J)),D1=$DATA(A(J,J))#2
|
---|
942 | . . IF J'>N,J'>M KVALUE R(J,J) SET:D1 R(J,J)=V1
|
---|
943 | . . QUIT
|
---|
944 | . ;
|
---|
945 | . SET V1=$GET(A(K,J)),D1=$DATA(A(K,J))#2
|
---|
946 | . SET V2=$GET(A(J,K)),D2=$DATA(A(J,K))#2
|
---|
947 | . IF K'>M,J'>N KVALUE R(K,J) SET:D2 R(K,J)=V2
|
---|
948 | . IF J'>M,K'>N KVALUE R(J,K) SET:D1 R(J,K)=V1
|
---|
949 | . QUIT
|
---|
950 | QUIT 1
|
---|
951 | ;===
|
---|
952 | ;
|
---|
953 | ;
|
---|
954 | MTXUNIT(R,N,SPARSE) ;
|
---|
955 | ;" Create a unit matrix R[N,N]
|
---|
956 | IF $GET(N)<1 QUIT 0
|
---|
957 | ;
|
---|
958 | NEW ROW,COL
|
---|
959 | FOR ROW=1:1:N FOR COL=1:1:N DO
|
---|
960 | . KVALUE R(ROW,COL)
|
---|
961 | . IF $GET(SPARSE) QUIT:ROW'=COL
|
---|
962 | . SET R(ROW,COL)=$SELECT(ROW=COL:1,1:0)
|
---|
963 | . QUIT
|
---|
964 | QUIT 1
|
---|
965 | ;===
|
---|
966 | ;
|
---|
967 | ;
|
---|
968 | PI() Quit 3.14159265358979
|
---|
969 | ;===
|
---|
970 | ;
|
---|
971 | ;
|
---|
972 | PRODUCE(IN,SPEC,MAX) ;
|
---|
973 | NEW VALUE,AGAIN,P1,P2,I,COUNT
|
---|
974 | SET VALUE=IN,COUNT=0
|
---|
975 | FOR DO QUIT:'AGAIN
|
---|
976 | . SET AGAIN=0
|
---|
977 | . SET I=""
|
---|
978 | . FOR SET I=$ORDER(SPEC(I)) QUIT:I="" DO QUIT:COUNT<0
|
---|
979 | . . QUIT:$GET(SPEC(I,1))=""
|
---|
980 | . . QUIT:'($DATA(SPEC(I,2))#2)
|
---|
981 | . . FOR QUIT:VALUE'[SPEC(I,1) DO QUIT:COUNT<0
|
---|
982 | . . . SET P1=$PIECE(VALUE,SPEC(I,1),1)
|
---|
983 | . . . SET P2=$PIECE(VALUE,SPEC(I,1),2,$LENGTH(VALUE))
|
---|
984 | . . . SET VALUE=P1_SPEC(I,2)_P2,AGAIN=1
|
---|
985 | . . . SET COUNT=COUNT+1
|
---|
986 | . . . IF $DATA(MAX),COUNT>MAX SET COUNT=-1,AGAIN=0
|
---|
987 | . . . QUIT
|
---|
988 | . . QUIT
|
---|
989 | . QUIT
|
---|
990 | QUIT VALUE
|
---|
991 | ;===
|
---|
992 | ;
|
---|
993 | ;
|
---|
994 | RADDEG(X) Quit X*180/3.14159265358979
|
---|
995 | ;===
|
---|
996 | ;
|
---|
997 | ;
|
---|
998 | REPLACE(IN,SPEC) ;
|
---|
999 | NEW L,MASK,K,I,LT,F,VALUE
|
---|
1000 | SET L=$LENGTH(IN),MASK=$JUSTIFY("",L)
|
---|
1001 | SET I="" FOR SET I=$ORDER(SPEC(I)) QUIT:I="" DO
|
---|
1002 | . QUIT:'($DATA(SPEC(I,1))#2)
|
---|
1003 | . QUIT:SPEC(I,1)=""
|
---|
1004 | . QUIT:'($DATA(SPEC(I,2))#2)
|
---|
1005 | . SET LT=$LENGTH(SPEC(I,1))
|
---|
1006 | . SET F=0 FOR SET F=$FIND(IN,SPEC(I,1),F) QUIT:F<1 DO
|
---|
1007 | . . QUIT:$EXTRACT(MASK,F-LT,F-1)["X"
|
---|
1008 | . . SET VALUE(F-LT)=SPEC(I,2)
|
---|
1009 | . . SET $EXTRACT(MASK,F-LT,F-1)=$TRANSLATE($JUSTIFY("",LT)," ","X")
|
---|
1010 | . . QUIT
|
---|
1011 | . QUIT
|
---|
1012 | SET VALUE="" FOR K=1:1:L DO
|
---|
1013 | . IF $EXTRACT(MASK,K)=" " SET VALUE=VALUE_$EXTRACT(IN,K) QUIT
|
---|
1014 | . SET:$DATA(VALUE(K)) VALUE=VALUE_VALUE(K)
|
---|
1015 | . QUIT
|
---|
1016 | QUIT VALUE
|
---|
1017 | ;===
|
---|
1018 | ;
|
---|
1019 | ;
|
---|
1020 | SEC(X,PREC) ;
|
---|
1021 | New L,LIM,K,SIGN,VALUE
|
---|
1022 | ;
|
---|
1023 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
1024 | ;" Winfried Gerum (8 June 1995)
|
---|
1025 | ;" Comment: The official description does not mention than
|
---|
1026 | ;" the function may also be called with the first
|
---|
1027 | ;" parameter in degrees, minutes and seconds.
|
---|
1028 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
1029 | ;;;
|
---|
1030 | ;
|
---|
1031 | Set PREC=$Get(PREC,11)
|
---|
1032 | Set X=X#(2*$%PI^MATH())
|
---|
1033 | Set (VALUE,L)=1,SIGN=-1
|
---|
1034 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
1035 | For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
1036 | . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
|
---|
1037 | . Quit
|
---|
1038 | If 'VALUE Quit "INFINITE"
|
---|
1039 | Quit 1/VALUE
|
---|
1040 | ;===
|
---|
1041 | ;
|
---|
1042 | ;
|
---|
1043 | SECH(X,PREC) ;;;Quit 1/$%COSH^MATH(X,PREC) ;" Number ~~
|
---|
1044 | ;" Winfried Gerum (8 June 1995)
|
---|
1045 | ;" Alan Frank (October 1995)
|
---|
1046 | Quit 1/$%COSH^MATH(X,$Get(PREC,11))
|
---|
1047 | ;;;
|
---|
1048 | ;===
|
---|
1049 | ;
|
---|
1050 | ;
|
---|
1051 | SIGN(X) Quit $SELECT(X<0:-1,X>0:1,1:0)
|
---|
1052 | ;===
|
---|
1053 | ;
|
---|
1054 | ;
|
---|
1055 | SIN(X,PREC) ;
|
---|
1056 | New L,LIM,K,SIGN,VALUE
|
---|
1057 | ;
|
---|
1058 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
1059 | ;" Winfried Gerum (8 June 1995)
|
---|
1060 | ;" Comment: The official description does not mention than
|
---|
1061 | ;" the function may also be called with the first
|
---|
1062 | ;" parameter in degrees, minutes and seconds.
|
---|
1063 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
1064 | ;;;
|
---|
1065 | ;
|
---|
1066 | Set PREC=$Get(PREC,11)
|
---|
1067 | Set X=X#(2*$%PI^MATH())
|
---|
1068 | Set (VALUE,L)=X,SIGN=-1
|
---|
1069 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
1070 | For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
1071 | . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
|
---|
1072 | . Quit
|
---|
1073 | Quit VALUE
|
---|
1074 | ;===
|
---|
1075 | ;
|
---|
1076 | ;
|
---|
1077 | SIN(X) ;
|
---|
1078 | ;;;" ;" Number ~~
|
---|
1079 | ;" Winfried Gerum (8 June 1995)
|
---|
1080 | ;" Comment: This version of the function is
|
---|
1081 | ;" optimized for speed, not for precision.
|
---|
1082 | ;" The 'precision' parameter is not supported,
|
---|
1083 | ;" and the precision is at best 1 in 10**-9.
|
---|
1084 | ;" Note that this function does not accept its
|
---|
1085 | ;" parameter in degrees, minutes and seconds.
|
---|
1086 | ;;;
|
---|
1087 | ;
|
---|
1088 | New A,N,PI,R,SIGN,XX
|
---|
1089 | ;
|
---|
1090 | ;" This approximation only works for 0 <= x <= pi/2
|
---|
1091 | ;" so reduce angle to correct quadrant.
|
---|
1092 | ;
|
---|
1093 | Set PI=$%PI^MATH(),X=X#(PI*2),SIGN=1
|
---|
1094 | Set:X>PI X=2*PI-X,SIGN=-1
|
---|
1095 | ;
|
---|
1096 | ;;;" Set:X*2<PI X=PI-X Set X=-PI/2+2 ;" Number ~~
|
---|
1097 | ;" Winfried Gerum (8 June 1995)
|
---|
1098 | Set:X*2<PI X=PI-X
|
---|
1099 | ;;;
|
---|
1100 | ;
|
---|
1101 | ;
|
---|
1102 | Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
|
---|
1103 | Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
|
---|
1104 | Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
|
---|
1105 | Quit R*SIGN
|
---|
1106 | ;===
|
---|
1107 | ;
|
---|
1108 | ;
|
---|
1109 | SINH(X,PREC) ;
|
---|
1110 | ;
|
---|
1111 | ;;;" New F,I,P,R,T,XX ;" Number ~~
|
---|
1112 | ;" Winfried Gerum (8 June 1995)
|
---|
1113 | ;" Eli Reidler (28 June 1996)
|
---|
1114 | New E,F,I,P,R,T,XX
|
---|
1115 | ;;;
|
---|
1116 | ;
|
---|
1117 | Set PREC=$Get(PREC,11)+1
|
---|
1118 | Set @("E=1E-"_PREC)
|
---|
1119 | Set XX=X*X,F=1,I=2,(P,R,T)=X
|
---|
1120 | 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
|
---|
1121 | Quit R
|
---|
1122 | ;===
|
---|
1123 | ;
|
---|
1124 | ;
|
---|
1125 | SQRT(X,PREC) ;
|
---|
1126 | If X<0 Set $Ecode=",M28,"
|
---|
1127 | If X=0 Quit 0
|
---|
1128 | ;
|
---|
1129 | ;;;" ;" Number ~~
|
---|
1130 | ;" Alan Frank (October 1995)
|
---|
1131 | Set PREC=$Get(PREC,11)
|
---|
1132 | ;;;
|
---|
1133 | ;
|
---|
1134 | ;
|
---|
1135 | ;;;" If X<1 Quit 1/$%SQRT^MATH(1/X) ;" Number ~~
|
---|
1136 | ;" Winfried Gerum (8 June 1995)
|
---|
1137 | If X<1 Quit 1/$%SQRT^MATH(1/X,PREC)
|
---|
1138 | ;;;
|
---|
1139 | ;
|
---|
1140 | New P,R,E
|
---|
1141 | Set PREC=$Get(PREC,11)+1
|
---|
1142 | ;
|
---|
1143 | ;;;" Set @(E="1E-"_PREC) ;" Number ~~
|
---|
1144 | ;" Winfried Gerum (8 June 1995)
|
---|
1145 | ;" Eli Reidler (28 June 1996)
|
---|
1146 | Set @("E=1E-"_PREC)
|
---|
1147 | ;;;
|
---|
1148 | ;
|
---|
1149 | Set R=X
|
---|
1150 | For Set P=R,R=X/R+R/2,P=P-R/R If -E<P,P<E Quit
|
---|
1151 | Quit R
|
---|
1152 | ;===
|
---|
1153 | ;
|
---|
1154 | ;
|
---|
1155 | TAN(X,PREC) ;
|
---|
1156 | New L,LIM,K,S,SIGN,VALUE
|
---|
1157 | ;
|
---|
1158 | ;;;" Set:X[":" X=$%DMSDEC^MATH(X,12) ;" Number ~~
|
---|
1159 | ;" Winfried Gerum (8 June 1995)
|
---|
1160 | ;" Comment: The official description does not mention than
|
---|
1161 | ;" the function may also be called with the first
|
---|
1162 | ;" parameter in degrees, minutes and seconds.
|
---|
1163 | Set:X[":" X=$%DMSDEC^MATH(X)
|
---|
1164 | ;;;
|
---|
1165 | ;
|
---|
1166 | Set PREC=$Get(PREC,11)
|
---|
1167 | Set X=X#(2*$%PI^MATH())
|
---|
1168 | Set (VALUE,L)=X,SIGN=-1
|
---|
1169 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
1170 | For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
1171 | . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
|
---|
1172 | . Quit
|
---|
1173 | Set S=VALUE
|
---|
1174 | Set X=X#(2*$%PI^MATH())
|
---|
1175 | Set (VALUE,L)=1,SIGN=-1
|
---|
1176 | Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
|
---|
1177 | For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
|
---|
1178 | . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
|
---|
1179 | . Quit
|
---|
1180 | If 'VALUE Quit "INFINITE"
|
---|
1181 | Quit S/VALUE
|
---|
1182 | ;===
|
---|
1183 | ;
|
---|
1184 | ;
|
---|
1185 | TANH(X,PREC) ;
|
---|
1186 | ;
|
---|
1187 | ;;;" ;" Number ~~
|
---|
1188 | ;" Alan Frank (October 1995)
|
---|
1189 | Set PREC=$Get(PREC,11)
|
---|
1190 | ;;;
|
---|
1191 | ;
|
---|
1192 | Quit $%SINH^MATH(X,PREC)/$%COSH^MATH(X,PREC)
|
---|
1193 | ;===
|
---|
1194 | ;
|
---|
1195 | ;
|
---|