source: cprs/branches/tmg-cprs/m_files/TMGMATH.m@ 1227

Last change on this file since 1227 was 796, checked in by Kevin Toppenberg, 14 years ago

Initial upload

File size: 28.7 KB
Line 
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
17ABS(X) Quit $Translate(+X,"-")
18 ;===
19 ;
20 ;
21ARCCOS(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 ;
48ARCCOS(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 ;
111ARCCOSH(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 ;
126ARCCOT(X,PREC) ;
127 Set PREC=$Get(PREC,11)
128 Set X=1/X
129 Quit $%ARCTAN^MATH(X,PREC)
130 ;===
131 ;
132 ;
133ARCCOTH(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 ;
148ARCCSC(X,PREC) ;
149 Set PREC=$Get(PREC,11)
150 Set X=1/X
151 Quit $%ARCSIN^MATH(X,PREC)
152 ;===
153 ;
154 ;
155ARCSEC(X,PREC) ;
156 Set PREC=$Get(PREC,11)
157 Set X=1/X
158 Quit $%ARCCOS^MATH(X,PREC)
159 ;===
160 ;
161 ;
162ARCSIN(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 ;
190ARCSIN(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 ;
232ARCSINH(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 ;
247ARCTAN(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 ;
309ARCTANH(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 ;
323CABS(Z) ;
324 New ZRE,ZIM
325 Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
326 Quit $%SQRT^MATH(ZRE*ZRE+(ZIM*ZIM))
327 ;===
328 ;
329 ;
330CADD(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 ;
338CCOS(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 ;
355CDIV(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 ;
366CEXP(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 ;
380CLOG(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 ;
400CMUL(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 ;
408COMPLEX(X) Quit +X_"%0"
409 ;===
410 ;
411 ;
412CONJUG(Z) ;
413 New ZIM,ZRE
414 Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
415 Quit ZRE_"%"_(-ZIM)
416 ;===
417 ;
418 ;
419COS(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 ;
441COS(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 ;
468COSH(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 ;
483COT(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 ;
512COTH(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 ;
528CPOWER(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 ;
571CSC(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 ;
599CSCH(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 ;
608CSIN(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 ;
626CSUB(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 ;
634DECDMS(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 ;
642DEGRAD(X) Quit X*3.14159265358979/180
643 ;===
644 ;
645 ;
646DMSDEC(X) ;
647 Quit $Piece(X,":")+($Piece(X,":",2)/60)+($Piece(X,":",3)/3600)
648 ;===
649 ;
650 ;
651E() Quit 2.71828182845905
652 ;===
653 ;
654 ;
655EXP(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 ;
665LOG(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 ;
686LOG10(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 ;
707MTXADD(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 ;
726MTXCOF(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 ;
738MTXCOPY(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 ;
753MTXDET(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 ;
776MTXEQU(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 ;
856MTXINV(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 ;
867MTXMUL(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 ;
890MTXSCA(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 ;
907MTXSUB(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 ;
931MTXTRP(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 ;
954MTXUNIT(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 ;
968PI() Quit 3.14159265358979
969 ;===
970 ;
971 ;
972PRODUCE(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 ;
994RADDEG(X) Quit X*180/3.14159265358979
995 ;===
996 ;
997 ;
998REPLACE(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 ;
1020SEC(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 ;
1043SECH(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 ;
1051SIGN(X) Quit $SELECT(X<0:-1,X>0:1,1:0)
1052 ;===
1053 ;
1054 ;
1055SIN(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 ;
1077SIN(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 ;
1109SINH(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 ;
1125SQRT(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 ;
1155TAN(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 ;
1185TANH(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 ;
Note: See TracBrowser for help on using the repository browser.