100 DIM A(57, 56), C(57), H(56), Z(56)
110 CLEAR
115 INPUT "Order of fit:", P$: P% = VAL(P$): IF (P% > 56 OR P% < 1) THEN BEEP: GOTO 115
120 GOSUB 245
125 FOR J% = 0 TO P%
130 GOSUB 235
135 C(0) = 1
140 FOR N% = 0 TO P%
145 N1% = N% + 1
150 IF N% = J% THEN N% = N% + 1
155 IF N% > P% THEN 190
160 B0 = -N%
165 FOR I% = 1 TO N1%
170 A(I%, J%) = B0 * C(I%) + C(I% - 1)
175 NEXT I%
180 A(0, J%) = B0 * C(0)
185 FOR I% = 0 TO N1%: C(I%) = A(I%, J%): NEXT I%
190 NEXT N%
195 GOSUB 220
200 FOR I% = 0 TO P%: A(I%, J%) = A(I%, J%) / C(0): NEXT I%
205 NEXT J%
210 GOSUB 275
215 BEEP: PRINT "finished": STOP
220 FOR I% = P% TO 1 STEP -1
225 C(I% - 1) = C(I%) * J% + C(I% - 1): NEXT I%
230 RETURN
235 FOR I% = 0 TO P% + 1: C(I%) = 0: NEXT I%
240 RETURN
245 PRINT : PRINT "INPUT DATA": PRINT "==========": PRINT
250 FOR I% = 0 TO P%
255 DISP: PRINT "Y value for abscissa value:", I%: INPUT Z(I%)
260 PRINT "Z("; I%; ")= "; Z(I%)
265 NEXT I%
270 RETURN
275 PRINT : PRINT "POLYNOMIAL FIT": PRINT "==============": PRINT
280 FOR I% = 0 TO P%: H(I%) = 0: NEXT I%
285 FOR J% = 0 TO P%
290 FOR I% = 0 TO P%
295 H(I%) = Z(J%) * A(I%, J%) + H(I%)
300 NEXT I%
305 NEXT J%
310 FOR I% = 0 TO P%: PRINT "H("; I%; ")= "; H(I%): NEXT I%
315 RETURN
320 END