# Can I tell Simplify to prefer nested polynomials?

I need to export expressions for numerical evaluations in a C codebase. Simplify often does a good job for this, and polynomial factorisation goes a long way.

However, occasionally I would come across things like

8 + 13 r^2 + 11 r^4 + 5 r^6 + r^8


as a subexpression. Obviously for floating point evaluation it is preferable to write this as

8 + r^2 (13 + r^2 (11 + r^2 (5 + r^2)))


Given a polynomial it is easy enough to write a function to do this conversion. The question is: how do I tell Mathematica that, after Simplify is done with its default rules, it should find all (sufficiently complicated by some metric) polynomial subexpression and turn each of them into the above "nested form"?

EDIT:

I wrote this question up in a hurry and perhaps didn't quite make the point clear. I'll admit I wasn't aware of HornerForm, but that sort of really isn't the issue. As mentioned before I've quickly cooked up a bit of code to perform the same basic function already so I perhaps neglected to check if there's a built-in one.

The real issue here is that HornerForm would undo any factorisation, that is, Simplify may well have produced a factor of

(1 + 3 r^3 + r^7)^2


which HornerForm would expand into

1 + r^3 (6 + r^3 (9 + r (2 + r^3 (6 + r^4))))


whereas I would rather prefer

(1 + 3 r^3 + r^7)^2 -> (HornerForm[1 + 3 r^3 + r^7, r])^2


which gives (1 + r^3 (3 + r^4))^2.

Similarly,

HornerForm[Sqrt[1 + 3 r^3 + r^7], r]


would just fail completely, where my desired output is Sqrt[1 + r^3 (3 + r^4)].

I suppose what I'm really looking for is to put each unfactorisable polynomial in the expression, on which Simplify was previously applied, through HornerForm.

Apologies for any confusion.

-
HornerForm[8 + 13 r^2 + 11 r^4 + 5 r^6 + r^8], see e.g. HornerForm of polynomials in terms of E^(i x) –  Artes Feb 14 '14 at 12:43
Is your goal optimizing an expression (for speed) for floating point evaluation? If so, this might be interesting. –  Szabolcs Feb 14 '14 at 18:23

Having now thought about it a bit further, the following seems to do the trick:

PolynomialSimplify[expr_, var_] :=
Simplify[expr] /. poly___Plus ? (PolynomialQ[#, var] &) :> HornerForm[poly, var]


There is probably a better way to achieve this though.

-

What about try to estimate the expense for numerical calculation using an estimatiExpense function, and use it as the ComplexityFunction?

estimateExpense[t_] :=
Apply[Times, 10^# & /@ Cases[t, _^n_?NumericQ :> n, Infinity]] +
10*Count[t, Times, Infinity] + LeafCount[t]

FullSimplify[8 + 13 r^2 + 11 r^4 + 5 r^6 + r^8,
ComplexityFunction -> estimateExpense]


(1 + r^3 (2 + r^4))^2

FullSimplify[8 + 13 r^2 + 11 r^4 + 5 r^6 + r^8,
ComplexityFunction -> estimateExpense]


8 + r^2 (13 + r^2 (11 + r^2 (5 + r^2)))

I just give a pretty random estimateExpense function. I am sure one can cook up more scientific ones based on real computational cost.

-
Thanks for this, it does seem to be an interesting idea. I should note, however, that my real concern here is not so much performance as floating point precision. If r is orders of magnitude away from 1 then 8 + 13 r^2 + 11 r^4 + 5 r^6 + r^8 is obviously a rather bad thing to be doing. –  Saran Feb 14 '14 at 14:12
I see. Thanks for explanation. I guess you may modify the ComplexityFunction to achieve that. I can understand your concern but I am not expert enough to write such a function to estimate numerical precision. –  Yi Wang Feb 14 '14 at 14:15