Polynomial of degree N

 

whit N from 2 to 20



Polynom N-order with N = 11
Pn3D_o11_006a.sca

Polynomial of degree N

whith dynamic calculation of the coefficients

The algorithm used to calculate the polynomial attractor of degree 2 (quadratic), from point Pn(x,y,z) to new point Pn+1(xNew, yNew, zNew) is given from formula:

    xNew = k[ 0] + k[ 1] * x + k[ 2] * y + k[ 3] * z + k[ 4] * x * x + k[ 5] * x * y + k[ 6] * x * z + k[ 7] * y * y + k[ 8] * y * z + k[ 9] * z * z;
    yNew = k[10] + k[11] * x + k[12] * y + k[13] * z + k[14] * x * x + k[15] * x * y + k[16] * x * z + k[17] * y * y + k[18] * y * z + k[19] * z * z;
    zNew = k[20] + k[21] * x + k[22] * y + k[23] * z + k[24] * x * x + k[25] * x * y + k[26] * x * z + k[27] * y * y + k[28] * y * z + k[29] * z * z;

where k[i] (with i form 0 to 29) are constants randomly calculated at the beginning, the variables (x,y,z) are the coordinates of the actual point Pn, and (xNew, yNew, zNew) will be the coordinates of new point Pn+1


for the 3rd order polynomial (cubic) it gets complicated, we pass from 30 to 60 coefficients:

    xNew= k[ 0] + k[ 1] * x + k[ 2] * y + k[ 3] * z + k[ 4] * x * x + k[ 5] * x * y + k[ 6] * x * z + k[ 7] * y * y + k[ 8] * y * z + k[ 9] * z * z + k[10] * x * x * x + k[11] * x * x * y + k[12] * x * x * z + k[13] * x * y * y + k[14] * x * y * z + k[15] * x * z * z + k[16] * y * y * y + k[17] * y * y * z + k[18] * y * z * z + k[19] * z * z * z
    yNew= k[20] + k[21] * x + k[22] * y + k[23] * z + k[24] * x * x + k[25] * x * y + k[26] * x * z + k[27] * y * y + k[28] * y * z + k[29] * z * z + k[30] * x * x * x + k[31] * x * x * y + k[32] * x * x * z + k[33] * x * y * y + k[34] * x * y * z + k[35] * x * z * z + k[36] * y * y * y + k[37] * y * y * z + k[38] * y * z * z + k[39] * z * z * z
    zNew= k[40] + k[41] * x + k[42] * y + k[43] * z + k[44] * x * x + k[45] * x * y + k[46] * x * z + k[47] * y * y + k[48] * y * z + k[49] * z * z + k[50] * x * x * x + k[51] * x * x * y + k[52] * x * x * z + k[53] * x * y * y + k[54] * x * y * z + k[55] * x * z * z + k[56] * y * y * y + k[57] * y * y * z + k[58] * y * z * z + k[59] * z * z * z

More in general the number of coefficients is given by the formula:
(order+1) * (order+2) * (order+3) / 2, where "order" is the degree of polynomial.
So for a 4rd order (quartic) we'll have 105 coefficients, and so on.

It would be nice to be able to calculate them dynamically...


We go back to the 2nd order polynomial (quadratic): if we use vectors k[i](x,y,x), we can rewrite it:

    xNew = k[0].x + k[1].x * x + k[2].x * y + k[3].x * z + k[4].x * x * x + k[5].x * x * y + k[6].x * x * z + k[7].x * y * y + k[8].x * y * z + k[9].x * z * z;
    yNew = k[0].y + k[1].y * x + k[2].y * y + k[3].y * z + k[4].y * x * x + k[5].y * x * y + k[6].y * x * z + k[7].y * y * y + k[8].y * y * z + k[9].y * z * z;
    zNew = k[0].z + k[1].z * x + k[2].z * y + k[3].z * z + k[4].z * x * x + k[5].z * x * y + k[6].z * x * z + k[7].z * y * y + k[8].z * y * z + k[9].z * z * z;

 

We can write more verbosely the formula:

    xNew = k[0].x + k[1].x * x^1 + k[2].x * y^1 + k[3].x * z^1 + k[4].x * x^2 + k[5].x * x^1 * y^1 + k[6].x * x^1 * z^1 + k[7].x * y^2 + k[8].x * y^1 * z^1 + k[9].x * z^2;
    yNew = k[0].y + k[1].y * x^1 + k[2].y * y^1 + k[3].y * z^1 + k[4].y * x^2 + k[5].y * x^1 * y^1 + k[6].y * x^1 * z^1 + k[7].y * y^2 + k[8].y * y^1 * z^1 + k[9].y * z^2;
    zNew = k[0].z + k[1].z * x^1 + k[2].z * y^1 + k[3].z * z^1 + k[4].z * x^2 + k[5].z * x^1 * y^1 + k[6].z * x^1 * z^1 + k[7].z * y^2 + k[8].z * y^1 * z^1 + k[9].z * z^2;


