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