20 + - * / ^ = SQR LOG EXP ABS SIN COS TAN INT RETURN
25 '
30 ' << MATHEMATICAL FUNCTIONS --- MATHF.BAS --- by Dr Russell Langley >>
35 '
40 CLEAR:DEF SEG:KEY OFF:ON ERROR GOTO 1610:DEFINT I-N:DIM X(401),Y(401),Z(401),T$(6),TC$(6),U$(9),UC$(9):HD$=" M A T H E M A T I C A L F U N C T I O N S ":VER$="(RL,8)"
65 LOCATE 3,8:PRINT "Press the NUMBER of your choice, and then press <ENTER> to run it.":PRINT TAB(22)"If trouble, try <Ctrl-Break> & GOTO 9.":PRINT STRING$(80,196):LOCATE 7,1:K=14
70 PRINT:PRINT TAB(K)"1. Derivative of any Function"CHR$(10)TAB(K)"2. Exponential Curve Fitting, Y = K + A * (B ^ X)"CHR$(10)TAB(K)"3. Factorials, Permutations, & Combinations"CHR$(10)TAB(K)"4. Function Plotter, Y = f(X)"
115 INPUT" (Y/N)";YN$:IF YN$="" THEN YN$="N":RETURN ELSE YN$=CHR$(ASC(YN$) AND 95):IF YN$="Y" OR YN$="N" THEN RETURN ELSE PRINT"WHAT? ";:GOTO 115
120 IF OP<>4 THEN PRINT TAB(46)"Press ANY KEY to continue";:Z$=INKEY$ ELSE PRINT TAB(68)"Press a key";:Z$=INKEY$
125 IF INKEY$ =""THEN 125 ELSE IF OP<>4 THEN LOCATE ,39:PRINT SPACE$(35):RETURN ELSE LOCATE 25,68:PRINT SPACE$(11);:RETURN
130 ' -- Read row I of A(I,J) for J=1 to M, in FF.
135 INPUT X$:IF X$="" THEN PRINT"WHAT";:GOTO 135 ELSE K=0
140 FOR J=1 TO M:Y$="":QP=0
145 K=K+1:IF K=LEN(X$) THEN Y$=Y$+RIGHT$(X$,1):IF J=M THEN 155 ELSE J=M:FLAG=1:GOTO 160
150 Z$=MID$(X$,K,1):IF Z$>" " THEN Y$=Y$+Z$:QP=1:GOTO 145 ELSE IF QP=0 THEN 145
155 A(I,J)=VAL(Y$)
160 NEXT J:IF FLAG=1 THEN FLAG=0:BEEP:PRINT"Not enough values. Please re-do whole line.":GOTO 135 ELSE RETURN
165 ' -- Line Counter
170 L=L+1:IF L>20 THEN GOSUB 120:L=1
175 RETURN
180 ' -- Permutation Sub
185 M=N1:E=0:IF X=0 THEN RETURN
190 FOR J=X TO 1 STEP -1:IF M>=NT THEN M=M/NT:E=E+T
195 M=M*N:N=N-N1:NEXT J:RETURN
200 '
205 ' -- FUNCTION PLOTTER
210 '
215 CLS:Z$="F U N C T I O N P L O T T E R":GOSUB 110
220 PRINT:PRINT:PRINT"After the `?', enter required function Y=f(X) in BASIC, starting with `Y='":PRINT"Any fractions must have 1 or more digits before the decimal point.":PRINT
225 IF T$(1)="+" THEN 255
230 DATA "+","-","*","/","^","=",SQR,LOG,EXP,ABS,SIN,COS,TAN,INT,RETURN
235 FOR I=1 TO 6:READ T$(I):NEXT I:FOR I=1 TO 9:READ U$(I):NEXT I
240 START=PEEK(&H30)+256*PEEK(&H31):AD1=PEEK(START)+256*PEEK(START+1):AD2=PEEK(AD1)+256*PEEK(AD1+1)+4:FOR I=1 TO 6:IF PEEK(AD2)<>32 THEN TC$(I)=TC$(I)+CHR$(PEEK(AD2)) ELSE I=I-1
245 AD2=AD2+1:NEXT I:AD2=AD2+1:FOR I=1 TO 9:IF PEEK(AD2)<>32 THEN UC$(I)=UC$(I)+CHR$(PEEK(AD2)):I=I-1
250 AD2=AD2+1:NEXT I
255 INPUT Z$:IF Z$="" THEN IF Z2$>"" THEN PRINT "Ok, I will use: "Z2$:PRINT:GOTO 285 ELSE BEEP:GOTO 255
260 FOR I=1 TO LEN(Z$):A=ASC(MID$(Z$,I,1)):IF A>96 THEN A=A-32:MID$(Z$,I,1)=CHR$(A)
265 NEXT I:PRINT:Z2$=Z$
270 FU$=Z$:AD=AD1+4:FOR I=1 TO LEN(FU$):FOR J=1 TO 6:IF MID$(FU$,I,1)=T$(J) THEN MID$(FU$,I,1)=TC$(J)
275 NEXT J:NEXT I:FOR J=1 TO 8:U=INSTR(FU$,U$(J)):IF U>0 THEN FU$=LEFT$(FU$,U-1)+UC$(J)+MID$(FU$,U+3)
280 NEXT J:FU$=FU$+":"+UC$(9)+CHR$(0):FOR I=1 TO LEN(FU$):POKE AD,ASC(MID$(FU$,I,1)):AD=AD+1:NEXT I
285 IF OP<>4 THEN RETURN
290 IT=1:INPUT"Smallest X ";SX:INPUT"Biggest X ";BX:IF BX<=SX THEN PRINT"Silly":GOTO 290
295 PRINT:PRINT W$:XI=(BX-SX)/400:X=SX:GOSUB 15:SY=Y:BY=SY:FOR I=0 TO 400:X(I)=SX+I*XI:X=X(I):GOSUB 15:Y(I)=Y:IF Y(I)<SY THEN SY=Y(I) ELSE IF Y(I)>BY THEN BY=Y(I)
300 NEXT I:B=BX:S=SX:H=5:IF ABS(BY-SY)<9.999E-06 THEN PRINT"Y = Constant,"BY:GOTO 435
305 '
310 ' -- Select Scales (H = No. of Axis Labels)
315 '
320 Q=1:KT=0
325 R=B-S:C=S
330 IF R<=1 THEN KT=KT+1:R=R*10:GOTO 330
335 IF R>10 THEN KT=KT-1:R=R/10:GOTO 335
340 IF Q>2 THEN 350 ELSE C=C*10^KT:IF C<0 AND C<>INT(C) THEN C=C-1
345 C=INT(C)/10^KT:R=(B-C)/H:KT=0:Q=Q+2:GOTO 330
350 F=INT(R):IF F<>R THEN F=F+1
355 IF R<0.5 THEN F=F-0.5
360 F=F/10^KT:IF Q=4 THEN 365 ELSE IF(B-S)/(H*F)<=0.8 THEN KT=1:Q=2:GOTO 325
365 S=C:D=F*INT(C/F):IF D<0 AND D<>C THEN D=D-F
370 IF D+H*F>=B THEN S=D
375 IT=IT+1:IF IT=2 THEN X=S:U=F:H=4:B=BY:S=SY:GOTO 320 ELSE Y=S:V=F
630 CLS:Z$="E X P O N E N T I A L C U R V E F I T T I N G":GOSUB 110:MS$="METHOD OF SEMIAVERAGES":PRINT STRING$(37-LEN(MS$)\2,32)MS$:PRINT:PRINT
635 INPUT"No. of equally spaced X-Values (any multiple of 3) ";N:Z=N/3:IF N<3 OR INT(Z)<>Z THEN 635 ELSE NN=N:M=Z
640 INPUT"First X-Value ";X(1):IF X(1)>100 THEN PRINT"Too big. Code X-Values by subtraction so that all X's are < 100.":GOTO 640 ELSE INPUT"Interval between successive X-Values ";XI:FOR I=2 TO N:X(I)=X(I-1)+XI:NEXT I
645 PRINT:IF N=3 THEN PRINT"Enter 3 Y-Values read from your hand-drawn smooth curve:"ELSE PRINT"ENTER";N;"Observed Y-Values corresponding to observed X-Values:"
650 FOR I=1 TO N
655 PRINT"Y(";MID$(STR$(I),2);") ";:INPUT Y$:IF Y$>" " THEN Y(I)=VAL(Y$):NEXT I ELSE BEEP:GOTO 655
660 K=0:FOR I=1 TO 3:S(I)=0:FOR J=1 TO M:K=K+1:S(I)=S(I)+Y(K):NEXT J:NEXT I:BM=(S(2)-S(3))/(S(1)-S(2)):B=BM^(1/M):A=(S(1)-S(2))*(1-B)/(1-BM)^2:K!=(S(1)*S(3)-S(2)^2)/(M*(S(1)+S(3)-2*S(2))):K!=INT(K!*10000+0.5)/10000
675 PRINT:PRINT"FUNCTION: Y ="K!;:IF A>0 THEN PRINT"+"AA;ELSE PRINT AA;
680 PRINT"* ("BB"^ X )":PRINT:PRINT "Estimates based on this equation:"
685 FOR I=1 TO N:Z(I)=K!+A*B^X(I):NEXT I:PRINT SPACE$(13)"X EST Y":FOR I=1 TO N:PRINT USING F1$;I;X(I);Z(I):IF INT(I/15)=I/15 THEN GOSUB 120
690 NEXT I:PRINT
695 PRINT"Do you want more estimates based on this equation";:GOSUB 115:IF YN$="N" THEN 720ELSE PRINT"Ok, use a slash '/' to signal end."
700 FOR I=N+1 TO 100:INPUT"X ";X$:IF X$="/"THEN NN=I-1:I=100 ELSE X(I)=VAL(X$):Z(I)=K!+A*B^X(I):Z(I)=INT(Z(I)*1000+0.5)/1000:LOCATE CSRLIN-1,15:PRINT "Y ="Z(I)
705 NEXT I:PRINT
710 '
715 ' -- Printout
720 PRINT"Do you want printout of these results";:GOSUB 115:IF YN$="N" THEN 770
725 LINE INPUT "PROBLEM ID ? ";ID$
730 BEEP:LOCATE 24,30:LINE INPUT;"Press <ENTER> when PRINTER is ON.";Z$:LOCATE 24,30:PRINT SPACE$(50);
735 LPRINT CHR$(14)"EXPONENTIAL CURVE FITTING"CHR$(20) '= W I D E letters on/off with Epson Dot Matrix Printers.
740 LPRINT SPACE$(12)MS$
745 LPRINT" ":LPRINT"Problem: "ID$:LPRINT" ":LPRINT"DATA X 0BS Y EST Y"
750 FOR I=1 TO N:LPRINT USING F2$;I;X(I);Y(I);Z(I):NEXT I:LPRINT " "
755 LPRINT"FUNCTION: Y ="K!;:IF A>0 THEN LPRINT"+"AA;ELSE LPRINT AA;
760 LPRINT"* ("BB"^ X )":LPRINT " ":IF NN>N THEN LPRINT"ALSO X EST Y"
765 FOR I=N+1 TO NN:LPRINT USING" ######.# ######.###";X(I);Z(I):NEXT I:LPRINT " ":LPRINT " "
770 PRINT"Do you want another Exponential Curve";:GOSUB 115:IF YN$="Y" THEN 630ELSE 55
775 '
780 ' -- POLYSOLVER (after B. P. Douglass)
785 '
790 CLS:CLEAR:DEFDBL A-Z:DEFINT I-N:Z$="P O L Y S O L V E R":GOSUB 110:PRINT:PRINT SPACE$(5)"This seeks N roots of: Y = A + B(1)*X +B(2)*X^2 + ... + B(N)*X^N":PRINT SPACE$(10)"(Complex or unstable roots will cause non-convergence)":PRINT
795 INPUT"Order N of your polynomial (1-10) ";N:NN=N
800 INPUT"A ";A(0):FOR I=1 TO N:PRINT"B("MID$(STR$(I),2)") ";:INPUT A(I):NEXT I:FOR I=0 TO N:D(I)=A(I):NEXT I
805 PRINT:PRINT"Approx roots will first be found using NEWTON & HORNER'S method.":PRINT"Enter trial value for each root in turn:":FOR J=1 TO NN:K=0
815 GOSUB 860:IF ABS(B(0))<9.99E-07 THEN 825ELSE IF C(1)=0 THEN C(1)=9.9999E-05
820 IF K>29 THEN 850ELSE K=K+1:X=X-B(0)/C(1):GOTO 815
825 X(J)=X:FOR I=1 TO N:A(I-1)=B(I):NEXT I:N=N-1:PRINT"Rough Root #"J"= "CSNG(X)" after "K" iterations.":NEXT J
830 PRINT:PRINT"Improved ANSWERS obtained by re-solving with above approx roots:":N=NN:FOR I=0 TO N:A(I)=D(I):NEXT I:FOR J=1 TO N:PRINT "WORKING";:LOCATE ,1:X=X(J)
835 GOSUB 860:IF ABS(B(0))<9.8E-08 THEN 845ELSE IF C(1)=0 THEN C(1)=9.999E-06
960 IF I/2=INT(I/2)THEN YE=YE+Y(I) ELSE YO=YO+Y(I)
965 NEXT I:GOTO 995
970 PRINT"Enter an Odd Number (3-101) of equally spaced X Values at which Y is known,":INPUT "including the Starting & Ending Values ";N:IF N<3 OR N/2=INT(N/2)THEN PRINT"That's not an odd number.":BEEP:GOTO 970
975 PRINT:N=N-1:W=(XN-X0)/N
980 FOR I=0 TO N:PRINT"At X = ";X0+I*W;TAB(20)"Y("MID$(STR$(I+1),2)") ";:INPUT Y(I)
985 IF I=0 OR I=N THEN 990 ELSE IF I/2=INT(I/2)THEN YE=YE+Y(I) ELSE YO=YO+Y(I)
990 NEXT I:PRINT
995 PRINT"INTEGRAL AREA = ";W/3*(Y(0)+2*YE+4*YO+Y(N))" (SIMPSON)":PRINT
1000 PRINT"Another Integral of this function";:GOSUB 115:IF YN$="Y" THEN 935ELSE 55
1005 '
1010 ' -- INTEGRATION by ROMBERG'S ALGORITHM (after B. P. Douglass)
1015 '
1020 CLEAR:ON ERROR GOTO 1610:GOSUB 220
1025 INPUT"Single or Double precision (S/D) ";Z$:IF Z$="D" OR X$="d" THEN DEFDBL R,S,W,X
1030 INPUT"Start integrating at what X value ";X0:INPUT"End integrating at what X value ";XN:IF XN<=X0 THEN PRINT"SILLY":GOTO 1030
1035 INPUT"How many iterations (2-15) ";N:IF N<2 OR N>15 THEN BEEP:GOTO 1035
1040 DIM R(N,N):W=XN-X0:X=X0:GOSUB 15:Y0=Y:X=XN:GOSUB 15:R(1,1)=W*(Y0+Y)/2:L=1
1045 FOR I=2 TO N:W=W/2:L=L+L:S=0:FOR K=1 TO L-1 STEP 2:X=X0+W*K:GOSUB 15:S=S+Y:NEXT K:R(I,1)=R(I-1,1)/2+W*S:PRINT"R("MID$(STR$(I),2)",1) ="R(I,1);
1050 EM=1:FOR J=2 TO I:EM=EM*4:R(I,J)=R(I,J-1)+(R(I,J-1)-R(I-1,J-1))/(EM-1):IF J>= I THEN I$=MID$(STR$(I),2):PRINT" R("I$","I$") ="R(I,J)
1055 NEXT J
1060 NEXT I:PRINT"FINAL INTEGRAL ="R(N,N)" (ROMBERG)":PRINT:LINE INPUT "Press <ENTER> to return to Math Functions Menu. ";Z$:GOTO 9
1065 '
1070 ' -- SIMULTANEOUS EQUATIONS
1075 '
1080 CLS:Z$=" SOLVING SIMULTANEOUS LINEAR EQUATIONS ":GOSUB 110
1090 PRINT" Given values of X's & Y for N equations, this program solves":PRINT" for the unknown B coefficients, if possible.":PRINT" The computations can be done with:"
1095 PRINT" Single precision (6 digit accuracy), or":PRINT" Double precision (16 digit accuracy, but slower).":PRINT
1100 INPUT"Which do you want (S/D) ";SD$
1105 IF SD$="S" OR SD$="s" THEN A0=9.99E-07:A3=0.000999999 ELSE IF SD$="D" OR SD$="d" THEN DEFDBL A,B,C,D,S,P,Y:A0=9.9999E-13:A3=0.001 ELSE PRINT"Please reply 'S' for Single, or 'D' for Double precision.":GOTO 1100
1110 PRINT
1115 INPUT"How many unknown coefficients, N (2-10) ";N:IF N<2 OR N>10 THEN 1115 ELSE M=N+1
1120 DIM A(N,M),B(N),C(N,N),D(N,N),IP(N),KR(N),KC(N),P(N)
1125 F4$="####.##### ":F5$="#####.#### ":F8$="########.# ":R$="Row##: ":D$="Determinant of X Matrix = "
1130 ' --Enter data from Keyboard.
1135 PRINT:PRINT"Enter"N"rows of data, each with X1 to X"MID$(STR$(N),2)" & Y, in Free Format:":I=1
1140 PRINT"Row"I;:GOSUB 135:I=I+1:IF I<=N THEN 1140
1145 CLS:IF QR=0 THEN PRINT"Data was read as ---" ELSE PRINT"Revised data is---"
1150 FOR I=1 TO N:PRINT USING R$;I;:FOR J=1 TO M:PRINT A(I,J);:NEXT J:PRINT:NEXT I:PRINT
1155 PRINT"Do you want to make any CHANGES to those values";:GOSUB 115:IF YN$="N" THEN 1220 ELSE INPUT"Which row";I
1160 IF I<1 OR I>N THEN PRINT"WHAT? ";:GOTO 1155
1165 ' --Edit a row.
1170 PRINT USING R$;I;:FOR J=1 TO M:PRINT A(I,J);:NEXT J:PRINT
1175 PRINT"Enter DATUM # (1-"MID$(STR$(M),2)") & RIGHT VALUE, in Free Format: ";
1180 INPUT "",X$:IF X$="" THEN 1155 ELSE QR=1:K=0
1185 FOR J=1 TO 2:Y$="":QP=0
1190 K=K+1:IF K=LEN(X$) THEN Y$=Y$+RIGHT$(X$,1):IF J=2 THEN 1200 ELSE J=2:FLAG=1:GOTO 1205
1195 Z$=MID$(X$,K,1):IF Z$>" " THEN Y$=Y$+Z$:QP=1:GOTO 1190 ELSE IF QP=0 THEN 1190
1200 Y(J)=VAL(Y$):IF Y(1)<1 OR Y(1)>M THEN PRINT"WHAT? ";:GOTO 1175
1205 NEXT J:IF FLAG=1 THEN FLAG=0:PRINT "No, please enter 2 VALUES on the same line.":GOTO 1175
1210 J=Y(1):A(I,J)=Y(2):GOTO 1145
1215 ' -- Calc now --
1220 COLOR 23,0:PRINT:PRINT"Working.":COLOR 7,0
1225 ' -- Copy A into Square C of order N.
1230 FOR I=1 TO N:FOR J=1 TO N:C(I,J)=A(I,J):NEXT J:NEXT I
1530 PRINT TAB(10)"This program computes Double Precision n!, nPx, or nCx,":PRINT TAB(10)"using a separate exponent which enables products to greatly":PRINT TAB(10)"exceed the normal upper limit of 1E38.":PRINT
1535 INPUT"Factorial, Permutation, Combination, or Menu (F/P/C/M)";X$:IF X$="" THEN 9 ELSE X$=CHR$(ASC(X$) AND 95):IF X$="M" THEN 9 ELSE IF INSTR("FPC",X$)=0 THEN BEEP:GOTO 1535
1540 INPUT"n";N:IF N<2 THEN BEEP:PRINT"n must be 2 or more!":GOTO 1540 ELSE IF X$<>"F" THEN 1550
1620 BEEP:BEEP:IF ERR=2 AND ERL=15 THEN PRINT" ERROR IN YOUR FORMULA: "Z2$:PRINT" Please re-enter it in correct BASIC (starting with `Y=' )":RESUME 255