and more prolixly... remembering also that a^0 = 1 we write:

    xNew = k[0].x * x^0 * y^0 * z^0 + k[1].x * x^1 * y^0 * z^0 + k[2].x * x^0 * y^1 * z^0 + k[3].x * x^0 * y^0 * z^1 + k[4].x * x^2 * y^0 * z^0 + k[5].x * x^1 * y^1 * z^0 + k[6].x * x^1 * y^0 * z^1 + k[7].x * x^0 * y^2 *z^0 + k[8].x * x^0 * y^1 * z^1 + k[9].x x^0 * y^0 * z^2;
    yNew = k[0].y * x^0 * y^0 * z^0 + k[1].y * x^1 * y^0 * z^0 + k[2].y * x^0 * y^1 * z^0 + k[3].y * x^0 * y^0 * z^1 + k[4].y * x^2 * y^0 * z^0 + k[5].y * x^1 * y^1 * z^0 + k[6].y * x^1 * y^0 * z^1 + k[7].y * x^0 * y^2 *z^0 + k[8].y * x^0 * y^1 * z^1 + k[9].y x^0 * y^0 * z^2;
    zNew = k[0].z * x^0 * y^0 * z^0 + k[1].z * x^1 * y^0 * z^0 + k[2].z * x^0 * y^1 * z^0 + k[3].z * x^0 * y^0 * z^1 + k[4].z * x^2 * y^0 * z^0 + k[5].z * x^1 * y^1 * z^0 + k[6].z * x^1 * y^0 * z^1 + k[7].z * x^0 * y^2 *z^0 + k[8].z * x^0 * y^1 * z^1 + k[9].z x^0 * y^0 * z^2;

 

now, if we set:

    float cf[10];
    cf[0] = x^0 * y^0 * z^0;
    cf[1] = x^1 * y^0 * z^0;
    cf[2] = x^0 * y^1 * z^0;
    ...
    cf[9] = x^0 * y^0 * z^2;

 

finally we can calculate the new point Pn+1 as summation of K * cf:

    vec3 pNew = vec3(0.0, 0.0, 0.0);
    for(int i; i<10; i++)
        pNew += k[i] * cf[i];


we need only to calculate dynamically the values of cf[0], cf[1], cf[2] ... etc.

So, we consider the degree of the polynomial (this case 2) and we calculate, for any variable (x, y, z), any power with exponents 0, 1 and 2.

    const int order = 2;
    float ex[order+1], ey[order+1], ez[order+1];

    ex[0] = ey[0] = ez[0] = 1.f; // x[0]=1 as x^0
    for(int i=1; i<=order; i++) {
        ex[i] = ex[i-1] * x;
        ey[i] = ey[i-1] * y;
        ez[i] = ez[i-1] * z;
    }
    // x[1] = x[0] * x = 1 * x --> x^1
    // and
    // x[2] = x[0] * x [1] * x --> x^2


usng also here the vectors, we rewrite the previous loop, with p(x,y,z) actual point:

    const int order = 2;
    static vec3 elv[order+1];

    elv[0] = vec3(1.f);
    for(int i=1; i<=order; i++)
        elv[i] = elv[i-1] * p;


and now we calculate the cf coefficients, dynamically... making sure that the sum of the exponents does not exceed the degree of the polynomial:

    int j = 0;
    // I use reverse loop because test on 0 is faster
    for(int x=order-1; x>=0; x--)
        for(int y=order-1; y>=0; y--)
            for(int z=order-1; z>=0; z--)
                if(x+y+z <= order)
                    cf[j++] = elv[x].x * elv[y].y * elv[z].z;

    // and now x^2*y^0*z^0... y^2, z^2 separately
    cf[j++] = elv[order].x;
    cf[j++] = elv[order].y;
    cf[j++] = elv[order].z;


so we arrived at the summation seen above

        int nCoeff = j;
        vec3 pNew = vec3(0.0, 0.0, 0.0);

        for(int i=0; i < nCoeff; i++) pNew += k[i] * cf[i];

the total number of coefficients can be calculate also from formula, previously seen:
nCoeff = ((order+1) * (order+2) * (order+3) / 2) / 3; .. (because we are using vectors of 3 components)

 

Now we write the whole formula, valid for any "order" of polynomial, so we can calculate also an degree of 20!, with 1771*3 = 5313 coefficients.

    // input : p, current (x,y,z) point
    // output: pNew, next calculate point
    void PowerN3D(vec3 &p, vec3 &pNew)
    {
        const int order = N; // where N is any int number
        vec3 elv[order+1];

        elv[0] = vec3(1.f);
        for(int i=1; i<=order; i++)
            elv[i] = elv[i-1] * p;

        int j = 0;
        // reverse index, only because compare with 0 is faster
        for(int x=order-1; x>=0; x--)
            for(int y=order-1; y>=0; y--)
                for(int z=order-1; z>=0; z--)
                    if(x+y+z <= order)
                        cf[j++] = elv[x].x * elv[y].y * elv[z].z;

        cf[j++] = elv[order].x;
        cf[j++] = elv[order].y;
        cf[j++] = elv[order].z;

        int nCoeff = j;
        pNew = vec3(0.f);

        for(int i=0; i < nCoeff; i++) pNew += kVal[i] * cf[i];
    }

This is a demonstrative function: need to declare and alloc cf
The source code is slightly different: I used std::vector and internal iterators, instead of static array, but it's a little more cryptic:

    void PowerN3D::Step(vec3 &p, vec3 &pNew)
    {
        elv[0] = vec3(1.f);
        for(int i=1; i<=order; i++) elv[i] = elv[i-1] * p;

        auto itCf = cf.begin();
        for(int x=order-1; x>=0; x--)
            for(int y=order-1; y>=0; y--)
                for(int z=order-1; z>=0; z--)
                    if(x+y+z <= order) *itCf++ = elv[x].x * elv[y].y * elv[z].z;

        *itCf++ = elv[order].x;
        *itCf++ = elv[order].y;
        *itCf++ = elv[order].z;

        pNew = vec3(0.f);

        itCf = cf.begin();
        for(auto &it : kVal) pNew += it * *itCf++;
    }