/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                          %%
%%                                                   Wroc{\l}aw, June 2019  %%
%%                                                                          %%
%%                                                                          %%
%%                                                                          %%
%%              LINEAR-TIME GEOMETRIC ALGORITHMS FOR EVALUATING             %%
%%                             B\'{E}ZIER CURVES                            %% 
%%                              (test program)                              %%
%%                                                                          %%
%%                           GNU C Compiler 7.4.0                           %%
%%                                                                          %%
%%                       Pawe{\l} Wo\'{z}ny, Filip Chudy                    %%
%%                                                                          %%
%%                                                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

*/

#include<stdio.h>
#include<stdlib.h>
#include<time.h>

#define MAXITER 500
#define CURVES 10000

float w_x[1000];
float w_y[1000];
float w_z[1000];
float weights[1000];
float q_x[1000];
float q_y[1000];
float q_z[1000];
float weightspace[1000 * 2];

float float_rand( float min, float max )
{
    float scale = rand() / (float) RAND_MAX; /* [0, 1.0] */
    return min + scale * ( max - min );      /* [min, max] */
}

void f_p_casteljau_2d(unsigned int n, float t)
{
    float t1 = 1.0 - t;
    for(int i = 0; i <= n; i++)
    {
        q_x[i] = w_x[i];
        q_y[i] = w_y[i];
    }
    for(int k = 1; k <= n; k++)
    {
        for(int i = 0; i <= n-k; i++)
        {
            q_x[i] = t1 * q_x[i] + t * q_x[i + 1];
            q_y[i] = t1 * q_y[i] + t * q_y[i + 1];
        }
    }
    return;
}

