Trailing-Edge
-
PDP-10 Archives
-
decuslib20-01
-
decus/20-0020/allrts.fct
There are 2 other files named allrts.fct in the archive. Click here to see a list.
100' NAME--ALLROOTS
110'
120' DESCRIPTION--COMPUTES INTEGRAL AND FRACTIONAL POWERS OF COMPLEX
130' NUMBERS.
140'
150' SOURCE--PETER W. MEIGS '69
160'
170' INSTRUCTIONS--THE PROGRAM WILL ASK FOR FOUR NUMBERS: A,B,M,N WHERE
180' A IS THE REAL PART OF THE NUMBER TO BE EXPONENTIATED AND B IS
190' THE IMAGINARY PART.
200' M/N FORMS THE POWER TO WHICH A+IB WILL BE RAISED WHERE M
210' IS THE NUMERATOR AND N IS THE DENOMINATOR OF THE NUMBER. FOR
220' EXAMPLE: 1,2,3,4 AS INPUT WILL RESULT IN RAISING 1+2*SQR(-1)
230' TO THE .75 POWER (3/4).
240'
250' NOTE--THIS PROGRAM WAS NOT DESIGNED FOR USE WITH NUMBERS OF
260' EITHER A LARGE OR SMALL MAGNITUDE. THOSE REASONABLY NEAR
270' UNITY WILL GIVE BEST RESULTS.
280' IT SHOULD ALSO BE NOTED THAT M/N SHOULD BE A REDUCED FRACTION.
290'
300'
310' * * * * * * * *MAIN PROGRAM* * * * * * * * * * * *
320'
330 PRINT"TYPE A,B,M,N FOR (A+I*B)^(M/N)";
340 INPUT A,B,N1,N
350 IF ABS(A)+ABS(B)<>0 THEN 380
360 PRINT "0,0 NOT ALLOWED--RETYPE";
370 GO TO 330
380 REM CHECK THE EXP PART.
390 IF N1*N=0 THEN 410
400 IF N1=INT(N1) THEN 430
410 PRINT "EXP MUST BE TWO NON-ZERO INTEGERS "
420 GO TO 330
430 IF N<>INT(N) THEN 410
440 IF N<10 THEN 460
450 PRINT"NOTE THAT "N;"LINES OF OUTPUT WILL BE PRINTED."
460 LET P=3.14159265358979323
470 LET M=SQR(A*A+B*B)
480 LET A1=A/M
490 LET B1=B/M
500 IF ABS(A1)<1E-6 THEN 530
510 LET T1=ATN(SQR(1-A1*A1)/A1)
520 GO TO 540
530 LET T1=P/2
540 IF ABS(B1*B1-1)<1E-6 THEN 570
550 LET T2=ATN(B1/SQR(1-B1*B1))
560 GO TO 580
570 LET T2=P/2
580 LET A2=M*COS(T1)
590 LET B2=M*SIN(T2)
600 IF ABS(A-A2)<1E-6 THEN 630
610 LET T1=P-T1
620 GO TO 580
630 IF ABS(B-B2)<1E-6 THEN 660
640 LET T2=-T2
650 GO TO 580
660 IF ABS(T1-T2)<1E-6 THEN 800
670 IF ABS(-T1-T2)>1E-6 THEN 700
680 LET T1=-T1
690 GO TO 800
700 IF ABS(P-T2-T1)>1E-6 THEN 730
710 LET T2=P-T2
720 GO TO 800
730 IF ABS(T1-T2-P)>1E-6 THEN 770
740 LET T2=-P-T2
750 LET T1=-T1
760 GO TO 800
770 PRINT"THERE HAS BEEN AN ERROR. PLEASE NOTIFY KIEWIT."
780 PRINT T1/P;T2/P
790 STOP
800 LET M1=M^(N1/N)
810 PRINT
820 PRINT
830 PRINT"ALL VALUES OF (",A,"+I*",B,")^(",N1,"/";
840 PRINT N,") ARE :"
850 PRINT
860 FOR K=0 TO N-1
870 LET A3=M1*COS(N1/N*(T1+2*P*K))
880 LET A4=M1*SIN(N1/N*(T1+2*P*K))
890 LET A3=INT(A3*1E5+.5)/1E5
900 LET A4=INT(A4*1E5+.5)/1E5
910 PRINT "RE="A3,"IM="A4
920 NEXT K
930 END