10   ' The minimax program attempts to find the Chebychev polynomial
   20   ' that approaches a given function accurately enough by a brute
   30   ' force attack.
   40   ' The target polynomial is represented by its zero's and those
   50   ' are adjusted in such a way that the minmax result is obtained.
   60   ' Note: A trailing P often means 'prime' in the sense of derivative
   70   ' Copyright (c) 1993-2000 Albert van der Horst. All rights reserved.
   80   '
   90   ' ###################################################
  100   ' ################## configuration ##################
  110   ' ###################################################
  120   '
  130   ' -------------- User supplied Configuration --------
  140   '
  150   ' Lower and higher bounds of the region to be approximated
  160   Low#=1:High#=2:goto *END1
  170   '
  180   ' Name of function to be approximated
  190   ' Must be a function, because Point destroys all strings
  200   fnTargetName
  210   return("exp")
  220   '
  230   ' Function to be approximated
  240   fnTarget(X)
  250      Result=exp(X)
  260   return(Result)
  270   '
  280   ' Analytic first derivative of fnTarget
  290   fnTargetP(X)
  300      Result=exp(X)
  310   return(Result)
  320   '
  330   ' Analytic second derivative of fnTarget
  340   fnTargetPP(X)
  350      Result=exp(X)
  360   return(Result)
  370   '
  380   *END1:' This label just to jump over something
  390   '
  400   ' Precision needed in digits after the decimal point
  410   '  or absolute precision
  420   '  Fill in the one you need, leave the others at zero
  430   ' Note : before decimal point you have what is left.
  440   '  If that is not enough, bad luck.
  450   ' Eps# itself in maximum precision
  460   Binprecision%=0
  470   DecPrecision%=10
  480   Eps#=0
  490   QualityWanted%=8:' #bits of agreement in the maxima of the polynomial
  500   '
  510   ' Order of approximation to start with
  520   ' (Affects running time at most)
  530   InitialOrder%=2
  540   '
  550   ' -------------- System supplied Configuration --------
  560   '
  561   HexOutput%=Binprecision%:' Remeber for output format
  570   if and{DecPrecision%=0,Binprecision%=0} then
  580      :DecPrecision%=int(1-log(Eps#)/log(10.))
  590   endif
  600   Aux%=int(max(DecPrecision%/4.8,Binprecision%/16)+0.999999)
  610   ' We need at least 2 words spare precision for the polynomial
  620   ' coefficients plus one word for ourselves
  630   ' Increase as needed.
  640   *Precision::point Aux%+4:'Increase if precision insufficient
  650   ' Do not do the following:
  660   ' word Aux%+2:' About 10 digits before decimal point
  670   ' Because non-trivial polynomials will no longer fit in P !!!!
  680   if Eps#=0. then Eps#=min(10.^(-DecPrecision%),2.^(-Binprecision%))
  690   if DecPrecision%=0 then DecPrecision%=Binprecision%*3\10
  700   EpsD=Eps#/2^QualityWanted%:'Deviation eps for iterations
  710   '
  720   MaxOrder%=50:' maximum order of polynomial admitted
  730   Order%=InitialOrder%:' May later be more
  740   ' X(0..Order%) is actually used
  750   dim X(MaxOrder%+1),Y(MaxOrder%+1):' p(x) passes through X(I),Y(I) I<Order
  760   ' Xex(0..Order%+1) is actually used, Xex(I)<X(I)<Xex(I+1)
  770   dim Interval(MaxOrder%+1):' Length of intervals
  780   dim Xex(MaxOrder%+1),Yex(MaxOrder%+1):' Extrema of the error function
  790   P=poly(0.0):' Declare a polynomial, to be tracking X(),Y()
  800   PP=poly(0.0):' Derivative
  810   PPP=poly(0.0):' Second derivative
  820   ' PP and PPP are its first and second derivatives
  830   '
  840   ' Speeds of convergence : .15 works always
  850   ' .25 rather safe, .5 fast but does not always work.
  860   Fudge=0.15
  870   Quality=0:QL=0.
  880   goto *Main
  890   '
  900   ' ###################################################
  910   ' ################## subroutines ####################
  920   ' ###################################################
  930   '
  940   ' -------------- Calculate P from X,Y ------------------
  950   '
  960   ' They say we should never do that, because X,Y is a better
  970   ' representation. However, we are using fixed point arithmetic here,
  980   ' with generous guard decimals and we need the derivatives.
  990   ' Secondly, P is our goal, so it better be reasonable.
 1000   ' This is adapted from 'Numerical recipes in C' section 3.5
 1010   '
 1020   ' `polcoe' makes the global `P' the interpolating polynomial
 1030   ' on the point X(0..Order%),Y(0..Order%)
 1040   *Polcoe
 1050   ' Uses two local polynomials M (master) and its derivative MP
 1060   '
 1070   M=poly(1.0)
 1080   for I%=0 to Order%
 1090      M*=(_X-X(I%))
 1100   next I%
 1110   '
 1120   MP=diff(M)
 1130   P=poly(0)
 1140   for I%=0 to Order%
 1150      P+=Y(I%)*(M\(_X-X(I%)))/val(MP,X(I%))
 1160   next I%
 1170   return
 1180   '
 1190   *TestPolcoe
 1200   print " For the following : expect x^2+.5*x -1."
 1210   Order%=2
 1220   X(0)=1:Y(0)=0.5
 1230   X(1)=2:Y(1)=4
 1240   X(2)=3:Y(2)=9.5
 1250   gosub *Polcoe
 1260   print P
 1270   print 0,X(0),Y(0),val(P,X(0))
 1280   print 1,X(1),Y(1),val(P,X(1))
 1290   print 2,X(2),Y(2),val(P,X(2))
 1300   return
 1310   '
 1320   ' -------------- Initialise X,Y ------------------------
 1330   '
 1340   ' Initialise the arrays X(0..Order%) Y(0..Order%) such that
 1350   ' [Low,High] is divided in Order%+2 intervals plus two half
 1360   ' intervals at the ends, at the zero's of the Order%+1-th
 1365   ' Chebychev polynomial
 1370   *InitXY
 1380    local I%
 1390    for I%=0 to Order%
 1400       X(I%)=High#+Low#+(High#-Low#)*cos(pi((Order%-I%+0.5)/(Order%+1)))
 1410       X(I%)/=2
 1420       Y(I%)=fnTarget(X(I%))
 1430    next I%
 1440   return
 1450   '
 1460   *TestInitXY
 1470    gosub *InitXY
 1480    gosub *Polcoe
 1490    N%=100
 1500    for I%=0 to 100
 1510       X=Low#+I%*(High#-Low#)/100
 1520       print X,val(P,X)-fnTarget(X),val(PP,X)-fnTargetP(X),val(PPP,X)-fnTarget(X)
 1530    next I%
 1540   return
 1550   '
 1560   ' -------------- Initialise X of extrema ---------------
 1570   '
 1580   ' Initialise the array Xex(0..Order%+1) such that they are
 1590   ' plausible starting positions for the extrema in the intervals
 1600   ' defined by *InitXY. X(i-1)<Xex(i)<X(i) inasfar X(i) defined.
 1610   *InitXex
 1620    local I%
 1630    Xex(0)=Low#
 1640    Xex(Order%+1)=High#
 1650    for I%=1 to Order%
 1660       Xex(I%)=(X(I%-1)+X(I%))/2
 1670    next I%
 1680   return
 1690   '
 1700   '
 1710   ' -------------- Refine extrema ------------------------
 1720   '
 1730   ' Refine the array Xex such that they actually point to the
 1740   ' position of the extrema of the intervals.
 1750   ' As a side effect fill in the extrema themselves into
 1760   ' Yex(0..Order%+1). Uses iteration to arrive at the zero's of
 1770   ' the derivative of the target function minus the
 1780   ' approximation polynomial. (Failing that it uses a direct
 1790   ' 'trisection' approach for the maximum of the target
 1800   ' function minus the approximation polynomial directly).
 1810   ' Assumes that P,P,PPP are set in accordance with X,Y.
 1820   ' Returns the index where an iteration did not converge, else -1
 1830   fnRefineExtrema
 1840   local I%,Iter%
 1850   local X,Y:' Point where we are iterating
 1860   local XL,XH,YL,YH:' Bracket containing the maximum
 1870   local X1,X2,Y1,Y2:' Auxiliaries for trisection
 1880   '
 1890   for I%=1 to Order%
 1900     Iter%=0
 1910     XL=X(I%-1):YL=val(P,XL)-fnTarget(XL)
 1920     XH=X(I%):YH=val(P,XH)-fnTarget(XH)
 1930     X=Xex(I%)
 1940     YP=val(PP,X)-fnTargetP(X)
 1950     while abs(YP)>EpsD
 1960        ' print Iter%,X,YP,val(P,X)-fnTarget(X)
 1970        X-=YP/(val(PPP,X)-fnTargetPP(X))
 1980        YP=val(PP,X)-fnTargetP(X)
 1990        if or{XL>X,X>XH} then ' Use Trisection
 2000           :X1=XL+(XH-XL)/3:Y1=val(P,X1)-fnTarget(X1)
 2010           :X2=XL+2*(XH-XL)/3:Y2=val(P,X2)-fnTarget(X2)
 2020           :print "tri",using(1,5),XL,X1,X2,XH
 2030           :if abs(Y1)<abs(Y2) then XL=X1:YL=Y1:else XH=X2:YH=Y2:endif
 2040           :X=(XL+XH)/2:YP=val(PP,X)-fnTargetP(X)
 2050           :YP=val(PP,X)-fnTargetP(X)
 2060        endif
 2070        Iter%+=1
 2080        if Iter%>10*DecPrecision% then return(I%)
 2090     wend
 2100     Xex(I%)=X
 2110     Yex(I%)=val(P,X)-fnTarget(X)
 2120   next I%
 2130   ' The next two extrema are almost for free
 2140   Yex(0)=val(P,Low#)-fnTarget(Low#)
 2150   Yex(Order%+1)=val(P,High#)-fnTarget(High#)
 2160   return(-1)
 2170   '
 2180   ' Assumes : RefineExtrema executed, this is not a proper test
 2190   *TestRefineExtrema
 2200    print "Extrema: i, y, x"
 2210    for I%=0 to Order%+1
 2220       print I%,using(1,DecPrecision%+4),Yex(I%),using(1,5),Xex(I%)
 2230    next I%
 2240   return
 2250   '
 2260   ' -------------- Fill the array with intervals ---------
 2270   '
 2280   ' Fill the array Interval(0..Order%+1) from X(0..Order%)
 2290   *FillInterval
 2300   local I%
 2310   Interval(0)=X(0)-Low#
 2320   for I%=1 to Order%
 2330     Interval(I%)=X(I%)-X(I%-1)
 2340   next I%
 2350   Interval(Order%+1)=High#-X(Order%)
 2360   return
 2370   '
 2380   ' --------- Fill the array with the interpolation points ---------
 2390   '
 2400   ' Fill the array X(0..Order%) from Interval(0..Order%+1)
 2410   *FillXFrom
 2420   local I%
 2430   X(0)=Low#+Interval(0)
 2440   for I%=1 to Order%
 2450     X(I%)=X(I%-1)+Interval(I%)
 2460   next I%
 2470   return
 2480   '
 2490   *TestInterval
 2500   Order%=2
 2510   Interval(0)=0.1:Interval(1)=0.2
 2520   Interval(2)=0.3:Interval(3)=0.4
 2530   gosub *FillXFrom
 2540   gosub *FillInterval
 2550   print "Expect :",0.1,0.2,0.3,0.4
 2560   print Interval(0),Interval(1),Interval(2),Interval(3)
 2570   return
 2580   '
 2590   ' -------------- Adjust intervals ----------------------
 2600   '
 2610   ' Make the intervals from Interval(0..Order%+1) larger or smaller
 2620   ' as follows from the extrema.
 2630   ' Quality of the approximation is also calculated.
 2640   ' It indicate the number of digits precision in base 'e'.
 2650   ' Negative qualitiy is worse for more negative.
 2660   *AdjustIntervals
 2670   local I%
 2680   local Average,T
 2690   local QH:' QL is global!
 2700   '
 2710   Average=0.
 2720   for I%=0 to Order%+1
 2730      Average+=abs(Yex(I%))
 2740   next I%
 2750   Average/=Order%+2
 2760   '
 2770   ' Adjust proportion, at the expense of absolute value
 2780   QL=10000000000.
 2790   QH=0.
 2800   for I%=0 to Order%+1
 2810     T=abs(Yex(I%))
 2820     QL=min(T,QL):QH=max(T,QH)
 2830     Correction=(T/Average)^Fudge
 2840     Interval(I%)/=Correction
 2850   next I%
 2860   '
 2870   Quality=-log((QH-QL)/QL)/log(2)
 2880   gosub *NormalizeIntervals
 2890   return
 2900   '
 2910   ' -------------- Normalize Intervals ----------------------
 2920   '
 2930   *NormalizeIntervals
 2940   local Correction,Average,T
 2950   Correction=0.
 2960   for I%=0 to Order%+1
 2970     Correction+=Interval(I%)
 2980   next I%
 2990   Correction=(High#-Low#)/Correction
 3000   for I%=0 to Order%+1
 3010     Interval(I%)*=Correction
 3020   next I%
 3030   return
 3040   '
 3050   *TestAdjust
 3060   gosub *TestInterval
 3070   Yex(0)=-5:Yex(1)=20:Yex(2)=-10:Yex(3)=10
 3080   print "Expect interchange of first two intervals"
 3090   gosub *AdjustIntervals
 3100   print Interval(0),Interval(1),Interval(2),Interval(3)
 3110   return
 3120   ' -------------- One more Interval --------------------------
 3130   '
 3140   ' One interval is added in the middle of the range in order to
 3150   ' make it posssible to attain more precision.
 3160   *OneMoreInterval
 3170   local I%
 3180   gosub *FillInterval
 3190   Order%+=1
 3200   if Order%>MaxOrder% then print "Run out of Orders":end endif
 3210   '
 3220   for I%=Order%+1 to Order%\2 step -1
 3230     Interval(I%)=Interval(I%-1)
 3240   next I%
 3250   '
 3260   gosub *NormalizeIntervals
 3270   gosub *FillXFrom
 3280   gosub *InitXex
 3290   return
 3300   '
 3310   '
 3320   ' -------------- Interpolate ---------------------------
 3330   '
 3340   ' Interpolate finds the polynomial P that agrees with fnTarget
 3350   ' at the points X(0..Order%) and its derivatives PP and PPP
 3360   *Interpolate
 3370   local I%
 3380   '
 3390   ' Fill Y()
 3400   for I%=0 to Order%
 3410    Y(I%)=fnTarget(X(I%))
 3420   next I%
 3430   '
 3440   ' Find P
 3450   gosub *Polcoe
 3460   for I%=0 to Order%
 3470      if abs(Y(I%)-val(P,X(I%)))>EpsD then
 3480      :print "Increase Precision. Type 'edit *precision'"
 3490      :edit *Precision
 3500      :end
 3510      endif
 3520   next I%
 3530   '
 3540   ' Derivatives
 3550   PP=diff(P)
 3560   PPP=diff(PP)
 3570   return
 3580   '
 3590   *TestInterpolate
 3600   local I%
 3610   Order%=2
 3620   X(0)=1:X(1)=2.5:X(2)=3:
 3630   gosub *Interpolate
 3640   print "Expect Zero 's"
 3650   for I%=0 to Order%
 3660   print val(P,X(I%))-fnTarget(X(I%))
 3670   next I%
 3680   return
 3690   '
 3700   ' -------------- Adjust X Coordinates -----------------
 3710   '
 3720   ' The X coordinates are adjusted according to the extrema that
 3730   ' have been found, via an excursion to an interval representation.
 3740   *AdjustX
 3750   gosub *FillInterval
 3760   gosub *AdjustIntervals
 3770   gosub *FillXFrom
 3780   return
 3790   '
 3800   ' -------------- save the results ---------------------
 3810   '
 3820   ' Auxiliary routine for SaveResult:
 3830   ' Output the coefficients to file #1 in decimal
 3840   *SaveInDecimal
 3850      print #1,"Coefficients in decimal     :",DecPrecision%;" digits + 2 spare digits"
 3860      for I%=0 to Order%
 3870        X=coeff(P,I%)
 3880        print #1,using(3,0),I%,using(3,0),X;".";
 3885         X=abs(X):' take care of negative X!
 3890        for J%=1 to DecPrecision%+2
 3900           X-=int(X)
 3910           X*=10
 3920           print #1,chr(asc("0")+int(X));
 3930        next J%
 3940        print #1
 3950      next I%
 3960   return
 3970   '
 3980   ' Auxiliary routine for SaveResult:
 3990   ' Output the coefficients to file #1 in hexadecimal
 4000   *SaveInHexadecimal
 4010      print P
 4020      print #1,"Coefficients in hexadecimal :",Binprecision%\4;" ";res;"/4 digits + 2 spare digits"
 4025      print #1,"Digits before the point in decimal (sorry)"
 4030      for I%=0 to Order%
 4040         X=coeff(P,I%)
 4050         print #1,using(3,0),I%,using(3,0),X;".";
 4055         X=abs(X):' take care of negative X!
 4060         for J%=1 to Binprecision%\4+3
 4070            X-=int(X)
 4080            X*=16
 4090            if int(X)<10 then
 4100               :print #1,chr(asc("0")+int(X));
 4110            :else
 4120               :print #1,chr(asc("A")+int(X)-10);
 4130            :endif
 4140         next J%
 4145      print #1
 4150      next I%
 4170   return
 4180   '
 4190   *SaveResult
 4200   OutFile=fnTargetName:+".COF"
 4210   print SavingResultsInF:'",Outfile,"'"
 4220   open OutFile for append as #1
 4230   print #1,"Approximation of            :",fnTargetName
 4235   print #1,"Interval                    :",Low#,High#
 4240   print #1,"Order of approximation      :",Order%
 4250   print #1,"Maximum absolute deviation  :",Eps#
 4260   print #1,"Quality of minmax (bits)    :",QualityWanted%
 4270   if HexOutput% then
 4280      :gosub *SaveInHexadecimal
 4290   :else
 4300      :gosub *SaveInDecimal
 4310   :endif
 4320   print #1,P
 4330   close #1
 4340   return
 4350   '
 4360   '
 4370   *Main:' Label to jump over the subroutines
 4380   ' ###################################################
 4390   ' ################## MAIN PROGRAM ###################
 4400   ' ###################################################
 4410   '
 4420   gosub *InitXY::' Filling in the Y's is superfluous
 4430   gosub *InitXex
 4440   gosub *Interpolate
 4450   Result%=fnRefineExtrema
 4460   if Result%<>-1 then print "No convergence":end
 4470   gosub *AdjustX
 4480   print "Quality:",using(4,1),Quality
 4490   gosub *TestRefineExtrema
 4500   'print "More":Dummy=input$(-1)
 4510   if and{QL>Eps#,Quality>0} then gosub *OneMoreInterval
 4520   if Quality<QualityWanted% then goto 4440
 4530   print "Quality:",using(3,1),Quality
 4540   gosub *SaveResult