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