# Context

The following is an issue I have come across many times and still do not have a satisfactory solution for.

The context of it is this. Let's say we have an equation, a linear ODE for the function f of the variable x:

equation = c2[x] f''[x] + c1[x] f'[x] + c0[x] f[x] + d0[x] == 0


each of the coefficients depends on x, but they can also depend on another function(s), let's say g[x], or on other parameter(s), let's say p.

We want to solve this equation numerically, therefore we discretize the variable x onto some grid, let's say

grid = {1.,2.,3.,4.};


The derivatives will of course become matrices, but this question is only about the coefficients, which become vectors.

At this point an example is probably useful. The equation above is not important, we can focus on one of the coefficients, which could for example take the form:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/ (g - x^3)) g^2 x^2;


(Note: I hope this is a sufficiently complicated example, but any simplifying structure here is unintended, the dependence on the variable x, the function(s) g and the parameters p can be arbitrary)

Of course x is known ahead of time (grid above), but p and g are not. Imagine a sequence of computations running many times in a loop, that contains the pieces:

(* ... *)
p = 3.; (* actually some computation *)
g = {0.714, 0.429,0.734, 0.921}; (* some other computation *)
cComputed = Block[{x = grid}, c]
(* use this vector c to construct and solve the discretized equation *)


{-47.6352, 78.8935, 1.77501*10^15, 3.19207*10^55}

Done this way, at every pass through the loop everything in c is computed from scratch. However since we know x ahead of time, we can already do a big part of the computations:

cPartlyComputed = g^2({1., 4., 9., 16.} Sinh[{0.5, 8., 40.5, 128.} g + {0.5, 2., 4.5, 8.}]
+ p^3  {1., 8., 27., 64.}/(g - {1., 8., 27., 64.}));


Note that I have not merely inserted {1.,2.,3.,4.} for x, but also computed x^2, 1/2 x^2 and 1/2 x^4.

## speed comparison

A quick test to show that this gives a significant speedup:

without = Function[{g, p, x}, g^2 x^2 ((p^3 x)/(g - x^3) + Sinh[1/2 (x^2 + g x^4)])];
with = Function[{g, p}, g^2 ({1., 4., 9., 16.} Sinh[g {0.5, 8., 40.5, 128.} + {0.5, 2., 4.5, 8.}] + (p^3 {1., 8., 27., 64.})/(g + {-1., -8., -27., -64.}))];

gval = RandomReal[{}, {4}];
pval = 3.;

Do[without[gval, pval, grid], {10^5}] // AbsoluteTiming//Print;
Do[with[gval, pval], {10^5}] // AbsoluteTiming//Print;


1.35754

0.719002

About a factor two speedup.

# The question

This is the question: how can we create the expression cPartlyComputed above automatically?

## naive solution

The naive solution illustrates the problem:

cTry = c /. x -> grid;


{1. g^2 ((1. p^3)/(-1. + g) + Sinh[1/2 (1. + 1. g)]), 4. g^2 ((2. p^3)/(-8. + g) + Sinh[1/2 (4. + 16. g)]), 9. g^2 ((3. p^3)/(-27. + g) + Sinh[1/2 (9. + 81. g)]), 16. g^2 ((4. p^3)/(-64. + g) + Sinh[1/2 (16. + 256. g)])}

The unknown vector g has entered into the elements of a list. When we have computed g and insert it into this, we will get a matrix, not the vector that we want.

# Best attempt so far

Since the essence of the problem seems to be the Listable property of Plus and Times when applied to symbolic arguments, the idea was to temporarily remove this property, only for symbolic arguments, or create our own Plus and Times identical to the built-in ones except for this property.

I tried the latter:

