100'  NAME--DGMINV
110'  
120'  DESCRIPTION--FINDS THE VALUE OF Z TO MAKE D(Z)=D0 WHERE
130'  D(Z) IS THE DIGAMMA FUNCTION OF Z, (Z>0). THIS IS THE INVERSE
140'  OF PROGRAM DIGAMA.
150'
160'  SOURCE--DEAN MYRON TRIBUS, THAYER SCHOOL OF ENGINEERING,
170'  DARTMOUTH COLLEGE, HANOVER, N.H. 03755
180'
190'  INSTRUCTIONS--ENTER DATA IN LINE 350 AND FOLLOWING THE VALUES OF
200'  D0 FOR WHICH YOU WISH TO KNOW Z.
210'
220'
230'  * *  *  *  *   MAIN PROGRAM   *  *  *  *  *  *  *  *  *  *
240'
250 DIM Z(16),D(16)
260 FOR I=0 TO 5 
270 READ J(I)
280 NEXT I
290 DATA 1, 1.2, 1.4, 1.6, 1.8, 2.0
300 FOR I=0 TO 5
310 FOR M=0 TO 8
320 READ C(M,I)
330 NEXT M
340 NEXT I
350 DATA -0.5772156,1.64493,-1.20205,1.08232,-1.03693,1.01734
360 DATA -1.00835,1.00408,-1.00201
370 DATA -0.2890399,1.26738,-.739016,.540832,-.425521,.344914
380 DATA -.283439, .234494,-.194666
390 DATA -.0613845,1.02536,-.494562,.3033376,-.20172,.138896
400 DATA -9.72778E-2,6.87339E-2,-4.87972E-2
410 DATA .1260475,.858432,-.351784,.185124,-.106281,.063462
420 DATA -3.86574E-2,2.38038E-2,-1.47473E-2
430 DATA .2849914,.736974,-.26193,.120409,-6.06949E-2,.031936
440 DATA -1.71866E-2,9.36679E-3,-5.14285E-3
450 DATA .4227843,.644934,-.202055,8.23229E-2,-3.69277E-2,.017343
460 DATA -8.34925E-3,4.07735E-3,-2.00838E-3
470 FOR I=0 TO 16
480 READ Z(I),D(I)
490 NEXT I
500 DATA .001,-1000.58,.01,-100.561,.02,-50.5448,.05,-20.4978
510 DATA .1,-10.4238,.2,-5.28904,.5,-1.96351,.75,-1.08586,1,-.577216
520 DATA 2,.422784,4,1.25612,6,1.70612,10,2.25175,20,2.97052,30,3.38444
530 DATA 50,3.90199,100,4.60016
540 READ D0
550 FOR I=1 TO 16
560 IF D0>D(I) THEN 610
570 LET Z1=Z(I-1)+(Z(I)-Z(I-1))*(D0-D(I-1))/(D(I)-D(I-1))
580 LET Z2=Z(I-1)
590 LET D2=D(I-1)
600 GO TO 620
610 NEXT I
620 LET Z=Z1
630 GO SUB 690
640 IF ABS(D-D0)>1E-4 THEN 670
650 PRINT "D="D, "Z="Z
660 GO TO 540
670 LET Z1=Z2+(D0-D2)*(Z-Z2)/(D-D2)
680 GO TO 620
690 REM START OF SUBROUTINE TO FIND DIGAMMA FUNCTION(USES SUB 770)
700 IF Z<0 THEN 760
710 IF Z<1 THEN 780
720 IF Z<2 THEN 820
730 IF Z<50 THEN 850
740 LET D=LOG(Z)-(1/(2*Z))-(1/(12*Z^2))+(1/(120*Z^4))-(1/(252*Z^6))
750 GO TO 900
760 PRINT "Z NEGATIVE"
770 GO TO 1080
780 LET J=1+Z
790 GO SUB 910
800 LET D=D-1/Z
810 GO TO 900
820 LET J=Z
830 GO SUB 910
840 GO TO 900
850 LET J=1+Z-INT(Z)
860 GO SUB 910
870 FOR I=1 TO INT(Z-1)
880 LET D=D+(1/(J+I-1))
890 NEXT I
900 RETURN
910 REMARK: THIS SUBROUTINE FINDS DIGAMMA OF J FOR 1<J<2
920 FOR K=0 TO 5
930 IF ABS(J-1-(K/5))>.1 THEN 960
940 LET I=K
950 GO TO 970
960 NEXT K
970 LET D=0
980 LET X=J-1-(I/5)
990 LET T=1
1000 FOR M=0 TO 8
1010 LET D=D+C(M,I)*T
1020 LET T=T*X
1030 IF ABS(T)<1E-8 THEN 1050
1040 NEXT M
1050 REM 
1060 RETURN
1070 DATA 1.5
1080 END