Annotation of freem/mlib/%ulmath.m, revision 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