A quartic polynomial in four variables, for numerical integration

Question

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[2]~A[6]
    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.

Answer 1

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.