ClearAll[CircleTimes, CirclePlus]
Attributes[CirclePlus] = Attributes[CircleTimes] = Complement[Attributes[Plus], {Listable, Protected}];
Verbatim[CircleTimes][a_] := a;
Verbatim[CirclePlus][a_] := a;
numericOrList[a_] := NumericQ[a] || Head[a] === List;
a_\[CircleTimes]b_ := a*b /; numericOrList[a] && numericOrList[b];
a_\[CirclePlus]b_ := a + b /; numericOrList[a] && numericOrList[b];
a_\[CircleTimes](b_\[CirclePlus]c_) :=
a\[CircleTimes]b\[CirclePlus]a\[CircleTimes]c /; numericOrList[a]
a_\[CircleTimes]x1_\[CirclePlus]a_\[CircleTimes]x2_ :=
a\[CircleTimes](x1 + x2) /; numericOrList[x1] && numericOrList[x2];

evaluate[expression_, var_, grid_] := Hold[Evaluate[
expression /. {Plus -> CirclePlus, Times -> CircleTimes} /. var -> grid
]] /. {CirclePlus -> Plus, CircleTimes -> Times}


## some explanation:

CircleTimes and CirclePlus are undefined built-in functions that display nicely. As mentioned we give them all the properties of Times and Plus except Listable (and Protected of course).

We check their arguments with numericOrList which is True for pure numbers or for lists, which in our context will always be lists of pure numbers. On such arguments we do want the Listable property, so all the definitions are to ensure that they reduce to Times and Plus in these cases.

The Verbatim is necessary to prevent an infinite recursion, see Carl Woll's answer to this question.

Finally the function evaluate uses these CircleTimes and CirclePlus together with the naive answer above, holding the result ready to be used in a later computation.

## performance

This seems to work:

evaluate[c,x,grid]


Hold[g^2 ({1., 4., 9., 16.} Sinh[ g {0.5, 8., 40.5, 128.} + {0.5, 2., 4.5, 8.}] + ( p^3 {1., 8., 27., 64.})/(g + {-1., -8., -27., -64.}))]

identical to the cPartlyComputed above that I did by hand.

However it is prohibitively slow, presumably because of all the pattern matching. On this other still very simple example:

c2=-5832. - 432. x^2 - 1728. f x^3 - 12. (32. + 72. d1f + 144. f^2 - 27. p^4) x^4
- 24. f (32. + 72. d1f - 27. p^4) x^5


it creates something that seems optimal:

evaluate[c2,x,grid]


Hold[d1f {-864., -13824., -69984., -221184.} + f (d1f {-1728., -55296., -419904., -1.76947*10^6} + p^4 {648., 20736., 157464., 663552.} + {-768., -24576., -186624., -786432.}) + f {-1728., -13824., -46656., -110592.} + f^2 {-1728., -27648., -139968., -442368.} + p^4 {324., 5184., 26244., 82944.} + {-6648., -13704., -40824., -111048.}]

