100'  NAME--DIGAMA
110'
120'  DESCRIPTION--FINDS THE DIGAMMA OF Z FOR Z>0.
130'
140'  SOURCE--NATIONAL BUREAU OF STANDARDS PUBLICATION AMS55.
150'  THE PROGRAM IS BASED ON THE USE OF EQUATIONS 6.3.5,6.3.6, AND 
160'  6.3.18, WHICHEVER IS THE CLOSEST. DERIVATIVES HIGHER THAN 
170'  THE FIRST ARE FOUND FROM EQUATION 6.4.10.
180'  PROGRAM BY DEAN MYRON TRIBUS, THAYER SCHOOL OF ENGINEERING,
190'  DARTMOUTH COLLEGE, HANOVER, N.H. 03755
200'
210'  INSTRUCTIONS--ENTER DATA IN LINES 360 AND FOLLOWING THE VALUES FOR
220'  WHICH YOU WISH TO KNOW THE DIGAMMA FUNCTION.
230'
240'
250'  * *  *  *  *  *  MAIN PRGRAM  *  *  *  *  *  *  * *  *  *  *
260'
270 FOR I=0 TO 5 
280 READ J(I)
290 NEXT I
300 DATA 1, 1.2, 1.4, 1.6, 1.8, 2.0
310 FOR I=0 TO 5
320 FOR M=0 TO 8
330 READ C(M,I)
340 NEXT M
350 NEXT I
360 DATA -0.5772156,1.64493,-1.20205,1.08232,-1.03693,1.01734
370 DATA -1.00835,1.00408,-1.00201
380 DATA -0.2890399,1.26738,-.739016,.540832,-.425521,.344914
390 DATA -.283439, .234494,-.194666
400 DATA -.0613845,1.02536,-.494562,.3033376,-.20172,.138896
410 DATA -9.72778E-2,6.87339E-2,-4.87972E-2
420 DATA .1260475,.858432,-.351784,.185124,-.106281,.063462
430 DATA -3.86574E-2,2.38038E-2,-1.47473E-2
440 DATA .2849914,.736974,-.26193,.120409,-6.06949E-2,.031936
450 DATA -1.71866E-2,9.36679E-3,-5.14285E-3
460 DATA .4227843,.644934,-.202055,8.23229E-2,-3.69277E-2,.017343
470 DATA -8.34925E-3,4.07735E-3,-2.00838E-3
480 READ Z
490 IF Z<0 THEN 550
500 IF Z<1 THEN 570
510 IF Z<2 THEN 610
520 IF Z<50 THEN 640
530 LET D=LOG(Z)-(1/(2*Z))-(1/(12*Z^2))+(1/(120*Z^4))-(1/(252*Z^6))
540 GO TO 690
550 PRINT "Z NEGATIVE"
560 GO TO 700
570 LET J=1+Z
580 GO SUB 710
590 LET D=D-1/Z
600 GO TO 690
610 LET J=Z
620 GO SUB 710
630 GO TO 690
640 LET J=1+Z-INT(Z)
650 GO SUB 710
660 FOR I=1 TO INT(Z-1)
670 LET D=D+(1/(J+I-1))
680 NEXT I
690 PRINT "Z="Z, "D="D
700 GO TO 480
710 REMARK: THIS SUBROUTINE FINDS DIGAMMA OF J FOR 1<J<2
720 FOR K=0 TO 5
730 IF ABS(J-1-(K/5))>.1 THEN 760
740 LET I=K
750 GO TO 770
760 NEXT K
770 LET D=0
780 LET X=J-1-(I/5)
790 LET T=1
800 FOR M=0 TO 8
810 LET D=D+C(M,I)*T
820 LET T=T*X
830 IF ABS(T)<1E-8 THEN 850
840 NEXT M
850 PRINT
860 RETURN
870 DATA 49,50,51
880END