Posts: 1,755
Threads: 112
Joined: Jan 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 SeptemberOctober issue of Datafile (V24 N5), namely:
HP71B Minimax Polynomial Fit
This is a 12page article including two programs, namely a 3liner to compute the
LeastSquares Nthdegree polynomial which best fits a given set of (x,y) datapoints
in the leastsquare sense, and a 50line program to compute the minimax Nthdegree
polynomial fit to an arbitrary set of datapoints (x,y), which can be entered from the
keyboard, automatically generated by evaluating an arbitrary userprovided 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 5thdegree 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 (xcoordinate 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.
Posts: 2,761
Threads: 100
Joined: Jul 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.
Posts: 1,755
Threads: 112
Joined: Jan 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 <= 1E6.
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 7thdegree polynomial:
P(x) = 0.99999661 * x  0.16664831 * x^{3} + 0.00830636 * x^{5}
 0.00018365 * x^{7}
with guaranteed absolute maximum error less than 1E6 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 nonzero 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),RESSIN(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 < 1E6, 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 7thdegree, oddpowersonly, minimax polynomial which has a guaranteed absolute maximum error e = 3.2E11 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.
Posts: 2,761
Threads: 100
Joined: Jul 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 7thdegree polynomial:
P(x) = 0.99999661 * x  0.16664831 * x^{3} + 0.00830636 * x^{5}  0.00018365 * x^{7}
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 * x^{3} + 0.00830636 * x^{5}  0.00018365 * x^{7}
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.
Posts: 1,792
Threads: 62
Joined: Jan 2005
Valentin 
Looks like an impressive effort! However, I'll have to delve into my HP71B 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 HP71B Surveying Pac manual to go with the ROM...
 KS
Posts: 1,755
Threads: 112
Joined: Jan 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 adhoc 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 crystalclear.
To save you the trouble of reading it, basically you transform the problem of finding a 7thdegree minimax polynomial for y=sin(x) in [0, Pi/4], say, for the related problem of finding a 2nddegree 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 specialform 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 (1E6, say) instead of 0, or else include a userdefined 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.
Posts: 1,755
Threads: 112
Joined: Jan 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 "HP71B" it goes without saying that it includes a Math ROM, because if not, it is *not* the real HP71B 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 HP71B did include most functionality now present in the Math ROM, including complex numbers definition and functionality, of which some scraps still remain in a barebones HP71B (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 HP71B deserves that name if and only if it includes a Math ROM.
Else, you should change its name, perhaps "HP71BL" for "Light" will do, though "HP71BM" would be better, or even "HP71BFU" ... ;)
BTW, don't miss my upcoming S&SCM #11, due today ! :)
Best regards from V.
Posts: 2,761
Threads: 100
Joined: Jul 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.58E11 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.76E11 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,32*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 5E10 and 7E11, respectively. A remarkable achievement for both programmer and HP71B! 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.6E11 and 9.8E11 respectively,
but I am curious to see how they would behave in a HP12CP. 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)*(34*sin^2(a/3))

Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Gerson:
Gerson posted:
"The slight difference observed is due, of course, to different machine precisions. Anyway, the maximum absolute errors are 5E10 and 7E11, respectively. A
remarkable achievement for both programmer and HP71B! 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 HP71 is a real pleasure but certainly this kind of program benefits a lot from being run in an emulator (say JeanFrançois Garnier's Emu71 or Hrastprogrammer's HP71X), 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 HP71B/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.