but takes 34 seconds to do so on my laptop (I don't understand why, the pattern matching is probably slow but that slow?) This is so slow that I cannot even test if it would at least work in a more general case.

# Requirements of the solution

So what I am looking for is this:

1. A function evaluate[expression_,var_,grid_] that replaces the symbolic variable var with the numerical list grid in expression, without messing it up by taking any other symbols in expression into the lists.
2. expression can be either a single expression or a (nested) list of expressions, each one being an arbitrary nonlinear function of one or several symbols. A more realistic example on which it should work is given below.
3. The result should contain no, or as few as possible, operations that can be carried out with just the knowledge of var.
4. While for a given problem this has to be done only once and the result will be used many times, it is still important that it be reasonably fast.
5. If the equation, and the functions and parameters, are at machine precision, then inserting packed arrays for the functions should ideally result in a packed array as output as well.
6. Ideally the code should be short and easy to understand.

It doesn't have to be related to my attempt above, perhaps there is a way to "hack" the built in Times and Plus to behave like this in a much faster way, or something else entirely.

## An actual example to test on

I've uploaded a more realistic case on which it should work. It depends on the functions f,g,h and their derivatives, which have been rewritten as f'[x] -> d1f etc( if it's easier this can of course be undone) and on the scalar parameter p.

# Update

There have been two answers so far, by Carl Woll and xzczd, which are both close to what I'm looking for but not quite. A quick comparison on the example linked above: with a grid size of 4, Carl Woll's solution takes about 0.013 seconds to produce something with a LeafCount of 4573, while xzczd's answer takes about 0.25 seconds to produce something with a LeafCount of 3409*. So the former is a lot faster but also computes less. Computations with the result of the latter are then about 25% faster than with the former.

However both leave computations undone. With Carl Woll's answer I think it's because the rules aren't repeated, resulting in expressions like f {-1728., -13824., -46656., -110592.} + f {0., -51104., -61502., -3.270*^6} that are not evaluated. (it is also becoming harder to read the code) With xzczd's answer the problem is that additions of numbers or a number and a list aren't computed, because for some reason including that it becomes incredibly slow.

Therefore I will set a bounty for a solution that does do really all the computations in a reasonable time (I don't think it should take more than 0.5 seconds for the example linked above, possibly much shorter). And if it's still easily readable a year from now that would be great!

[* Actually to obtain this I've added to his last two patterns also the ones with patt2 instead of patt]

• @xzczd Hmm that is strange, I expected it wouldn't matter much when compiled, but how can it be slower? You're doing less. Even if not compiled the first is faster in this case. I am not sure if it is true in an actual application though, as in the linked example, and for grid sizes up to several hundred, and at arbitrary precision (didn't think this last thing was relevant but it might be now). – Jansen Oct 23 '18 at 11:58
• Oops, my previous (now deleted) test is wrong. The cPartlyComputed can't be defined that way. The conclusion doesn't change much though: Clear[g, p]; test1 = Compile[{{g, _Real, 1}, p, {x, _Real, 1}}, #] &@c; test2 = Hold@Compile[{{g, _Real, 1}, p}, c] /. OwnValues@c /. x -> grid // ReleaseHold; pvalue = 3.; gvalue = {0.714, 0.429, 0.734, 0.921}; Do[test1[gvalue, pvalue, grid], {10^6}]; // AbsoluteTiming // Print; (*{3.192,Null}*) Do[test2[gvalue, pvalue], {10^6}]; // AbsoluteTiming // Print;(*{3.193,Null}*) As we can see the influence of inlining of grid is quite limited. – xzczd Oct 23 '18 at 12:11
• I cannot access the uploaded file, but according to my personal limited experience, the inlining of constant almost never influence the speed much. Still, I suspect this is a XY problem and the timing should be shortened in other ways. – xzczd Oct 23 '18 at 12:18
• @xzczd This is not a fair comparison though, in your test2 you have only replaced x with the grid, but not pre-computed for instance x^2, which is what I'm looking to do. It should definitely be faster to precompute as much as possible. It might depend on the expression by how much of course, but in my experience it's significant. – Jansen Oct 23 '18 at 12:26
• It does help: test3 = Hold@Compile[{{g, _Real, 1}, p}, c] /. OwnValues@c /. x :> grid /. HoldPattern[grid^a_.] :> RuleCondition[grid^a] // ReleaseHold; Do[test3[gvalue, pvalue], {10^6}]; // AbsoluteTiming // Print;(* {2.76234,Null} *) But still, I doubt if this will have a significant influence in real situation. – xzczd Oct 23 '18 at 12:35

## 4 Answers

You can set x to be numeric, and then search for numeric subexpressions to be replaced. For example:

evaluate[params_, expr_, v_, grid_] := Module[{x},
NumericQ[x]=True;
ReplaceAll[
Activate @ Inactive[Function][params, expr /. v->x],
h_?NumericQ :> RuleCondition[h /. x->grid]
]
]


For your first example:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/ (g - x^3)) g^2 x^2;
evaluate[{g, p}, c, x, {1, 2, 3, 4}]


Function[{g, p}, g^2 {1, 4, 9, 16} ((p^3 {1, 2, 3, 4})/(g + {-1, -8, -27, -64}) + Sinh[1/2 ({1, 4, 9, 16} + g {1, 16, 81, 256})])]

For your second example:

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4 -24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f, d1f, p}, c2, x, {1, 2, 3, 4}]


Function[{f, d1f, p}, -5832. + {-432., -1728., -3888., -6912.} - 1728. f {1, 8, 27, 64} - 12. (32. + 72. d1f + 144. f^2 - 27. p^4) {1, 16, 81, 256} - 24. f (32. + 72. d1f - 27. p^4) {1, 32, 243, 1024}]

You'll notice that multiplicands and addends that are both numeric don't necessarily get combined in this approach, but multiplying and adding packed arrays is really fast, and is unlikely to be where most of the computational time is spent.

Still, it is not too difficult to combine these multiplicands and addends, and the above code could be modified to do so:

evaluate[params_, expr_, v_, grid_] := Module[{x, foo, tmp, rule1, rule2},
SetAttributes[foo, NumericFunction];
NumericQ[x]=True;
rule1 = RuleDelayed[
HoldPattern[Plus[h__]],
Total @ KeyValueMap[
If[#1, foo[Plus@@#2], Plus@@(#2/.rule1)]&,
GroupBy[{h},NumericQ]
]
];
rule2 = h_?NumericQ z_:1 :> If[z===1, foo[h], foo[h] ReplaceAll[rule2] @ z];
Activate[
Inactive[Function][params, expr] /.
v->x /.
rule1 /.
rule2
] /. foo[z_] :> RuleCondition[z /. {foo -> Identity, x -> grid}]
]


Your second example:

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4-24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f,d1f,p}, c2, x, {1,2,3,4}]


Function[{f, d1f, p}, f {-1728., -13824., -46656., -110592.} + (p^4 (-27.) + 32. + d1f 72. + f^2 144.) {-12., -192., -972., -3072.} + f (p^4 (-27.) + 32. + d1f 72.) {-24., -768., -5832., -24576.} + {-6264., -7560., -9720., -12744.}]

Update

My first couple answers traverse the expression in a top-down approach, isolating the x-dependent parts along the way. Your CirclePlus/CircleTimes approach (as well as @xzxcd) traverse the expression in a bottom-up approach.

Now, there are two things that my approach didn't do that you wanted:

1. Convert $$f(x) (a + b + \ldots)$$ to $$a f(x) + b f(x) + \ldots$$

2. Convert $$f(x) a + g(x) a + h(x) b + \dots$$ to $$(f(x) + g(x)) a + h(x) b + \ldots$$

This is difficult to do in a single-pass top-down approach, but straightforward in single-pass bottom-up approach. So, rather than processing the expression with two or more passes, I will switch to using a bottom-up approach as well.

Now, the reason for the slowness of the addition rules in the other rules is because they rely on the flat and orderless attributes of CirclePlus. Consider the following definition:

ClearAll[CirclePlus]
SetAttributes[CirclePlus, {Flat,Orderless}]

a_?NumericQ ⊕ b_?NumericQ := a+b


Let's count how many times NumericQ is called when evaluating 1 ⊕ x ⊕ 2

i = 0;
TracePrint[1 ⊕ x ⊕ 2,_NumericQ, TraceAction :> (i++&)]
i


3⊕x

70

We see that NumericQ is called 70 times! This is why including such binary rules causes the other answers to be so slow. Rather than using binary definitions and relying on the Flat and Orderless attributes, it is much faster to use definitions that don't rely on these attributes. Here is an example:

Clear[CirclePlus]

Verbatim[CirclePlus][a__] /; Count[{a}, _?NumericQ]>1 := CirclePlus @@ GroupBy[{a}, NumericQ, Apply[Plus]]


Then:

i = 0;
TracePrint[CirclePlus[1, x, 2], _NumericQ, TraceAction -> (i++&)]
i


3⊕x

7

Much better. Note that the excessive NumericQ calls are exacerbated by the Orderless attribute. I will first present a broad overview of my updated answer.

1. Replace Times and Plus with CircleTimes and CirclePlus.

2. Have CircleTimes separate numeric parts (assuming x is numeric) from non-numeric parts, wrapping the numeric parts in a n wrapper, and the non-numeric parts in a s wrapper.

3. Define distributive and product rules for the n wrapper.

4. Have CirclePlus factor non-numeric parts.

Here is the code:

ClearAll[CircleTimes, CirclePlus, nWrapper]

evaluate[params_, expr_, v_, grid_] := Module[{x, res},
NumericQ[x]=True;
res = Inactive[Function][params, expr] /. v->x /. {Times->CircleTimes, Plus->CirclePlus};
ReplaceAll[
Activate[res /. sWrapper -> Identity],
{
nWrapper[a_] :> RuleCondition @ ReplaceAll[a, {x -> grid, nWrapper -> Identity}],
n_?NumericQ :> RuleCondition @ ReplaceAll[n, x->grid]
}
]
]

CircleTimes[a__] := Times @@ KeyValueMap[
If[#, nWrapper[#2], sWrapper[#2]]&,
GroupBy[{a}, NumericQ, Apply[Times]]
]

CirclePlus[a__] := Plus @@ KeyValueMap[
Times[#1, nWrapper[#2]]&,
GroupBy[
Flatten @ Replace[{a}, Verbatim[Plus][b__] :> {b}, {1}],
sPart -> nPart,
Apply[Plus]
]
]

sPart[_?NumericQ] := 1
sPart[u_sWrapper] := u
sPart[t_Times] := Replace[t, _?NumericQ->1, {1}]
sPart[u_] := sWrapper[u]

nPart[t_?NumericQ] := nWrapper[t]
nPart[u_sWrapper] := 1
nPart[t_Times] := Replace[t, Except[_?NumericQ]->1, {1}]
nPart[u_] := 1

SetAttributes[nWrapper, {Flat, NumericFunction}]
nWrapper[1] = 1;
nWrapper /: f_nWrapper a_Plus := Distribute @ Unevaluated[f a]
nWrapper /: f_nWrapper sWrapper[a_Plus b_.] := Distribute @ Unevaluated[f a] If[b===1, 1, sWrapper[b]]
nWrapper /: Verbatim[nWrapper][f_] Verbatim[nWrapper][g_] := nWrapper[f g]


In your examples, numbers and lists now distribute, producing your expected outputs:

c = (Sinh[1/2 (g x^4+x^2)]+p^3 x/(g-x^3)) g^2 x^2;
evaluate[{g,p}, c, x, {1, 2, 3, 4}]


Function[{g, p}, g^2 ((p^3 {1, 8, 27, 64})/( g + {-1, -8, -27, -64}) + {1, 4, 9, 16} Sinh[{1/2, 2, 9/2, 8} + g {1/2, 8, 81/2, 128}])]

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4 -24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f, d1f, p}, c2, x, {1, 2, 3, 4}]


Function[{f, d1f, p}, f {-1728., -13824., -46656., -110592.} + f^2 {-1728., -27648., -139968., -442368.} + d1f {-864., -13824., -69984., -221184.} + p^4 {324., 5184., 26244., 82944.} + f (d1f {-1728., -55296., -419904., -1.76947*10^6} + {-768., -24576., \ -186624., -786432.} + p^4 {648., 20736., 157464., 663552.}) + {-6648., -13704., -40824., -111048.}]

Here's a view of the "collected" expression before x is replaced with a list:

evaluate[{f, d1f, p}, c2, x, x]


Function[{f, d1f, p}, f (-1728. x^3) + f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + f (d1f (-1728. x^5) - 768. x^5 + p^4 (648. x^5)) + (-5832. - 432. x^2 - 384. x^4)]

Notice how all of the x-dependent coefficients are grouped together. Basically, the function evaluate acts like Collect, but without restructuring. For instance, compare the above to:

evaluate[{f, d1f, p}, Expand @ c2, x, x]


Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + d1f f (-1728. x^5) + f p^4 (648. x^5) + (-5832. - 432. x^2 - 384. x^4) + f (-1728. x^3 - 768. x^5)]

The function evaluate does not do any factoring or expanding of the parameters. On the other hand, if you used Collect, you would get a restructuring dependent on the order of parameters used:

Activate[
Inactive[Function][{f, d1f, p}, Collect[c2, {f, d1f, p}, w]]
] /. w[s_] -> s


Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + (-5832. - 432. x^2 - 384. x^4) + f (d1f (-1728. x^5) + p^4 (648. x^5) + (-1728. x^3 - 768. x^5))]

