100'  NAME--FIB
110'
120'  DESCRIPTION--LOCATES AND EVALUATES THE MAXIMUM OF A UNIMODAL
130'  FUNCTION OF ONE VARIABLE WITHIN A SPECIFIED INTERVAL
140'  OF THAT VARIABLE.
150'
160'  SOURCE--HOWARD ROBERTS THAYER SCHOOL OF ENGINEERING, 
170'  REVISED 4/17/68 BY MALCOLM LEWIS. THE PROGRAM USES THE FIBONACCI
180'  SEARCH TECHNIQUE OUTLINES IN CHAPTER 4 OF "OPTIMIZATION", BY
190'  A. 0. CONVERSE.
200'
210'  INSTRUCTIONS--THE USER MUST FURNISH:
220'    1) INTERVAL TO BE EXAMINED (A,B) 
230'    2) MINIMUM ALLOWABLE SEPARATION OF EVALUATIONS (E) 
240'    3) MAXIMUM ALLOWABLE FINAL INTERVAL SIZE (H) 
250'      (ENTER FOUR POINTS IN LINE 820) 
260'    4) FUNCTION TO BE MAXIMIZED
270'      (ENTER IN LINE 1370-1390) 
280'
290'  THE PROGRAM WILL RETURN THE FOLLOWING INFORMATION:
300'    1) F1,X1,F2,X2 FOR EACH INTERVAL (OPTIONAL) 
310'    2) NUMBER OF EXPERIMENTS REQUIRED TO ESTABLISH THE POSITION
320'      OF F(MAX)  TO SPECIFIED TOLERANCES.
330'    3) INTERVAL IN WHICH FUNCTION MAXIMUM LIES (X,Y) 
340'    4) VALUE OF FUNCTION AT MIDDLE OF FINAL INTERVAL
350'  NOTE:FOR GENERAL INFORMATION ON NOMENCLATURE AND COMPUTATIONAL
360'  PROCEDURE CONTINUE LISTING UNTIL LINE 760.
370'
380'
390'  *  *  *  *  *  *  *  *  MAIN PROGRAM  *  *  *  *  *  *  *  *  *  *  *
400'
410 ' GENERAL INFORMATION:
420 '     NOMENCLATURE:
430 '          A=LOWER LIMIT OF INITIAL INTERVAL
440 '          B=UPPER LIMIT OF INITIAL INTERVAL
450 '          C=WIDTH OF INITIAL INTERVAL
460 '          E=MINIMUM ALLOWED SEPARATION OF EVALUATIONS
470 '          F=CURRENT VALUE OF FUNCTION
480 '          F1=EVALUATION #1 FOR CURRENT INTERVAL
490 '          F2=EVALUATION #2 FOR CURRENT INTERVAL
500 '          G1,G2,G3=FIBONACCI NUMBERS IN CURRENT USE
510 '          H=MAXIMUM ALLOWABLE FINAL INTERVAL SIZE
520 '          H1=COMPUTED WIDTH OF FINAL INTERVAL
530 '          I=COUNTER IN "L2 COMPUTATION" SUBROUTINE
540 '          J=COMPUTED TOTAL NUMBER OF EXPERIMENTS REQUIRED
550 '          K$=INPUT SWITCH ON INTERMEDIATE PRINTOUT OPTION
560 '          L1,L2=INTERVAL FRACTION POSITION INDICES
570 '          N=COUNTER IN MAIN SEARCH ROUTINE
580 '          X=CURRENT VALUE OF INDEPENDENT VARIABLE
590 '          X1,X2=LOCATIONS OF F1, F2
600 '          Z=INTERMEDIATE VARIABLE IN MAIN SEARCH ROUTINE
610 '     
620 '     COMPUTATION PROCEDURE:
630 '          F1,F2 ARE PLACED AT X1>X2.  F1,F2 ARE COMPARED AND
640 '            THE INTERVAL (A,B) IS DIMINISHED  BY REDEFINING
650 '            A OR B ACCORDING TO THE FOLLOWING SCHEME:
660 '
670 '   IF F1<F2 THEN             *  IF F1>F2 THEN
680 '       REDEFINE: B=X1        *      REDEFINE: A=X2
690 '                 X1=X2       *                X2=X1
700 '                 F1=F2       *                F2=F1
710 '                 X2=B-L2*C   *                X1=A+L2*C
720 '
730 '          THE ENTIRE CYCLE IS REPEATED J TIMES, AT WHICH POINT
740 '            THE INTERVAL WIDTH IS LESS THAN H.
750'
760'
770'  NOTE THIS IS THE END OF THE GENERAL INFORMATION.
780 ' ****** INPUT ******
790 READ A,B   'START AND END OF INTERVAL
800 READ E,H   ' MINIMUM SEPARATION, MAXIMUM FINAL INTERVAL
810 '
820 DATA 0,10,.1,.5   'A,B,E,H
830 '
840 PRINT "DO YOU WANT INTERMEDIATE PRINTOUT, YES OR NO";
850 INPUT K$
860 '
870 '****** INITIALIZATION ******
880 LET C=B-A   ' SIZE OF INITIAL INTERVAL
890 GOSUB 1420   ' GO COMPUTE L2
900 LET L1=1
910 IF K$ = "NO" THEN 930   ' IF "NO", NO INTERMEDIATE PRINTOUT
920 PRINT "    X1","    F1","    X2","    F2"
930 '
940 LET X=A+L2*C   ' SET X AT DIST(L2*C) TO RIGHT OF POINT A.
950 LET X1=X   ' START X1 AND X2 AT THE SAME POSITION
960 LET X2=X
970 GOSUB 1370   ' EVALUATE F AT X=X1
980 LET F1=F   ' SET F1=F(X1)
990'
1000 ' ****** MAIN SEARCH ROUTINE *******
1010 FOR N=2 TO J   ' PERFORM (J-2) EXPERIMENTS
1020 IF A<X2 THEN 1090    ' FOR A>X2, DECREMENT X. (REMEMBER: X2>X1)
1030 LET X=A+L2*C   ' FOR A>X2, INCREMENT X:  X(I)=A+L(I)*(B-A)
1040 LET X2=X1   ' MOVE X2 TO RIGHT TO X1
1050 LET X1=X   ' MOVE X1 TO NEW HIGHER X-VALUE
1060 GOSUB 1370   ' GO EVALUATE F AT X=X1
1070 LET F1=F   ' SET F1=F(X1)
1080 GO TO 1140
1090 LET X=B-L2*C   ' FOR A<X2, DECREMENT X:  X(I)=B-L(I)*(B-A)
1100 LET X1=X2   ' MOVE X1 TO LEFT TO X2
1110 LET X2=X   ' MOVE X2 TO NEW LOWER X-VALUE
1120 GOSUB 1370   ' GO EVALUATE F AT X=X2
1130 LET F2=F   ' SET F2=F(X2)
1140 IF K$ = "NO" THEN 1160   ' IF "NO", NO INTERMEDIATE PRINTOUT
1150 PRINT X1,F1,X2,F2
1160 '
1170 IF F1<F2 THEN 1210   ' FOR F1>F2, MOVE SEARCH REGION TO RIGHT
1180 LET A=X2   ' RAISE THE LOWER BOUND ON THE SEARCH REGION
1190 LET F2=F1   ' SAVE NEW MAX VALUE OF F(X)
1200 GO TO 1240
1210 LET F1=F2   ' FOR F1<F2, MOVE REGION TO LEFT; SAVE NEW F(MAX)
1220 LET B=X1   ' LOWER THE UPPER BOUND ON THE SEARCH REGION
1230' RECURSIVELY COMPUTE L2
1240 LET Z=L2   ' TEMPORARILY SAVE  OLD VALUE OF L2
1250 LET L2=L1-L2          ' SET L(N+1)=L(N-1)-L(N)
1260 LET L1=Z   ' LET L1= OLD L2-VALUE
1270 NEXT N
1280'
1290 ' ******* FINAL OUTPUT *******
1300 LET X = (A+B)/2  ' COMPUTE MIDPOINT OF FINAL INTERVAL
1310 GO SUB 1370  ' EVALUATE THE FINAL VALUE OF F(X)
1320 PRINT " OPTIMUM VALUE LIES IN THE INTERVAL ";A;" TO ";B
1330 PRINT " VALUE OF FUNCTION AT X = ";X;" IS = ";F
1340 GO TO 1610   '   END OF PROGRAM
1350 ' 
1360 '
1370 ' ********** DEFINE YOUR FUNCTION HERE **********
1380 LET F = -3*X^2 + 21.6*X + 1 ' SAMPLE FUNCTION; REPLACE WITH YOURS
1390 '
1400 RETURN
1410 '
1420 ' ******* SUBR TO COMPUTE L2 *******
1430 '
1440 LET G1 = 1 ' FIRST FIBONACCI NUMBER
1450 LET G2 = 1 ' SECOND   "        "
1460 LET I = 2 ' THE ITERATION STARTS AT I = 2
1470 LET G3 = G2 + G1 ' FIBONACCI RECURSION FORMULA
1480 LET H1 = (C + G1*E)/G3 ' WIDTH OF FINAL INTERVAL FOR I TRIALS
1490 IF H1 <= H THEN 1550 ' IF H1 IS SMALL ENOUGH, GO COMPUTE L2
1500 LET I = I+1 ' H1 TOO LARGE; TRY N = I+1
1510 LET G1 = G2 ' NOW G(I-2) = G(I-1)
1520 LET G2 = G3  ' NOW G(I-1) = G(I)
1530 IF I >= 50 THEN 1550 ' IF N >= 50 COMPUTE L2 REGARDLESS OF  H1
1540 GO TO 1470 ' INCREMENT FIBONACCI SERIES
1550 LET L2 = G2/G3 + ((-1)^I)*E/(G3*C)  ' COMPUTE L2
1560 LET J = I ' SET TOTAL NUMBER OF EXPERIMENTS EQUAL TO I
1570 PRINT
1580 PRINT " # OF EXPERIMENTS = ";J
1590 PRINT " ORIGINAL INTERVAL = ";A;" TO ";B
1600 RETURN
1610 END