void f_p_new_2d(unsigned int n, float t)
{
    float h = 1.0;
    float u = 1.0 - t;
    unsigned int n1 = n + 1;
    q_x[0] = w_x[0];
    q_y[0] = w_y[0];
    if(t <= 0.5)
    {
        u = t / u;
        for(int k = 1; k <= n; k++)
        {
            h = h * u * (n1 - k);
            h = h / (k + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
        }
    }
    else
    {
        u = u / t;
        for(int k = 1; k <= n; k++)
        {
            h = h * (n1 - k);
            h = h / (k * u + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
        }
    }
    return;
}

void f_p_casteljau_3d(unsigned int n, float t)
{
    float t1 = 1.0 - t;
    for(int i = 0; i <= n; i++)
    {
        q_x[i] = w_x[i];
        q_y[i] = w_y[i];
        q_z[i] = w_z[i];
    }
    for(int k = 1; k <= n; k++)
    {
        for(int i = 0; i <= n-k; i++)
        {
            q_x[i] = t1 * q_x[i] + t * q_x[i + 1];
            q_y[i] = t1 * q_y[i] + t * q_y[i + 1];
            q_z[i] = t1 * q_z[i] + t * q_z[i + 1];
        }
    }
    return;
}

void f_p_new_3d(unsigned int n, float t)
{
    float h = 1.0;
    float u = 1.0 - t;
    unsigned int n1 = n + 1;
    q_x[0] = w_x[0];
    q_y[0] = w_y[0];
    q_z[0] = w_z[0];
    if(t <= 0.5)
    {
        u = t / u;
        for(int k = 1; k <= n; k++)
        {
            h = h * u * (n1 - k);
            h = h / (k + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
            q_z[0] = h1 * q_z[0] + h * w_z[k];
        }
    }
    else
    {
        u = u / t;
        for(int k = 1; k <= n; k++)
        {
            h = h * (n1 - k);
            h = h / (k * u + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
            q_z[0] = h1 * q_z[0] + h * w_z[k];
        }
    }
    return;
}

void f_r_casteljau_2d(unsigned int n, float t)
{
    float t1 = 1.0 - t;
    for(int i = 0; i <= n; i++)
    {
        q_x[i] = w_x[i];
        q_y[i] = w_y[i];
        weightspace[i] = weights[i];
    }
    for(int k = 1; k <= n; k++)
    {
        for(int i = 0; i <= n-k; i++)
        {
            float u = t1 * weightspace[i]; 
            float v = t * weightspace[i + 1]; 
            weightspace[i] = u + v;
            u /= weightspace[i];
            v = 1.0 - u;
            q_x[i] = u * q_x[i] + v * q_x[i + 1];
            q_y[i] = u * q_y[i] + v * q_y[i + 1];
        }
    }
    return;
}

void f_r_new_2d(unsigned int n, float t)
{
    float h = 1.0;
    float u = 1.0 - t;
    unsigned int n1 = n + 1;
    q_x[0] = w_x[0];
    q_y[0] = w_y[0];
    if(t <= 0.5)
    {
        u = t / u;
        for(int k = 1; k <= n; k++)
        {
            h = h * u * (n1 - k) * weights[k];
            h = h / (k * weights[k - 1] + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
        }
    }
    else
    {
        u = u / t;
        for(int k = 1; k <= n; k++)
        {
            h = h * (n1 - k) * weights[k];
            h = h / (k * u * weights[k - 1] + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
        }
    }
    return;
}

void f_r_casteljau_3d(unsigned int n, float t)
{
    float t1 = 1.0 - t;
    for(int i = 0; i <= n; i++)
    {
        q_x[i] = w_x[i];
        q_y[i] = w_y[i];
        q_z[i] = w_z[i];
        weightspace[i] = weights[i];
    }
    for(int k = 1; k <= n; k++)
    {
        for(int i = 0; i <= n-k; i++)
        {
            float u = t1 * weightspace[i]; 
            float v = t * weightspace[i + 1]; 
            weightspace[i] = u + v;
            u /= weightspace[i];
            v = 1.0 - u;
            q_x[i] = u * q_x[i] + v * q_x[i + 1];
            q_y[i] = u * q_y[i] + v * q_y[i + 1];
            q_z[i] = u * q_z[i] + v * q_z[i + 1];
        }
    }
    return;
}

void f_r_new_3d(unsigned int n, float t)
{
    float h = 1.0;
    float u = 1.0 - t;
    unsigned int n1 = n + 1;
    q_x[0] = w_x[0];
    q_y[0] = w_y[0];
    q_z[0] = w_z[0];
    if(t <= 0.5)
    {
        u = t / u;
        for(int k = 1; k <= n; k++)
        {
            h = h * u * (n1 - k) * weights[k];
            h = h / (k * weights[k - 1] + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
            q_z[0] = h1 * q_z[0] + h * w_z[k];
        }
    }
    else
    {
        u = u / t;
        for(int k = 1; k <= n; k++)
        {
            h = h * (n1 - k) * weights[k];
            h = h / (k * u * weights[k - 1] + h);
            float h1 = 1 - h;
            q_x[0] = h1 * q_x[0] + h * w_x[k];
            q_y[0] = h1 * q_y[0] + h * w_y[k];
            q_z[0] = h1 * q_z[0] + h * w_z[k];
        }
    }
    return;
}

void main()
{
    clock_t start_t, end_t;
    for(int n = 1; n <= 20; n++)
    {
        if(n <= 6 || n % 5 == 0)
        {
            float t;
            double fp_old_2d = 0.0, fp_new_2d = 0.0;
            double fp_old_3d = 0.0, fp_new_3d = 0.0;
            double fr_old_2d = 0.0, fr_new_2d = 0.0;
            double fr_old_3d = 0.0, fr_new_3d = 0.0;
            //printf("==== n = %d ====\n", n);
            for(int v = 0; v < CURVES; v++)
            {
                for(int i = 0; i <= n; i++) w_x[i] = float_rand(-1.0, 1.0);
                for(int i = 0; i <= n; i++) w_y[i] = float_rand(-1.0, 1.0);
                for(int i = 0; i <= n; i++) w_z[i] = float_rand(-1.0, 1.0);
                for(int i = 0; i <= n; i++) weights[i] = float_rand(0.01, 1.0);
    
                for(int i = 0; i <= MAXITER; i++)
                {
                    t = (1. * i) / MAXITER;
                
                    start_t = clock();
                    f_p_casteljau_2d(n, t);
                    end_t = clock();
                    fp_old_2d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
    
                    start_t = clock();
                    f_p_new_2d(n, t);
                    end_t = clock();
                    fp_new_2d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            
                    start_t = clock();
                    f_p_casteljau_3d(n, t);
                    end_t = clock();
                    fp_old_3d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
    
                    start_t = clock();
                    f_p_new_3d(n, t);
                    end_t = clock();
                    fp_new_3d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            
                    start_t = clock();
                    f_r_casteljau_2d(n, t);
                    end_t = clock();
                    fr_old_2d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
    
                    start_t = clock();
                    f_r_new_2d(n, t);
                    end_t = clock();
                    fr_new_2d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            
                    start_t = clock();
                    f_r_casteljau_3d(n, t);
                    end_t = clock();
                    fr_old_3d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
        
                    start_t = clock();
                    f_r_new_3d(n, t);
                    end_t = clock();
                    fr_new_3d += (double)(end_t - start_t) / CLOCKS_PER_SEC;
                }
            }
            printf("%d & 2", n);

            if(fp_new_2d < fp_old_2d)
            {
                printf(" & \\textbf{%1.3lf} & %1.3lf", fp_new_2d, fp_old_2d);
            }
            else
            {
                printf(" & %1.3lf & \\textbf{%1.3lf}", fp_new_2d, fp_old_2d);
            }
            if(fr_new_2d < fr_old_2d)
            {
                printf(" & \\textbf{%1.3lf} & %1.3lf", fr_new_2d, fr_old_2d);
            }
            else
            {
                printf(" & %1.3lf & \\textbf{%1.3lf}", fr_new_2d, fr_old_2d);
            }
            printf("\\\\\n & 3");
            if(fp_new_3d < fp_old_3d)
            {
                printf(" & \\textbf{%1.3lf} & %1.3lf", fp_new_3d, fp_old_3d);
            }
            else
            {
                printf(" & %1.3lf & \\textbf{%1.3lf}", fp_new_3d, fp_old_3d);
            }
            if(fr_new_3d < fr_old_3d)
            {
                printf(" & \\textbf{%1.3lf} & %1.3lf", fr_new_3d, fr_old_3d);
            }
            else
            {
                printf(" & %1.3lf & \\textbf{%1.3lf}", fr_new_3d, fr_old_3d);
            }
            printf(" \\\\ \\hline\n");
        }
    }
    return;
}