vs:

Activate[
Inactive[Function][{f, d1f, p}, Collect[c2, {p, d1f, f}, w]]
] /. w[s_] -> s


Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4 + f (-1728. x^5)) + p^4 (324. x^4 + f (648. x^5)) + (-5832. - 432. x^2 - 384. x^4) + f (-1728. x^3 - 768. x^5)]

However, for your big example, the kind of restructuring done by Collect would be bad, as it will expand all powers, creating many more terms to evaluate.

At any rate, let's see how evaluate does on your big example:

evaluate[
{d1f,d1g,d1h,d2f,d2g,d2h,f,g,h,p},
exampleCoefficients,
x,
{1,2,3,4}
]; //AbsoluteTiming


{0.048297, Null}

• Wow, didn't know NumericQ[x] = True; NumericQ[1 + x] gives True! – xzczd Oct 23 '18 at 17:47
• Thanks, that is very fast, but as you say it leaves a lot of computations undone. On the bigger example I linked to it takes only about 0.01 seconds but the result is still full of multiplications and additions of numerical lists and factors. How would you modify it to resolve these too? – Jansen Oct 24 '18 at 9:56
• @Jansen See update. – Carl Woll Oct 24 '18 at 17:03
• Thanks, but now it is very slow. Even on a piece of the linked expression it doesn't manage. It might be for the same reason as @xzcd's answer, because if I remove your rule1 it becomes very fast (but of course incomplete). – Jansen Oct 25 '18 at 12:28
• @Jansen Should be faster now. – Carl Woll Oct 25 '18 at 16:16

