# A quartic polynomial in four variables, for numerical integration

About half year ago, I derived a integral of a piecewise polynomial, which was a complex formula. Below is one of that formula: Obviously, the above formula has relationship with four varibles $u_i,\cdots u_{i+3}$.

In the begining, I implement it in Wolfram Mathematica with the help of Compile[] technique and corresponding option like CompilationTarget -> "C". The screenshot as shown below: Recently, I would like to rewrite it to C to achieve more higher performance.

#include <math.h>
#include <stdio.h>
//calc_A()
double calc_A(int i, const double *U) {
double u0, u1, u2, u3;
double res;
u0 = U[i];
u1 = U[i+1];
u2 = U[i+2];
u3 = U[i+3];
//implementation of the explicit formula
res = 3 * pow(u1, 4) + 3 * pow(u1, 3)*(u2 - 2 * u3) +
pow(u1, 2)*(3 * pow(u2, 2) - 7 * u2*u3 + 3 * pow(u3, 2)) +
u1*u2*(3 * pow(u2, 2) - 8 * u2*u3 + 4 * pow(u3, 2)) +
pow(u2, 2)*(3 * pow(u2, 2) - 9 * u2*u3 + 8 * pow(u3, 2)) +
pow(u0, 2)*(8 * pow(u1, 2) + 3 * pow(u2, 2) + 4 * u1*(u2 - 5 * u3) - 10 * u2*u3 +
15 * pow(u3, 2)) - u0*(9 * pow(u1, 3) + pow(u1, 2)*(8 * u2 - 19 * u3) +
u1*(7 * pow(u2, 2) - 22 * u2*u3 + 10 * pow(u3, 2)) +
u2*(6 * pow(u2, 2) - 19 * u2*u3 + 20 * pow(u3, 2)));
return res;
}

int main() {
double U[] = { 0., 0., 0., 0., 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1., 1., 1. };
int i;
double res;
//calculate A~A
for (i = 2; i <= 6; i++) {
res = calc_A(i, U);
printf("%f\n", res);
}
system("pause");
return 0;
}


where, $\mathbf U=\{u_0,\cdots,u_m\}$

TEST According to my experience, the code for res is inefficient, but I don't know how to deal with it by myself.

• What? res is a variable of type double. Not sure what you mean by that question. I don't see any human optimizations here. The compiler may make some optimizations of the formulae but that should be nothing for you to worry about. – Martin York Jan 12 '17 at 12:51
• What exactly is the question here? Are you simply asking how to speed up your C code? pow is not a fast function; that is where your optimization efforts need to be focused. But first, you need to actually ensure that the code runs too slowly, rather than simply presuming that it is inefficient. What is your desired speed range? How far off are you currently? – Cody Gray Jan 12 '17 at 13:00
• @LokiAstari My question is how to implement the explicit polynomial formual in C. For example, $f(x)=a_3x^3+a_2x^2+a_1x+a_0=x(x(a_3x+a_2)+a_1)+a_0$ by the horner algorithm. But for the long formual, do you think whether it is a good method to apply horner algorithm? – xyz Jan 12 '17 at 13:05
• @CodyGray I ask this question just for improving the C code, because I have no experience in C. – xyz Jan 12 '17 at 13:09

## 1 Answer

The compiler (when -fastmath or equivalent is active) will do a lot of optimizations for you, benchmark the current code and see if it's fast enough for your purposes. This means that the relative amount of time spent in that code is too small for speedups in that part to yield better results.

First step is to replace the pow calls with explicit (cached) multiplication. pow expects 2 floating point values and is a complex formula though it will most likely have a fast path for integer exponents (which after inlining will result in the same but I have little faith in compiler optimizations):

double calc_A(int i, const double *U) {
double u0, u1, u2, u3;
double res;
u0 = U[i];
u1 = U[i+1];
u2 = U[i+2];
u3 = U[i+3];

double u0_2 = u0*u0;
double u1_2 = u1*u1;
double u1_3 = u1_2*u1;
double u1_4 = u1_2*u1_2;
//...
res = 3 * u1_4 + 3 * u1_3*(u2 - 2 * u3) +
u1_2*(3 * u2_2 - 7 * u2*u3 + 3 * u3_2) +
u1*u2*(3 * u2_2 - 8 * u2*u3 + 4 * u3_2) +
//...
}


Looking at the use case (calling it for every index in an array) it would pay to SIMD the function if your cpu supports double precision vector math (modern cpus almost certainly do).

This will mean using intrinsics for the math which will be less clear code but will probably be faster.

• Thanks for your help:). I got your idea, that is, $x^2=x*x,\quad x^3=(x^2)*x \quad x^4=(x^2)*x$ – xyz Jan 12 '17 at 14:41