Annotation of freem/mlib/%ulmath.m, revision 1.1.1.1
1.1 snw 1: %ulmath ; MATH library - version 0.5.0.1
2: ;
3: ; Unless otherwise noted, the code below
4: ; was approved in document X11/95-11
5: ;
6: ; CORRECTION INFORMATION
7: ; Due to the overall routine size this has not been provided
8: ; It is available from Ed de Moel who provided these functions
9: ;
10: ABS(X) Quit $Translate(+X,"-")
11: ;===
12: ;
13: ;
14: %ARCCOS(X) ;
15: ;
16: New A,N,R,SIGN,XX
17: If X<-1 Set $Ecode=",M28,"
18: If X>1 Set $Ecode=",M28,"
19: Set SIGN=1 Set:X<0 X=-X,SIGN=-1
20: Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
21: Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
22: Set A(6)=0.0066700901,A(7)=-0.0012624911
23: Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
24: ;
25: Set R=$$SQRT(1-X,11)*R
26: ;;;
27: ;
28: Quit R*SIGN
29: ;===
30: ;
31: ;
32: ARCCOS(X,PREC) ;
33: I '$D(PREC) Q $$%ARCCOS(X)
34: ;
35: New L,LIM,K,SIG,SIGS,VALUE
36: ;;;
37: ;
38: If X<-1 Set $Ecode=",M28,"
39: If X>1 Set $Ecode=",M28,"
40: Set PREC=$Get(PREC,11)
41: ;
42: If $Translate(X,"-")=1 Quit 0
43: ;;;
44: ;
45: Set SIG=$Select(X<0:-1,1:1),VALUE=1-(X*X)
46: ;
47: Set X=$$SQRT(VALUE,PREC)
48: ;;;
49: ;
50: If $Translate(X,"-")=1 Do Quit VALUE
51: . ;;;
52: . ;
53: . Set VALUE=$$PI()/2*X
54: . Quit
55: ;
56: If X>0.9 Do Quit VALUE
57: . ;;;
58: . ;
59: . Set SIGS=$Select(X<0:-1,1:1)
60: . Set VALUE=1/(1/X/X-1)
61: . ;
62: . Set X=$$SQRT(VALUE,PREC)
63: . ;;;
64: . ;
65: . ;
66: . Set VALUE=$$ARCTAN(X,PREC)*SIGS
67: . ;;;
68: ;
69: . Quit
70: Set (VALUE,L)=X
71: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
72: For K=3:2 Do Quit:($Translate(L,"-")<LIM)
73: . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
74: . Quit
75: Quit $Select(SIG<0:$$PI()-VALUE,1:VALUE)
76: ;===
77: ;
78: ;
79: ARCCOSH(X,PREC) ;
80: If X<1 Set $Ecode=",M28,"
81: New SQ
82: ;
83: Set PREC=$Get(PREC,11)
84: ;;;
85: ;
86: Set SQ=$$SQRT(X*X-1,PREC)
87: Quit $$LOG(X+SQ,PREC)
88: ;===
89: ;
90: ;
91: ARCCOT(X,PREC) ;
92: Set PREC=$Get(PREC,11)
93: Set X=1/X
94: Quit $$ARCTAN(X,PREC)
95: ;===
96: ;
97: ;
98: ARCCOTH(X,PREC) ;
99: New L1,L2
100: ;
101: Set PREC=$Get(PREC,11)
102: ;;;
103: ;
104: Set L1=$$LOG(X+1,PREC)
105: Set L2=$$LOG(X-1,PREC)
106: Quit L1-L2/2
107: ;===
108: ;
109: ;
110: ARCCSC(X,PREC) ;
111: Set PREC=$Get(PREC,11)
112: Set X=1/X
113: Quit $$ARCSIN(X,PREC)
114: ;===
115: ;
116: ;
117: ARCSEC(X,PREC) ;
118: Set PREC=$Get(PREC,11)
119: Set X=1/X
120: Quit $$ARCCOS(X,PREC)
121: ;===
122: ;
123: ;
124: %ARCSIN(X) ;
125: ;;; ; Number ~~
126: ; Winfried Gerum (8 June 1995)
127: ;
128: New A,N,R,SIGN,XX
129: If X<-1 Set $Ecode=",M28,"
130: If X>1 Set $Ecode=",M28,"
131: Set SIGN=1 Set:X<0 X=-X,SIGN=-1
132: Set A(0)=1.5707963050,A(1)=-0.2145988016,A(2)=0.0889789874
133: Set A(3)=-0.0501743046,A(4)=0.0308918810,A(5)=-0.0170881256
134: Set A(6)=0.0066700901,A(7)=-0.0012624911
135: Set R=A(0),XX=1 For N=1:1:7 Set XX=XX*X,R=A(N)*XX+R
136: ;
137: Set R=$$SQRT(1-X,11)*R
138: ;;;
139: ;
140: Set R=$$PI()/2-R
141: Quit R*SIGN
142: ;===
143: ;
144: ;
145: ARCSIN(X,PREC) ;
146: I '$D(PREC) Q $$%ARCSIN(X)
147: New L,LIM,K,SIGS,VALUE
148: Set PREC=$Get(PREC,11)
149: ;
150: If $Translate(X,"-")=1 Do Quit VALUE
151: . ;;;
152: . ;
153: . Set VALUE=$$PI()/2*X
154: . Quit
155: ;
156: If X>0.99999 Do Quit VALUE
157: . ;;;
158: . ;
159: . Set SIGS=$Select(X<0:-1,1:1)
160: . Set VALUE=1/(1/X/X-1)
161: . ;
162: . ;;;
163: . ;
164: . Set VALUE=$$ARCTAN(X,PREC)*SIGS
165: . ;;;
166: . ;
167: . Quit
168: Set (VALUE,L)=X
169: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
170: For K=3:2 Do Quit:($Translate(L,"-")<LIM)
171: . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
172: . Quit
173: Quit VALUE
174: ;===
175: ;
176: ;
177: ARCSINH(X,PREC) ;
178: If X<1 Set $Ecode=",M28,"
179: New SQ
180: ;
181: Set PREC=$Get(PREC,11)
182: ;;;
183: ;
184: Set SQ=$$SQRT(X*X+1,PREC)
185: Quit $$LOG(X+SQ,PREC)
186: ;===
187: ;
188: ;
189: ARCTAN(X,PREC) ;
190: New FOLD,HI,L,LIM,LO,K,SIGN,SIGS,SIGT,VALUE
191: Set PREC=$Get(PREC,11)
192: Set LO=0.0000000001,HI=9999999999
193: Set SIGT=$Select(X<0:-1,1:1),X=$Translate(X,"-")
194: Set X=$Select(X<LO:LO,X>HI:HI,1:X)
195: ;
196: Set FOLD=$Select(X'<1:0,1:1)
197: ;;;
198: ;
199: Set X=$Select(FOLD:1/X,1:X)
200: Set L=X,VALUE=$$PI()/2-(1/X),SIGN=1
201: ;
202: If X<1.3 Do Quit VALUE
203: . ;;;
204: . ;
205: . Set X=$Select(FOLD:1/X,1:X),VALUE=1/((1/X/X)+1)
206: . ;
207: . Set X=$$SQRT(VALUE,PREC)
208: . ;;;
209: . ;
210: . If $Translate(X,"-")=1 Do Quit
211: . . Set VALUE=$$PI()/2*X
212: . . Quit
213: . If X>0.9 Do Quit
214: . . Set SIGS=$Select(X<0:-1,1:1)
215: . . Set VALUE=1/(1/X/X-1)
216: . . Set X=$$SQRT(VALUE)
217: . . Set VALUE=$$ARCTAN(X,10)
218: . . Set VALUE=VALUE*SIGS
219: . . Quit
220: . Set (VALUE,L)=X
221: . Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
222: . For K=3:2 Do Quit:($Translate(L,"-")<LIM)
223: . . Set L=L*X*X*(K-2)/(K-1)*(K-2)/K,VALUE=VALUE+L
224: . . Quit
225: . Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
226: . Quit
227: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
228: For K=3:2 Do Quit:($Translate(1/L,"-")<LIM)
229: . ;
230: . Set L=L*X*X,VALUE=VALUE+(1/(K*L)*SIGN)
231: . ;;;
232: . ;
233: . Set SIGN=SIGN*-1
234: . Quit
235: Set VALUE=$Select(FOLD:$$PI()/2-VALUE,1:VALUE)
236: Set VALUE=$Select(SIGT<1:-VALUE,1:VALUE)
237: Quit VALUE
238: ;===
239: ;
240: ;
241: ARCTANH(X,PREC) ;
242: If X<-1 Set $Ecode=",M28,"
243: If X>1 Set $Ecode=",M28,"
244: ;
245: Set PREC=$Get(PREC,11)
246: ;;;
247: ;
248: Quit $$LOG(1+X/(1-X),PREC)/2
249: ;===
250: ;
251: ;
252: CABS(Z) ;
253: New ZRE,ZIM
254: Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
255: Quit $$SQRT(ZRE*ZRE+(ZIM*ZIM))
256: ;===
257: ;
258: ;
259: CADD(X,Y) ;
260: New XRE,XIM,YRE,YIM
261: Set XRE=+X,XIM=+$Piece(X,"%",2)
262: Set YRE=+Y,YIM=+$Piece(Y,"%",2)
263: Quit XRE+YRE_"%"_(XIM+YIM)
264: ;===
265: ;
266: ;
267: CCOS(Z,PREC) ;
268: New E1,E2,IA
269: ;
270: Set PREC=$Get(PREC,11)
271: ;;;
272: ;
273: Set IA=$$CMUL(Z,"0%1")
274: Set E1=$$CEXP(IA,PREC)
275: Set IA=-IA_"%"_(-$Piece(IA,"%",2))
276: Set E2=$$CEXP(IA,PREC)
277: Set IA=$$CADD(E1,E2)
278: Quit $$CMUL(IA,"0.5%0")
279: ;===
280: ;
281: ;
282: CDIV(X,Y) ;
283: New D,IM,RE,XIM,XRE,YIM,YRE
284: Set XRE=+X,XIM=+$Piece(X,"%",2)
285: Set YRE=+Y,YIM=+$Piece(Y,"%",2)
286: Set D=YRE*YRE+(YIM*YIM)
287: Set RE=XRE*YRE+(XIM*YIM)/D
288: Set IM=XIM*YRE-(XRE*YIM)/D
289: Quit RE_"%"_IM
290: ;===
291: ;
292: ;
293: CEXP(Z,PREC) ;
294: New R,ZIM,ZRE
295: ;
296: Set PREC=$Get(PREC,11)
297: ;;;
298: ;
299: Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
300: Set R=$$EXP(ZRE,PREC)
301: Quit R*$$COS(ZIM,PREC)_"%"_(R*$$SIN(ZIM,PREC))
302: ;===
303: ;
304: ;
305: CLOG(Z,PREC) ;
306: New ABS,ARG,ZIM,ZRE
307: ;
308: Set PREC=$Get(PREC,11)
309: ;;;
310: ;
311: Set ABS=$$CABS(Z)
312: Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
313: ;
314: Set ARG=$$ARCTAN(ZIM/ZRE,PREC)
315: ;;;
316: ;
317: Quit $$LOG(ABS,PREC)_"%"_ARG
318: ;===
319: ;
320: ;
321: CMUL(X,Y) ;
322: New XIM,XRE,YIM,YRE
323: Set XRE=+X,XIM=+$Piece(X,"%",2)
324: Set YRE=+Y,YIM=+$Piece(Y,"%",2)
325: Quit XRE*YRE-(XIM*YIM)_"%"_(XRE*YIM+(XIM*YRE))
326: ;===
327: ;
328: ;
329: COMPLEX(X) Quit +X_"%0"
330: ;===
331: ;
332: ;
333: CONJUG(Z) ;
334: New ZIM,ZRE
335: Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
336: Quit ZRE_"%"_(-ZIM)
337: ;===
338: ;
339: ;
340: COS(X,PREC) ;
341: I '$D(PREC) Q $$%COS(X)
342: New L,LIM,K,SIGN,VALUE
343: ;
344: Set:X[":" X=$$DMSDEC(X)
345: ;;;
346: ;
347: Set PREC=$Get(PREC,11)
348: Set X=X#(2*$$PI())
349: Set (VALUE,L)=1,SIGN=-1
350: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
351: For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
352: . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
353: . Quit
354: Quit VALUE
355: ;===
356: ;
357: ;
358: %COS(X) ;
359: ;;; ; Number ~~
360: ; Winfried Gerum (8 June 1995)
361: ;
362: New A,N,PI,R,SIGN,XX
363: ;
364: ; This approximation only works for 0 <= x <= pi/2
365: ; so reduce angle to correct quadrant.
366: ;
367: Set PI=$$PI(),X=X#(PI*2),SIGN=1
368: Set:X>PI X=2*PI-X
369: Set:X*2>PI X=PI-X,SIGN=-1
370: ;
371: Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
372: Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
373: Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
374: Quit R*SIGN
375: ;===
376: ;
377: ;
378: COSH(X,PREC) ;
379: ;
380: New E,F,I,P,R,T,XX
381: ;;;
382: ;
383: Set PREC=$Get(PREC,11)+1
384: Set @("E=1E-"_PREC)
385: Set XX=X*X,F=1,(P,R,T)=1,I=1
386: 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
387: Quit R
388: ;===
389: ;
390: ;
391: COT(X,PREC) ;
392: New C,L,LIM,K,SIGN,VALUE
393: Set:X[":" X=$$DMSDEC(X)
394: ;;;
395: ;
396: Set PREC=$Get(PREC,11)
397: Set (VALUE,L)=1,SIGN=-1
398: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
399: For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
400: . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
401: . Quit
402: Set C=VALUE
403: Set X=X#(2*$$PI())
404: Set (VALUE,L)=X,SIGN=-1
405: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
406: For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
407: . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
408: . Quit
409: If 'VALUE Quit "INFINITE"
410: Quit VALUE=C/VALUE
411: ;===
412: ;
413: ;
414: COTH(X,PREC) ;
415: New SINH
416: If 'X Quit "INFINITE"
417: ;
418: Set PREC=$Get(PREC,11)
419: ;;;
420: ;
421: Set SINH=$$SINH(X,PREC)
422: If 'SINH Quit "INFINITE"
423: Quit $$COSH(X,PREC)/SINH
424: ;===
425: ;
426: ;
427: CPOWER(Z,N,PREC) ;
428: New AR,NIM,NRE,PHI,PI,R,RHO,TH,ZIM,ZRE
429: ;
430: Set PREC=$Get(PREC,11)
431: ;;;
432: ;
433: Set ZRE=+Z,ZIM=+$Piece(Z,"%",2)
434: Set NRE=+N,NIM=+$Piece(N,"%",2)
435: If 'ZRE,'ZIM,'NRE,'NIM Set $Ecode=",M28,"
436: ;
437: If 'ZRE,'ZIM Quit "0%0"
438: ;;;
439: ;
440: Set PI=$$PI()
441: ;
442: Set R=$$SQRT(ZRE*ZRE+(ZIM*ZIM),PREC)
443: ;;;
444: ;
445: ;
446: If ZRE Set TH=$$ARCTAN(ZIM/ZRE,PREC)
447: ;;;
448: ;
449: Else Set TH=$SELECT(ZIM>0:PI/2,1:-PI/2)
450: ;;;
451: ;
452: Set RHO=$$LOG(R,PREC)
453: Set AR=$$EXP(RHO*NRE-(TH*NIM),PREC)
454: Set PHI=RHO*NIM+(NRE*TH)
455: Quit AR*$$COS(PHI,PREC)_"%"_(AR*$$SIN(PHI,PREC))
456: ;===
457: ;
458: ;
459: CSC(X,PREC) ;
460: New L,LIM,K,SIGN,VALUE
461: ;
462: ; Winfried Gerum (8 June 1995)
463: Set:X[":" X=$$DMSDEC(X)
464: ;;;
465: ;
466: Set PREC=$Get(PREC,11)
467: ;;;
468: ;
469: Set X=X#(2*$$PI())
470: Set (VALUE,L)=X,SIGN=-1
471: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
472: For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
473: . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
474: . Quit
475: If 'VALUE Quit "INFINITE"
476: Quit 1/VALUE
477: ;===
478: ;
479: ;
480: ;
481: CSCH(X,PREC)
482: Quit 1/$$SINH(X,$Get(PREC,11))
483: ;;;
484: ;
485: ;===
486: ;
487: ;
488: CSIN(Z,PREC) ;
489: New IA,E1,E2
490: ;
491: Set PREC=$Get(PREC,11)
492: ;;;
493: ;
494: Set IA=$$CMUL(Z,"0%1")
495: Set E1=$$CEXP(IA,PREC)
496: Set IA=-IA_"%"_(-$Piece(IA,"%",2))
497: Set E2=$$CEXP(IA,PREC)
498: Set IA=$$CSUB(E1,E2)
499: Set IA=$$CMUL(IA,"0.5%0")
500: Quit $$CMUL("0%-1",IA)
501: ;===
502: ;
503: ;
504: CSUB(X,Y) ;
505: New XIM,XRE,YIM,YRE
506: Set XRE=+X,XIM=+$Piece(X,"%",2)
507: Set YRE=+Y,YIM=+$Piece(Y,"%",2)
508: Quit XRE-YRE_"%"_(XIM-YIM)
509: ;===
510: ;
511: ;
512: DECDMS(X,PREC) ;
513: Set PREC=$Get(PREC,5)
514: Set X=X#360*3600
515: Set X=+$Justify(X,0,$Select((PREC-$Length(X\1))'<0:PREC-$Length(X\1),1:0))
516: Quit X\3600_":"_(X\60#60)_":"_(X#60)
517: ;===
518: ;
519: ;
520: DEGRAD(X) Quit X*3.14159265358979/180
521: ;===
522: ;
523: ;
524: DMSDEC(X) ;
525: Quit $Piece(X,":")+($Piece(X,":",2)/60)+($Piece(X,":",3)/3600)
526: ;===
527: ;
528: ;
529: E() Quit 2.71828182845905
530: ;===
531: ;
532: ;
533: EXP(X,PREC) ;
534: New L,LIM,K,VALUE
535: Set PREC=$Get(PREC,11)
536: Set L=X,VALUE=X+1
537: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
538: For K=2:1 Set L=L*X/K,VALUE=VALUE+L Quit:($Translate(L,"-")<LIM)
539: Quit VALUE
540: ;===
541: ;
542: ;
543: LOG(X,PREC) ;
544: New L,LIM,M,N,K,VALUE
545: If X'>0 Set $Ecode=",M28,"
546: Set PREC=$Get(PREC,11)
547: Set M=1
548: ;
549: For N=0:1 Quit:(X/M)<10 Set M=M*10
550: ;;;
551: ;
552: If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
553: Set X=X/M
554: Set X=(X-1)/(X+1),(VALUE,L)=X
555: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
556: For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
557: Set VALUE=VALUE*2+(N*2.30258509298749)
558: Quit VALUE
559: ;===
560: ;
561: ;
562: LOG10(X,PREC) ;
563: New L,LIM,M,N,K,VALUE
564: If X'>0 Set $Ecode=",M28,"
565: Set PREC=$Get(PREC,11)
566: Set M=1
567: ;
568: For N=0:1 Quit:(X/M)<10 Set M=M*10
569: ;;;
570: ;
571: If X<1 For N=0:-1 Quit:(X/M)>0.1 Set M=M*0.1
572: Set X=X/M
573: Set X=(X-1)/(X+1),(VALUE,L)=X
574: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
575: For K=3:2 Set L=L*X*X,M=L/K,VALUE=M+VALUE Set:M<0 M=-M Quit:M<LIM
576: Set VALUE=VALUE*2+(N*2.30258509298749)
577: Quit VALUE/2.30258509298749
578: ;===
579: ;
580: ;
581: MTXADD(A,B,R,ROWS,COLS) ;
582: ; Add A[ROWS,COLS] to B[ROWS,COLS],
583: ; result goes to R[ROWS,COLS]
584: IF $DATA(A)<10 QUIT 0
585: IF $DATA(B)<10 QUIT 0
586: IF $GET(ROWS)<1 QUIT 0
587: IF $GET(COLS)<1 QUIT 0
588: ;
589: NEW ROW,COL,ANY
590: FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
591: . KVALUE R(ROW,COL) SET ANY=0
592: . SET:$DATA(A(ROW,COL))#2 ANY=1
593: . SET:$DATA(B(ROW,COL))#2 ANY=1
594: . SET:ANY R(ROW,COL)=$GET(A(ROW,COL))+$GET(B(ROW,COL))
595: . QUIT
596: QUIT 1
597: ;===
598: ;
599: ;
600: MTXCOF(A,I,K,N) ;
601: ; Compute cofactor for element [i,k]
602: ; in matrix A[N,N]
603: NEW T,R,C,RR,CC
604: SET CC=0 FOR C=1:1:N DO:C'=K
605: . SET CC=CC+1,RR=0
606: . FOR R=1:1:N SET:R'=I RR=RR+1,T(RR,CC)=$GET(A(R,C))
607: . QUIT
608: QUIT $$MTXDET(.T,N-1)
609: ;===
610: ;
611: ;
612: MTXCOPY(A,R,ROWS,COLS) ;
613: ; Copy A[ROWS,COLS] to R[ROWS,COLS]
614: IF $DATA(A)<10 QUIT 0
615: IF $GET(ROWS)<1 QUIT 0
616: IF $GET(COLS)<1 QUIT 0
617: ;
618: NEW ROW,COL
619: FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
620: . KVALUE R(ROW,COL)
621: . SET:$DATA(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)
622: . QUIT
623: QUIT 1
624: ;===
625: ;
626: ;
627: MTXDET(A,N) ;
628: ; Compute determinant of matrix A[N,N]
629: IF $DATA(A)<10 QUIT ""
630: IF N<1 QUIT ""
631: ;
632: ; First the simple cases
633: ;
634: IF N=1 QUIT $GET(A(1,1))
635: IF N=2 QUIT $GET(A(1,1))*$GET(A(2,2))-($GET(A(1,2))*$GET(A(2,1)))
636: ;
637: NEW DET,I,SIGN
638: ;
639: ; Det A = sum (k=1:n) element (i,k) * cofactor [i,k]
640: ;
641: SET DET=0,SIGN=1
642: FOR I=1:1:N DO
643: . SET DET=$GET(A(1,I))*$$MTXCOF(.A,1,I,N)*SIGN+DET
644: . SET SIGN=-SIGN
645: . QUIT
646: QUIT DET
647: ;===
648: ;
649: ;
650: MTXEQU(A,B,R,N,M) ;
651: ; Solve matrix equation A [M,M] * R [M,N] = B [M,N]
652: IF $GET(M)<1 QUIT ""
653: IF $GET(N)<1 QUIT ""
654: IF '$$MTXDET(.A) QUIT 0
655: ;
656: NEW I,I1,J,J1,J2,K,L,T,T1,T2,TEMP,X
657: ;
658: SET X=$$MTXCOPY(.A,.T,N,N)
659: SET X=$$MTXCOPY(.B,.R,N,M)
660: ;
661: ; Reduction of matrix A
662: ; Steps of reduction are counted by index K
663: ;
664: FOR K=1:1:N-1 DO
665: . ;
666: . ; Search for largest coefficient of T
667: . ; (denoted by TEMP)
668: . ; in first column of reduced system
669: . ;
670: . SET TEMP=0,J2=K
671: . FOR J1=K:1:N DO
672: . . QUIT:$TRANSLATE($GET(T(J1,K)),"-")>$TRANSLATE(TEMP,"-")
673: . . SET TEMP=T(J1,K),J2=J1
674: . . QUIT
675: . ;
676: . ; Exchange row number K with row number J2,
677: . ; if necessary
678: . ;
679: . DO:J2'=K
680: . . ;
681: . . FOR J=K:1:N DO
682: . . . SET T1=$GET(T(K,J)),T2=$GET(T(J2,J))
683: . . . KILL T(K,J),T(J2,J)
684: . . . IF T1'="" SET T(J2,J)=T1
685: . . . IF T2'="" SET T(K,J)=T2
686: . . . QUIT
687: . . FOR J=1:1:M DO
688: . . . SET T1=$GET(R(K,J)),T2=$GET(R(J2,J))
689: . . . KILL R(K,J),R(J2,J)
690: . . . IF T1'="" SET R(J2,J)=T1
691: . . . IF T2'="" SET R(K,J)=T2
692: . . . QUIT
693: . . QUIT
694: . ;
695: . ; Actual reduction
696: . ;
697: . FOR I=K+1:1:N DO
698: . . FOR J=K+1:1:N DO
699: . . . QUIT:'$GET(T(K,K))
700: . . . SET T(I,J)=-$GET(T(K,J))*$GET(T(I,K))/T(K,K)+$GET(T(I,J))
701: . . . QUIT
702: . . FOR J=1:1:M DO
703: . . . QUIT:'$GET(T(K,K))
704: . . . SET R(I,J)=-$GET(R(K,J))*$GET(T(I,K))/T(K,K)+$GET(R(I,J))
705: . . . QUIT
706: . . QUIT
707: . QUIT
708: ;
709: ; Backsubstitution
710: ;
711: FOR J=1:1:M DO
712: . IF $GET(T(N,N)) SET R(N,J)=$GET(R(N,J))/T(N,N)
713: . IF N-1>0 FOR I1=1:1:N-1 DO
714: . . SET I=N-I1
715: . . FOR L=I+1:1:N DO
716: . . . SET R(I,J)=-$GET(T(I,L))*$GET(R(L,J))+$GET(R(I,J))
717: . . . QUIT
718: . . IF $GET(T(I,I)) SET R(I,J)=$GET(R(I,J))/$GET(T(I,I))
719: . . QUIT
720: . QUIT
721: QUIT $$MTXDET(.R)
722: ;===
723: ;
724: ;
725: MTXINV(A,R,N) ;
726: ; Invert A[N,N], result goes to R[N,N]
727: IF $DATA(A)<10 QUIT 0
728: IF $GET(N)<1 QUIT 0
729: ;
730: NEW T,X
731: SET X=$$MTXUNIT(.T,N)
732: QUIT $$MTXEQU(.A,.T,.R,N,N)
733: ;===
734: ;
735: ;
736: MTXMUL(A,B,R,M,L,N) ;
737: ; Multiply A[M,L] by B[L,N], result goes to R[M,N]
738: IF $DATA(A)<10 QUIT 0
739: IF $DATA(B)<10 QUIT 0
740: IF $GET(L)<1 QUIT 0
741: IF $GET(M)<1 QUIT 0
742: IF $GET(N)<1 QUIT 0
743: ;
744: NEW I,J,K,SUM,ANY
745: FOR I=1:1:M FOR J=1:1:N DO
746: . SET (SUM,ANY)=0
747: . KVALUE R(I,J)
748: . FOR K=1:1:L DO
749: . . SET:$DATA(A(I,K))#2 ANY=1
750: . . SET:$DATA(B(K,J))#2 ANY=1
751: . . SET SUM=$GET(A(I,K))*$GET(B(K,J))+SUM
752: . . QUIT
753: . SET:ANY R(I,J)=SUM
754: . QUIT
755: QUIT 1
756: ;===
757: ;
758: ;
759: MTXSCA(A,R,ROWS,COLS,S) ;
760: ; Multiply A[ROWS,COLS] with the scalar S,
761: ; result goes to R[ROWS,COLS]
762: IF $DATA(A)<10 QUIT 0
763: IF $GET(ROWS)<1 QUIT 0
764: IF $GET(COLS)<1 QUIT 0
765: IF '($DATA(S)#2) QUIT 0
766: ;
767: NEW ROW,COL
768: FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
769: . KVALUE R(ROW,COL)
770: . SET:$DATA(A(ROW,COL))#2 R(ROW,COL)=A(ROW,COL)*S
771: . QUIT
772: QUIT 1
773: ;===
774: ;
775: ;
776: MTXSUB(A,B,R,ROWS,COLS) ;
777: ; Subtract B[ROWS,COLS] from A[ROWS,COLS],
778: ; result goes to R[ROWS,COLS]
779: IF $DATA(A)<10 QUIT 0
780: IF $DATA(B)<10 QUIT 0
781: IF $GET(ROWS)<1 QUIT 0
782: IF $GET(COLS)<1 QUIT 0
783: ;
784: NEW ROW,COL,ANY
785: FOR ROW=1:1:ROWS FOR COL=1:1:COLS DO
786: . KVALUE R(ROW,COL) SET ANY=0
787: . SET:$DATA(A(ROW,COL))#2 ANY=1
788: . SET:$DATA(B(ROW,COL))#2 ANY=1
789: . ;
790: . SET:ANY R(ROW,COL)=$GET(A(ROW,COL))-$GET(B(ROW,COL))
791: . ;;;
792: . ;
793: . QUIT
794: QUIT 1
795: ;===
796: ;
797: ;
798: MTXTRP(A,R,M,N) ;
799: ; Transpose A[M,N], result goes to R[N,M]
800: IF $DATA(A)<10 QUIT 0
801: IF $GET(M)<1 QUIT 0
802: IF $GET(N)<1 QUIT 0
803: ;
804: NEW I,J,K,D1,V1,D2,V2
805: FOR I=1:1:M+N-1 FOR J=1:1:I+1\2 DO
806: . SET K=I-J+1
807: . IF K=J DO QUIT
808: . . SET V1=$GET(A(J,J)),D1=$DATA(A(J,J))#2
809: . . IF J'>N,J'>M KVALUE R(J,J) SET:D1 R(J,J)=V1
810: . . QUIT
811: . ;
812: . SET V1=$GET(A(K,J)),D1=$DATA(A(K,J))#2
813: . SET V2=$GET(A(J,K)),D2=$DATA(A(J,K))#2
814: . IF K'>M,J'>N KVALUE R(K,J) SET:D2 R(K,J)=V2
815: . IF J'>M,K'>N KVALUE R(J,K) SET:D1 R(J,K)=V1
816: . QUIT
817: QUIT 1
818: ;===
819: ;
820: ;
821: MTXUNIT(R,N,SPARSE) ;
822: ; Create a unit matrix R[N,N]
823: IF $GET(N)<1 QUIT 0
824: ;
825: NEW ROW,COL
826: FOR ROW=1:1:N FOR COL=1:1:N DO
827: . KVALUE R(ROW,COL)
828: . IF $GET(SPARSE) QUIT:ROW'=COL
829: . SET R(ROW,COL)=$SELECT(ROW=COL:1,1:0)
830: . QUIT
831: QUIT 1
832: ;===
833: ;
834: ;
835: PI() Quit 3.14159265358979
836: ;===
837: ;
838: ;
839: RADDEG(X) Quit X*180/3.14159265358979
840: ;===
841: ;
842: ;
843: SEC(X,PREC) ;
844: New L,LIM,K,SIGN,VALUE
845: ;
846: Set:X[":" X=$$DMSDEC(X)
847: ;;;
848: ;
849: Set PREC=$Get(PREC,11)
850: Set X=X#(2*$$PI())
851: Set (VALUE,L)=1,SIGN=-1
852: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
853: For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
854: . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
855: . Quit
856: If 'VALUE Quit "INFINITE"
857: Quit 1/VALUE
858: ;===
859: ;
860: ;
861: SECH(X,PREC)
862: Quit 1/$$COSH(X,$Get(PREC,11))
863: ;;;
864: ;===
865: ;
866: ;
867: SIGN(X) Quit $SELECT(X<0:-1,X>0:1,1:0)
868: ;===
869: ;
870: ;
871: SIN(X,PREC) ;
872: I '$D(PREC) Q $$%SIN(X)
873: New L,LIM,K,SIGN,VALUE
874: ;
875: Set:X[":" X=$$DMSDEC(X)
876: ;;;
877: ;
878: Set PREC=$Get(PREC,11)
879: Set X=X#(2*$$PI())
880: Set (VALUE,L)=X,SIGN=-1
881: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
882: For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
883: . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
884: . Quit
885: Quit VALUE
886: ;===
887: ;
888: ;
889: %SIN(X) ;
890: ;
891: New A,N,PI,R,SIGN,XX
892: ;
893: ; This approximation only works for 0 <= x <= pi/2
894: ; so reduce angle to correct quadrant.
895: ;
896: Set PI=$$PI(),X=X#(PI*2),SIGN=1
897: Set:X>PI X=2*PI-X,SIGN=-1
898: ;
899: Set:X*2<PI X=PI-X
900: ;;;
901: ;
902: ;
903: Set XX=X*X,A(1)=-0.4999999963,A(2)=0.0416666418
904: Set A(3)=-0.0013888397,A(4)=0.0000247609,A(5)=-0.0000002605
905: Set (X,R)=1 For N=1:1:5 Set X=X*XX,R=A(N)*X+R
906: Quit R*SIGN
907: ;===
908: ;
909: ;
910: SINH(X,PREC) ;
911: ;
912: New E,F,I,P,R,T,XX
913: ;;;
914: ;
915: Set PREC=$Get(PREC,11)+1
916: Set @("E=1E-"_PREC)
917: Set XX=X*X,F=1,I=2,(P,R,T)=X
918: 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
919: Quit R
920: ;===
921: ;
922: ;
923: SQRT(X,PREC) ;
924: If X<0 Set $Ecode=",M28,"
925: If X=0 Quit 0
926: ;
927: Set PREC=$Get(PREC,11)
928: ;;;
929: ;
930: ;
931: If X<1 Quit 1/$$SQRT(1/X,PREC)
932: ;;;
933: ;
934: New P,R,E
935: Set PREC=$Get(PREC,11)+1
936: ;
937: Set @("E=1E-"_PREC)
938: ;;;
939: ;
940: Set R=X
941: For Set P=R,R=X/R+R/2,P=P-R/R If -E<P,P<E Quit
942: Quit R
943: ;===
944: ;
945: ;
946: TAN(X,PREC) ;
947: New L,LIM,K,S,SIGN,VALUE
948: ;
949: Set:X[":" X=$$DMSDEC(X)
950: ;;;
951: ;
952: Set PREC=$Get(PREC,11)
953: Set X=X#(2*$$PI())
954: Set (VALUE,L)=X,SIGN=-1
955: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
956: For K=3:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
957: . Set L=L/(K-1)*X/K*X,VALUE=VALUE+(SIGN*L)
958: . Quit
959: Set S=VALUE
960: Set X=X#(2*$$PI())
961: Set (VALUE,L)=1,SIGN=-1
962: Set LIM=$Select((PREC+3)'>11:PREC+3,1:11),@("LIM=1E-"_LIM)
963: For K=2:2 Do Quit:($Translate(L,"-")<LIM) Set SIGN=SIGN*-1
964: . Set L=L*X*X/(K-1*K),VALUE=VALUE+(SIGN*L)
965: . Quit
966: If 'VALUE Quit "INFINITE"
967: Quit S/VALUE
968: ;===
969: ;
970: ;
971: TANH(X,PREC) ;
972: ;
973: Set PREC=$Get(PREC,11)
974: ;;;
975: ;
976: Quit $$SINH(X,PREC)/$$COSH(X,PREC)
977: ;===
978:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to %ulmath.m CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Coherent Logic Development] / freem / mlib