The matching for certain pattern e.g. 1⊕{1, 2} seems to be extremely slow, so I exclude them from the matcher:

ClearAll[CircleTimes, CirclePlus]
Attributes[CirclePlus] =
Attributes[CircleTimes] = Complement[Attributes[Plus], {Listable, Protected}];
Verbatim[CircleTimes][a_] := a;
Verbatim[CirclePlus][a_] := a;

With[{patt =(*__?NumericQ|*){__?NumericQ}, patt2 = __?NumericQ},
(a : patt)⊗(b : patt) := a b;
(a : patt)⊗(b : patt2) := a b;
(a : patt2)⊗(b : patt2) := a b;
(a : patt)⊕(b : patt) := a + b;
(*(a:patt)⊕(b:patt2):=a+b;*)
(a : patt)⊗(b_⊕c_) := a⊗b⊕a⊗c;
a_⊗(x1 : patt)⊕a_⊗(x2 : patt) := a⊗(x1 + x2);]

evaluate[expression_, var_, grid_] :=
Hold[Evaluate[
expression /. {Plus -> CirclePlus, Times -> CircleTimes} /.
var -> grid]] /. {CirclePlus -> Plus, CircleTimes -> Times}

test = evaluate[c2, x, grid]; // AbsoluteTiming
(* {0., Null} *)


