Minimax Polynomial Fit - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Minimax Polynomial Fit (/thread-79477.html) Minimax Polynomial Fit - Valentin Albillo - 10-10-2005 Hi all, Busy as always, just a quick post to let you know that a new article from mine appears published in the latest September-October issue of Datafile (V24 N5), namely: HP-71B Minimax Polynomial Fit This is a 12-page article including two programs, namely a 3-liner to compute the Least-Squares Nth-degree polynomial which best fits a given set of (x,y) datapoints in the least-square sense, and a 50-line program to compute the minimax Nth-degree polynomial fit to an arbitrary set of datapoints (x,y), which can be entered from the keyboard, automatically generated by evaluating an arbitrary user-provided function, or read/stored from/to a data file. The user can specify the degree of the polynomial or else specify instead a maximum absolute error and let the program select the lowest degree which will achieve it. Once the coefficients have been computed, the user can evaluate the resulting polynomial at any given arguments right from the keyboard (even as part of more complicated expressions), and can also scan any given interval for maximum errors in it, at whatever step size desired. For example, when asked to automatically produce the minimax polynomial of lowest degree which fits y = Gamma(x) in the interval [1,2] with maximum absolute error less than 0.00005 (i.e.: 4 correct decimal places), my program produces the 5th-degree minimax polynomial: ``` P(x) = A0 + A1*x + A2*x^2 + ... + A5*x^5 ``` with ``` A0 = 3.73459956 A1 = -6.72738793 A2 = 6.54439953 A3 = -3.37095678 A4 = 0.91923039 A5 = -0.09992158 ``` Once computed, we can use the automatic scan option to check the [1,2] interval at 0.01 steps to confirm that it has absolute maximum error app. 0.00003812 (<0.00005) at x=1.06 ``` >FNE(1, 2, 0.01) 1.06000000 (x-coordinate where the max. err. occurs) 0.00003812 (absolute maximum error in [1,2]) ``` and can evaluate it for any argument, say Ln(Pi), like this: ``` >FIX 5 @ FNF(LN(PI)) 0.93479 (P(Ln(Pi)) = predicted value of Gamma(Ln(Pi))) ``` while the correct value is ``` GAMMA(LN(PI)) -> 0.93480 ```i.e., we sure got 4 correct decimal places, more nearly five. Any and all comments are very welcome. Best regards from V. Re: Minimax Polynomial Fit - Gerson W. Barbosa - 10-10-2005 Quote: The user can specify the degree of the polynomial or else specify instead a maximum absolute error and let the program select the lowest degree which will achieve it. Hello Valentin, Thanks for letting us know of your new article. Just one question: can the user also specify polynomials with odd or even powers only? This would be particularly useful when fitting odd and even functions, as less coefficients would be needed for a given accuracy. Best regards, Gerson. Re: Minimax Polynomial Fit - Valentin Albillo - 10-11-2005 Hi, Gerson: Gerson posted: "Can the user also specify polynomials with odd or even powers only? This would be particularly useful when fitting odd and even functions, as less coefficients would be needed for a given accuracy." First of all, thanks for your interest in my article. As for your question, there's no need of a separate option to achieve that effect. You just simply specify an interval where it shows that the function is odd or even, and the minimax error condition itself will guarantee that only odd or even powers are used. For instance, let's say that you want to fit a minimax polynomial of lowest possible degree to y=sin(x) for values of x not exceeding Pi/2 and using just odd powers of x (as sin(-x) = -sin(x), so sin(x) is an odd function). Also, you want a guaranteed maximum error <= 1E-6. As I've stated, you only need to specify an interval where the 'oddness' of the functions shows, in this case [-Pi/2, Pi/2] will do. For this particular function, interval, and maximum error, my program returns the following 7th-degree polynomial: P(x) = 0.99999661 * x - 0.16664831 * x3 + 0.00830636 * x5 - 0.00018365 * x7 with guaranteed absolute maximum error less than 1E-6 in [-Pi/2, Pi/2]. As you can see, only odd powers of x were returned, all the remaining coefficients being zero within machine rounding limits so our minimax approximation only requires four non-zero coefficients. Here's a short table of sin(x) vs. P(x): ```>FIX 7 @ FOR X=-PI/2 TO PI/2 STEP PI/20 @ DISP X,SIN(X),FNF(X),RES-SIN(X) @ NEXT X X SIN(X) P(X) ERROR ------------------------------------------------------------------------------- -1.5707963 -1.0000000 -0.9999994 0.0000006 -1.4137167 -0.9876883 -0.9876887 -0.0000004 -1.2566371 -0.9510565 -0.9510560 0.0000005 -1.0995574 -0.8910065 -0.8910061 0.0000004 -0.9424778 -0.8090170 -0.8090173 -0.0000003 -0.7853982 -0.7071068 -0.7071074 -0.0000006 -0.6283185 -0.5877853 -0.5877856 -0.0000003 -0.4712389 -0.4539905 -0.4539903 0.0000002 -0.3141593 -0.3090170 -0.3090164 0.0000006 -0.1570796 -0.1564345 -0.1564340 0.0000005 0.0000000 0.0000000 0.0000000 0.0000000 0.1570796 0.1564345 0.1564340 -0.0000005 0.3141593 0.3090170 0.3090164 -0.0000006 0.4712389 0.4539905 0.4539903 -0.0000002 0.6283185 0.5877853 0.5877856 0.0000003 0.7853982 0.7071068 0.7071074 0.0000006 0.9424778 0.8090170 0.8090173 0.0000003 1.0995574 0.8910065 0.8910061 -0.0000004 1.2566371 0.9510565 0.9510560 -0.0000005 1.4137167 0.9876883 0.9876887 0.0000004 1.5707963 1.0000000 0.9999994 -0.0000006 ``` and, as you can see, the absolute error never exceeds 0.0000006 < 1E-6, so the results are accurate to six decimal places plus or minus 1 ulp. This exact same approach, but specifying the [-Pi/6, Pi/6] interval instead, yields a 7th-degree, odd-powers-only, minimax polynomial which has a guaranteed absolute maximum error e = 3.2E-11 in that interval, while using only four coefficients. Combined with range reduction, this could form the basis of a simple yet very accurate implementation of trigonometric functions in machines lacking them, as you know. Similar approaches would allow implementing most any specialized function, say Gamma(x), Lambert's W(x), etc. Best regards from V. Edited: 11 Oct 2005, 5:55 a.m. Re: Minimax Polynomial Fit - Gerson W. Barbosa - 10-11-2005 Quote: As I've stated, you only need to specify an interval where the 'oddness' of the functions shows, in this case [-Pi/2, Pi/2] will do. For this particular function, interval, and maximum error, my program returns the following 7th-degree polynomial: P(x) = 0.99999661 * x - 0.16664831 * x3 + 0.00830636 * x5 - 0.00018365 * x7 Hi Valentin, Thanks for your detailed and comprehensive reply. However, I still have another question, if you don't mind: is it possible to force the first coefficient to 1 and then recalculate the other coefficients? Just rounding the first coefficient to 1 would cause the accuracy to decrease, as we can see in the table below: ``` P(x) = x - 0.16664831 * x3 + 0.00830636 * x5 - 0.00018365 * x7 ``` ``` X SIN(X) P(X) ERROR ------------------------------------------------------------------------------- 0.0000000 0.0000000 0.0000000 0.0000000 0.1570796 0.1564345 0.1564345 0.0000001 0.3141593 0.3090170 0.3090175 0.0000005 0.4712389 0.4539905 0.4539919 0.0000014 0.6283185 0.5877853 0.5877877 0.0000024 0.7853982 0.7071068 0.7071100 0.0000032 0.9424778 0.8090170 0.8090205 0.0000035 1.0995574 0.8910065 0.8910099 0.0000033 1.2566371 0.9510565 0.9510602 0.0000037 1.4137167 0.9876883 0.9876935 0.0000051 1.5707963 1.0000000 1.0000046 0.0000046 ``` Of course the new polynomial with recalculated coefficients wouldn't be as accurate as the original one, but would be better than this. Setting the first coefficient to 1 is important to reduce the number of steps or registers necessary to implement these polynomials on a RPN calculator. Best regards, Gerson. Impressive! [Re: Minimax Polynomial Fit] - Karl Schneider - 10-11-2005 Valentin -- Looks like an impressive effort! However, I'll have to delve into my HP-71B manuals in order to try out your work. I saw no mention of this program needing a Math ROM, but I have one with the manual. Now all I need is a HP-71B Surveying Pac manual to go with the ROM... -- KS Re: Minimax Polynomial Fit - Valentin Albillo - 10-13-2005 Hi again, Gerson: Gerson wrote: "is it possible to force the first coefficient to 1 and then recalculate the other coefficients?" Yes, there are two ways of doing it: Modifying the program itself: as one of the coefficients would be fixed, you'll actually need one equation less to solve, and so one less datapoint in your tuples. A very simple ad-hoc modification to the SPOLMM subprogram, slightly affecting the part where it is building the matrix of the system to solve, and then including the fixed coefficient in the returned coefficient's array would suffice. This would be faster and most accurate, but requires modifying the program as it is now. You can read the reference you give in your own article, which I did, where it is explictily stated how to proceed to achieve exactly what you want, even giving a fully worked example to make it crystal-clear. To save you the trouble of reading it, basically you transform the problem of finding a 7th-degree minimax polynomial for y=sin(x) in [0, Pi/4], say, for the related problem of finding a 2nd-degree polynomial for y=f(x) in [0, (Pi/4)^2], where f(x) is defined thus: ``` f(x) = sin(sqrt(x))/x^(3/2)-1/x ``` this way you get a polynomial which directly translates to your desired special-form one for the original y=sin(x) function in the original [0, Pi/4] interval. You can achieve that with my program by simply selecting [D]efined function, then answering the "f(x)=" prompt like this: ``` f(x) = SIN(SQR(X))/(X*SQR(X))-1/X ``` and specifying the [0, (Pi/4)^2] interval, as stated. You'll get the results you want. Note:As f(x) is defined, it would generate an error when evaluated at x=0. You can either specify the interval using a small number (1E-6, say) instead of 0, or else include a user-defined function, like this: ``` DEF FNY(X) @ IF X=0 THEN FNY=-1/6 ELSE FNY=SIN(SQR(X))/(X*SQR(X))-1/X ``` then answer the "f(x)=" prompt with f(x)=FNY(X) Best regards from V. Edited: 13 Oct 2005, 6:37 a.m. Re: Impressive! [Re: Minimax Polynomial Fit] - Valentin Albillo - 10-13-2005 Hi, Karl: Karl posted: "I saw no mention of this program needing a Math ROM, but I have one with the manual." Thanks for your continued interest in my humble productions. Yes, I didn't mention that a Math ROM is needed, but of course it is. If you've followed my posts since long, you'll know that when I write "HP-71B" it goes without saying that it includes a Math ROM, because if not, it is *not* the real HP-71B originally designed by its wise creators, but the maimed thingy they were forced to release when it was clear the System ROMs were going to exceed the allotted 64 Kb of ROM space. The original HP-71B did include most functionality now present in the Math ROM, including complex numbers definition and functionality, of which some scraps still remain in a bare-bones HP-71B (try to enter the program line: 10 A=(2,3)+SIN((4,5)) and see it being duly accepted, no syntax complaints), and so for me an HP-71B deserves that name if and only if it includes a Math ROM. Else, you should change its name, perhaps "HP-71BL" for "Light" will do, though "HP-71BM" would be better, or even "HP-71BFU" ... ;-) BTW, don't miss my upcoming S&SCM #11, due today ! :-) Best regards from V. Re: Minimax Polynomial Fit - Gerson W. Barbosa - 10-16-2005 Hi Valentin, Thanks for your step by step explanations, they've really helped me. And thanks for reminding me to read my own reference too :-) At the time I didn't give too much importance to these details because I didn't have any minimax tool available. As I said, I used curve fitting and did some minor adjustments by hand to get the best results I could. I got a 9.58E-11 maximum error with the three coefficients I obtained for SIN this way. The coefficients I obtained using a true minimax polynomial yielded a maximum 9.76E-11 error. Not bad! (Of course, my manual solution took much more time and would be impossible to be done for higher degree polynomials). Using the procedure you showed me, I found the odd power coefficients for a 13th degree polynomial to approximate ATAN (of course, the coefficient indeces refer to the auxiliary 5th degree polynomial): ```8 DEF FNY(Y) @ IF X=0 THEN FNY=-1/3 ELSE FNY=ATAN(SQR(X))/(X*SQR(X))-1/X 9 END DEF X1,X2,N=0,3-2*SQR(2),20 . . . H= 0.00000000050 A0=-0.33333333283 A1= 0.19999978075 A2=-0.14284193118 A3= 0.11072111863 A4=-0.08629967887 A5= 0.05061006653 ``` Just for a comparison, here are the same coefficients calculated by an expensive mathematical software, using the tecnique shown in the reference: ```a0=-0.3333333328040560 a1= 0.1999997758846710 a2=-0.1428416640124080 a3= 0.1107161126831530 a4=-0.0862630685349892 a5= 0.0505192312223211 ``` The slight difference observed is due, of course, to different machine precisions. Anyway, the maximum absolute errors are 5E-10 and 7E-11, respectively. A remarkable achievement for both programmer and HP-71B! You program is a joy to use, I wish I only had a real 71B! Best regards, Gerson. P.S.: Below are my approximations to ATAN and SIN. Maximum errors observed in a spread sheet were 6.6E-11 and 9.8E-11 respectively, but I am curious to see how they would behave in a HP-12CP. Any simulator available? ```--------------------------------------------------------------------- a = -0.3333333329 b = 0.1999997763 c = -0.142841662 d = 0.1107161 e = -0.0862632 f = 0.05052 atan(x) = x*(1+x^2*(a+x^2*(b+x^2*(c+x^2*(d+x^2*(e+f*x^2)))))), [0 <= x <= sqrt(2)-1] atan(a) = 2*atan(sqr(1/(a^2)+1)-1/a) --------------------------------------------------------------------- a = -0.1666666652 b = 0.008333217 c = -0.000197282 sin(x) = x*(1+x^2*(a+x^2*(b+c*x^2))), [0 <= x <= pi/6] sin(a) = sin(a/3)*(3-4*sin^2(a/3)) --------------------------------------------------------------------- ``` Re: Minimax Polynomial Fit - Valentin Albillo - 10-17-2005 Hi, Gerson: Gerson posted: "The slight difference observed is due, of course, to different machine precisions. Anyway, the maximum absolute errors are 5E-10 and 7E-11, respectively. A remarkable achievement for both programmer and HP-71B! You program is a joy to use, I wish I only had a real 71B!" I'm glad you liked my program and found it useful. Owning a physical HP-71 is a real pleasure but certainly this kind of program benefits a lot from being run in an emulator (say Jean-François Garnier's Emu71 or Hrastprogrammer's HP-71X), because it's so math intensive. This way you get results in seconds and can try different approaches with utmost ease, making the 'fitting experience' rather enjoyable, IMHO. Also, it's a real testimony to HP-71B/Math capabilities that this kind of functionality, which usually requires advanced software to accomplish, can be implemented in a mere 50 lines of code, and moreover, providing three different input methods (including mass storage) and two modes of output. A basic implementation would surely fit in 30 lines of code or so. And it's minimaxing, no less ! :-) Best regards from V.