Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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"?


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.


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.

share|improve this question
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 at 12:43
Is your goal optimizing an expression (for speed) for floating point evaluation? If so, this might be interesting. –  Szabolcs Feb 14 at 18:23
add comment

2 Answers

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.

share|improve this answer
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 Tunyasuvunakool Feb 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 at 14:15
add comment

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.

share|improve this answer
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.