Still, the result is very close to the optimized one.

• Thanks, but this still becomes very slow very fast for some reason. For instance on -5832. - 432. x^2 - 1728. f x^3 - 384. x^4 - 864. d1f x^4 - 1728. f^2 x^4 + 324. p^4 x^4 - 768. f x^5 - 1728. d1f f x^5 + 648. f p^4 x^5 - 85.3333 x^6 - 384. d1f x^6 - 432. d1f^2 x^6 + 144. p^4 x^6 it takes 4.7 seconds for me. (with a grid of size 2, only 5.1 with a grid of size 50, so clearly it is wasting a lot of time somewhere). The example is a small piece from what I linked, if you know somewhere else that I can upload it to where you can access it I'd be happy to. – Jansen Oct 23 '18 at 15:44
• @jansen You may consider Github. – xzczd Oct 23 '18 at 16:51
• Ok, I've changed the link to Github. – Jansen Oct 24 '18 at 10:30
• @Jansen Check my update. – xzczd Oct 24 '18 at 12:40
• Thanks, this works quite well, on the bigger example it takes about 0.25 seconds and gives a speedup of about a factor 2 (machine precision) or 4 (arbitrary precision). I really don't understand why a pattern with ⊕ is very slow when the exact same pattern with ⊗ is very fast? I've checked that if you remove all definitions except the one you commented out it is still the case, so it can't be an interplay between those definitions. It has to be something to do with the expression it's applied to then. – Jansen Oct 24 '18 at 14:18

