Dear Mr. Fateman, hi all,
> > I'd really appreciate a result like
> > A*x*(1+B*X*(1+C*X*(1+D*X))) , e.g., instead of x*(A+x*(F+x*(G+x*H))),
> > which is of the kind that horner() gives me.
> You can extract all of the coefficients, A, B, C, ... of a polynomial
> and do anything you wish with them.
> You may find hipow and bothcoef to be useful.
thank you for your hints, I'll exercise on that subject this weekend.
IMHO, the scheme having a "1" in each bracketed term is most suitable
for the human eye in order to check for potential numerical problems
caused by cancellation effects... therefore, I thought "there has to
be a command for that, I'm just too stupid to find it"...
> If you believe that your floating point numbers deserve to be believed
> to more digits when converting to rationals, you
> can set ratepsilon to some value of your choice. It is, by default,
> 2.0*10^(-8).
thank you again for that hint; the "info on Maxima" command on my
systems yields:
Maxima version: 5.24.0 // Maxima build date: 20:39 4/5/2011 //
Host type: i686-pc-mingw32
Lisp implementation type: GNU Common Lisp (GCL) // Lisp
implementation version: GCL 2.6.8
+ wxMaxima 11.04.0
on that system, the command ev(%,numer) with % referring to an exact
expanded version of the LPs generated by the recurrence relations
gives me, e.g.:
lp_20_x = 131460.6941413879*x^20 - 640449.5355606079*x^18 +
1324172.688388825*x^16
-1513340.215301514*x^14 + 1043287.572669983*x^12 -
444238.5793304443*x^10
+114889.2877578735*x^8 - 17020.63522338867*x^6
+ 1276.54764175415*x^4
-37.00138092041016*x^2 + 0.17619705200195
subst(0.5,x,%) -> lp_20_x=-0.048358381067374
subst(1.0,x,%) -> lp_20_x=1.0
if ratepsilon has its default value, the command horner(%) acting on
the series given above yields:
rat: replaced 0.17619705200195 by 2425/13763 = 0.17619705006176
rat: replaced -37.0013809204102 by -26789/724 = -37.0013812154696
rat: replaced 1276.54764175415 by 53615/42 = 1276.547619047619
rat: replaced -17020.6352233887 by -1259527/74 = -17020.6351351351
rat: replaced 114889.2877578735 by 804225/7 = 114889.2857142857
rat: replaced -444238.579330444 by -3109670/7 = -444238.571428571
rat: replaced 1043287.572669983 by 7303013/7 = 1043287.571428572
rat: replaced -1513340.21530151 by -7566701/5 = -1513340.2
rat: replaced 1324172.688388825 by 3972518/3 = 1324172.666666667
rat: replaced -640449.535560608 by -8325844/13 = -640449.538461539
rat: replaced 131460.6941413879 by 1708989/13 = 131460.6923076923
(%o53)
lp_20_x=(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(66157938890433180*x^2-322307913370583280)
+666393374505218360)-761592432798846012)+525037212143763540)-223564228557048600)
+57818334971650500)-8565679363007130)+642425944129450)-18621042338535)+88671628500)
/(503252628060)
which quite obviously causes numerical problems:
subst(0.5,x,%) -> lp_20_x=-0.048358476478776
subst(1.0,x,%) -> lp_20_x=0.9935268531601
if I set ratepsilon:2.0*10^(-18), the result is:
rat: replaced 0.17619705200195 by 46189/262144 = 0.17619705200195
rat: replaced -37.0013809204102 by -4849845/131072 =
-37.0013809204102
rat: replaced 1276.54764175415 by 334639305/262144 =
1276.54764175415
rat: replaced -17020.6352233887 by -557732175/32768 =
-17020.6352233887
rat: replaced 114889.2877578735 by 15058768725/131072 =
114889.2877578735
rat: replaced -444238.579330444 by -29113619535/65536 =
-444238.579330444
rat: replaced 1043287.572669983 by 136745788725/131072 =
1043287.572669983
rat: replaced -1513340.21530151 by -49589132175/32768 =
-1513340.21530151
rat: replaced 1324172.688388825 by 165647382454/125095 =
1324172.688388825
rat: replaced -640449.535560608 by -83945001525/131072 =
-640449.535560608
rat: replaced 131460.6941413879 by 34461632205/262144 =
131460.6941413879
(%o75)
lp_20_x=(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(x^2*(4310977880684475*x^2-21002199931539750)
+43423467426021376)-49626819915453000)+34212428881107750)-14567872942923300)
+3767553347307750)-558156051453000)+41861703858975)-1213382720550)+5778012955)/(32792903680)
and no numerical problems occur:
subst(0.5,x,%) -> lp_20_x=-0.048358381067374
subst(1.0,x,%) -> lp_20_x=1.0
BTW: a large value of ratepsilon seems not to cause any damage if
keepfloat is set to true before having horner() act on the series with
DP float coefficients.
> The paper you are presumably referring to is
> http://www.cs.berkeley.edu/~fateman/papers/horner.pdf
Most presumably due to the larger number of hits acquired since
publication, only the following file appeared in Google's first 20
results to my search request:
http://www.cs.berkeley.edu/~fateman/legendre-horner.pdf
which is dates back to Dec. 4. 2007. The URL you gave (thank you for
providing it; I'm sorry that I did not include the URL...) seems to
reference be the most recent version of the same paper. But I think
the main claim of both versions is identical.
> (x^2*(x^2*(x^2*(x^2*(x^2*(x^2*
> (x^2*(x^2*(x^2*(34461632205*x^2-167890003050)+347123925225)-396713057400)+273491577450)-
> 116454478140)+30117537450)-4461857400)+334639305)-9699690)+46189)/(262144)
>
> I don't know if you want to convert these integers to floats; they
> should all work fine in double precision.
within Maxima, the integer representation is best, of course. Even
when transferring the results obtained with Maxima to the real world
of a programming language like Fortran, Java..., DP floats could
exactly represent integers in your formula. But I cannot spot possibly
occurring cancellation effects; that only works if I use the scheme
...*(1+CF*x^n*(1+... because almost ever my xes are limited to the
range [-1;1] and thus, if one of the CFs has an absolute value larger
than 0.99 or so, a red alert is triggered (of course depending on the
signs of x and CF as well as the power of x)...
> > these problems described by Mr. Fateman were caused by "rat" rather
> > than by "horner"?
> I hope not. I redid some of the calculations and I get exactly the same
> answers either way, so
------------
2nd response
------------
> My suspicion right now is that Macsyma is using SINGLE PRECISION FLOATS
> and that the paper is basically wrong because of that. So the problem
> is not caused by RAT or by HORNER, but by Macsyma converting rationals to the nearest
> single float instead of double float, and that's why HORNER is inaccurate.
I've included the somewhat lenthy output of my Maxima system in this
mail in order to demonstrate that at least on my Maxima system, the
value of "ratepsilon" is crucial when letting horner() act on series
with non-integer coefficients. I frequently encounter this scenario
when checking "modern" results against published "ancient" ones.
> Thanks for looking at the paper.
my pleasure; I'd like to see much more "beware of ogre" signs in the
world of numerical calculation...
Kind regards
Dragan Stoikovitch