You may use the Trott-Strzebonski in-place evaluation technique which I first learned from here: 19542.

With

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/(g - x^3)) g^2 x^2;
grid = {1., 2., 3., 4.};


Then

cPartlyComputed =
Function[{g, p}, Evaluate[c]] /.
{
f_[x, n_] :> With[{v = f[grid, n]}, v /; True],
x :> With[{v = grid}, v /; True]
}


Function as attribute HoldAll so it can be treated as a held expression. f_[x, n_] is getting the Power[x, _] entries and x picks up any left over x'es.

Hope this helps.

Here is a different attempt to work around the Listable attribute of Plus and Times.

doPE[c_, x_, grid_] :=
InternalInheritedBlock[{Plus, Times},
ClearAttributes[{Plus, Times}, {Flat, Listable}];
Hold @@ {
c /. x -> grid //.
{
a___ * (n_?NumericQ | n_List) * L_List :> a * Thread[n * L],
a___ + (n_?NumericQ | n_List) + L_List :> a + Thread[n + L],
a___ + f_*L_List + f_*M_List :> a + f Thread[L + M]
}
}
]


Test:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/(g - x^3)) g^2 x^2;
c2 = -5832. - 432. x^2 - 1728. f x^3 -
12. (32. + 72. d1f + 144. f^2 - 27. p^4) x^4 -
24. f (32. + 72. d1f - 27. p^4) x^5;
grid = {1., 2., 3., 4.};

doPE[c, x, grid]

doPE[c2, x, grid]

Hold[g^2 {1., 4., 9., 16.}
((p^3 {1., 2., 3., 4.})/(g + {-1., -8., -27., -64.}) +
Sinh[1/2 ({1., 4., 9., 16.} + g {1., 16., 81., 256.})])]

Hold[{-6264., -7560., -9720., -12744.} +
f {-1728., -13824., -46656., -110592.} +
f (32. + 72. d1f - 27. p^4) {-24., -768., -5832., -24576.} +
(32. + 72. d1f + 144. f^2 - 27. p^4) {-12., -192., -972., -3072.}]

• This is extremely slow though. On the linked expression it doesn't finish in 10 minutes. I think the problem may have the same origin as in @xzczd's answer: if I remove the addition rule it becomes very fast (but of course incomplete). – Jansen Oct 27 '18 at 9:28
• @Jansen Sorry, I didn't test it on anything larger. I'll see if I can make anything useful out of this. – Mr.Wizard Oct 27 '18 at 10:12
• @Jansen Please check my current answer. I don't think it's as good as others, but I hope it at least works now, and perhaps it will inspire a better solution that also uses InternalInheritedBlock. – Mr.Wizard Oct 27 '18 at 11:55