//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%                                                                          %%
//%%                                                        Wroc{\l}aw, 2023  %%
//%%                                                                          %%
//%%                                                                          %%
//%%                                                                          %%
//%%                     FAST EVALUATION OF DERIVATIVES OF                    %%
//%%                             B\'{E}ZIER CURVES                            %%
//%%                               (test program)                             %%
//%%                                                                          %%
//%%                                 G++ 11.3.0                               %%
//%%                                                                          %%
//%%                                                                          %%
//%%                     Filip Chudy, Pawe{\l} Wo\'{z}ny                      %%
//%%                                                                          %%
//%%                                                                          %%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#include <cmath>
#include <ctime>
#include <iostream>
#include <random>

class InvalidDerivation : public std::exception
{
    public:
        const char* what()
        {
            return "InvalidDerivation";
        }
};

class InvalidDimensions : public std::exception
{
    public:
        const char* what()
        {
            return "InvalidDimensions";
        }
};

void evaluate_bernstein(int n, double t, double* result)
{
    double t1 = 1.0 - t;
    double tt;
    if(t < 0.5)
    {
        tt = t / t1;
        result[0] = pow(t1, n);
        for(int j = 1; j <= n; j++)
            result[j] = (result[j - 1] * tt * (n - j + 1)) / j;
    }
    else
    {
        tt = t1 / t;
        result[n] = pow(t, n);
        for(int j = n; j > 0; j--)
            result[j - 1] =  (result[j] * tt * j) / (n - j + 1);
    }
}

void fill_bernstein_derivatives(int n, double t, const double* source, double* result)
{
    result[0] = -n * source[0] - source[1];
    for(int j = 1; j < n; j++)
        result[j] = (n - j + 1) * source[j - 1] + (2 * j - n) * source[j] 
                  - (j + 1) * source[j + 1];
    result[n] = source[n - 1] + n * source[n];
}


class Point2D
{
    public:
        double x, y;
        Point2D(double x, double y) : x(x), y(y) {}
        Point2D() : x(0.0), y(0.0) {}

        Point2D operator + (Point2D other) const
        {
            return {x + other.x, y + other.y};
        }
        Point2D operator - (Point2D other) const
        {
            return {x - other.x, y - other.y};
        }

        Point2D operator * (double scalar) const
        {
            return {x * scalar, y * scalar};
        }

        Point2D& operator = (const Point2D& pt)
        = default;
};

Point2D genAlg(unsigned int n, const Point2D* points, const double* basis)
{
    double h = 1.0;
    Point2D Q = points[0];
    for(int k = 1; k <= n; k++)
    {
        h = basis[k] * h;
        h = h / (h + basis[k - 1]);
        Q = Q * (1.0 - h) + points[k] * h;
    }
    return Q;
}

Point2D genAlgPosNeg(unsigned int n, const double* weights, const Point2D* points, const double* basis)
{
    double hw, hw_pos_sum = 0.0, hw_neg_sum = 0.0, hw_pos[n + 1], hw_neg[n + 1];
    Point2D points_pos[n + 1], points_neg[n + 1];
    int pos_cnt = 0, neg_cnt = 0;
    for(int k = 0; k <= n; k++)
    {
        hw = weights[k] * basis[k];
        if(hw < 0)
        {
            hw = -hw;
            hw_neg_sum += hw;
            hw_neg[neg_cnt] = hw;
            points_neg[neg_cnt] = points[k];
            neg_cnt += 1;
        }
        else if(hw > 0)
        {
            hw_pos_sum += hw;
            hw_pos[pos_cnt] = hw;
            points_pos[pos_cnt] = points[k];
            pos_cnt += 1;
        }

    }
    for(int k = 0; k < neg_cnt; k++) hw_neg[k] /= hw_neg_sum;
    for(int k = 0; k < pos_cnt; k++) hw_pos[k] /= hw_pos_sum;

    return (genAlg(pos_cnt - 1, points_pos, hw_pos) * hw_pos_sum 
        - genAlg(neg_cnt - 1, points_neg, hw_neg) * hw_neg_sum) * (1.0 / (hw_pos_sum - hw_neg_sum));

}

Point2D geomVector(unsigned int n, const double* weights, const Point2D* points, const double* basis)
{
    Point2D result = Point2D(0.0, 0.0);
    double remain_coefficient = weights[0] * basis[0];
    for(int k = 0; k < n; k++)
    {
        Point2D diff = points[k] - points[k + 1];
        result = result + diff * remain_coefficient;
        remain_coefficient += weights[k + 1] * basis[k + 1];
    }
    return result;
}


class Point3D
{
    public:
        double x, y, z;
        Point3D(double x, double y, double z) : x(x), y(y), z(z) {}
        Point3D() : x(0.0), y(0.0), z(0.0) {}

        Point3D operator + (Point3D other) const
        {
            return {x + other.x, y + other.y, z + other.z};
        }
        Point3D operator - (Point3D other) const
        {
            return {x - other.x, y - other.y, z - other.z};
        }

        Point3D operator * (double scalar) const
        {
            return {x * scalar, y * scalar, z * scalar};
        }

        Point3D& operator = (const Point3D& pt)
        = default;
};

Point3D genAlg(unsigned int n, const Point3D* points, const double* basis)
{
    double h = 1.0;
    Point3D Q = points[0];
    for(int k = 1; k <= n; k++)
    {
        h = basis[k] * h;
        h = h / (h + basis[k - 1]);
        Q = Q * (1.0 - h) + points[k] * h;
    }
    return Q;
}

Point3D genAlgPosNeg(unsigned int n, const double* weights, const Point3D* points, const double* basis)
{
    double hw, hw_pos_sum = 0.0, hw_neg_sum = 0.0, hw_pos[n + 1], hw_neg[n + 1];
    Point3D points_pos[n + 1], points_neg[n + 1];
    int pos_cnt = 0, neg_cnt = 0;
    for(int k = 0; k <= n; k++)
    {
        hw = weights[k] * basis[k];
        if(hw < 0)
        {
            hw = -hw;
            hw_neg_sum += hw;
            hw_neg[neg_cnt] = hw;
            points_neg[neg_cnt] = points[k];
            neg_cnt += 1;
        }
        else if(hw > 0)
        {
            hw_pos_sum += hw;
            hw_pos[pos_cnt] = hw;
            points_pos[pos_cnt] = points[k];
            pos_cnt += 1;
        }

    }
    for(int k = 0; k < neg_cnt; k++) hw_neg[k] /= hw_neg_sum;
    for(int k = 0; k < pos_cnt; k++) hw_pos[k] /= hw_pos_sum;

    return (genAlg(pos_cnt - 1, points_pos, hw_pos) * hw_pos_sum 
        - genAlg(neg_cnt - 1, points_neg, hw_neg) * hw_neg_sum) * (1.0 / (hw_pos_sum - hw_neg_sum));

}

Point3D geomVector(unsigned int n, const double* weights, const Point3D* points, const double* basis)
{
    Point3D result = Point3D(0.0, 0.0, 0.0);
    double remain_coefficient = weights[0] * basis[0];
    for(int k = 0; k < n; k++)
    {
        Point3D diff = points[k] - points[k + 1];
        result = result + diff * remain_coefficient;
        remain_coefficient += weights[k + 1] * basis[k + 1];
    }
    return result;
}


class Curve1D
{
    public:
        int n;
        double* points;

        Curve1D(int n, const double* old_points) : n(n)
        {
            points = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = old_points[i];
            }
        }

        explicit Curve1D(int n) : n(n)
        {
            points = new double[n + 1];
        }

        Curve1D() : n(1), points(nullptr) {}

        ~Curve1D()
        {
            delete points;
            points = nullptr;
        }

        Curve1D& operator = (const Curve1D& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
            }
            return *this;
        }

        [[nodiscard]] double evaluate_de_casteljau(double t) const
        {
            double aux[n + 1];
            double t1 = 1.0 - t;
            for(int i = 0; i <= n; i++) aux[i] = points[i];
            for(int k = 1; k <= n; k++)
                for(int i = 0; i <= n - k; i++)
                    aux[i] = aux[i] * t1 + aux[i + 1] * t;
            return aux[0];
        }

        [[nodiscard]] double evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            double q = points[0];
            if( t <= 0.5)
            {
                u = t / u;
                for(int k = 1; k <= n; k++)
                {
                    h = h * u * (n1 - k);
                    h = h / (k + h);
                    double h1 = 1 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            else
            {
                u = u / t;
                for(int k = 1; k <= n; k++)
                {
                    h = h * (n1 - k);
                    h = h / (k * u + h);
                    double h1 = 1.0 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            return q;
        }

        double evaluate_wozny_chudy_h(double t, const double* h) const
        {
            double q = points[0];
            for(int k = 1; k <= n; k++)
            {
                q = q * (1 - h[k]) + points[k] * h[k];
            }
            return q;
        }

        [[nodiscard]] Curve1D derive_lower_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve1D res = Curve1D(n - 1);
            for(int i = 0; i <= n - 1; i++)
            {
                res.points[i] = (points[i + 1] - points[i]) * n;
            }
            return res;
        }

        [[nodiscard]] Curve1D derive_keep_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve1D res = Curve1D(n);
            for(int i = 0; i <= n - 1; i++)
            {
                int i1 = i + 1;
                double delta = points[i1] - points[i];
                res.points[i] = res.points[i] + delta * (n - i);
                res.points[i1] = delta * i1;
            }
            return res;
        }

        void derive_lower_degree_inplace()
        {
            if (n == 0) throw InvalidDerivation();
            for(int i = 0; i <= n - 1; i++) points[i] = (points[i + 1] - points[i]) * n;
            n--;
        }
};

class Curve1DVol
{

public:
    int n;
    double* points;

    Curve1DVol(int n, double* old_points) : n(n), points(old_points) {}

    Curve1DVol& operator = (const Curve1DVol& c)
    = default;

    Curve1DVol& operator = (const Curve1D& c)
    {
        n = c.n;
        points = c.points;
        return *this;
    }

    [[nodiscard]] double evaluate_wozny_chudy(double t) const
    {
        double h = 1.0;
        double u = 1.0 - t;
        unsigned int n1 = n + 1;
        double q = points[0];
        if( t <= 0.5)
        {
            u = t / u;
            for(int k = 1; k <= n; k++)
            {
                h = h * u * (n1 - k);
                h = h / (k + h);
                double h1 = 1 - h;
                q = q * h1 + points[k] * h;
            }
        }
        else
        {
            u = u / t;
            for(int k = 1; k <= n; k++)
            {
                h = h * (n1 - k);
                h = h / (k * u + h);
                double h1 = 1.0 - h;
                q = q * h1 + points[k] * h;
            }
        }
        return q;
    }

    double evaluate_wozny_chudy_h(double t, const double* h) const
    {
        double q = points[0];
        for(int k = 1; k <= n; k++)
        {
            q = q * (1 - h[k]) + points[k] * h[k];
        }
        return q;
    }
};

class Curve2D
{
    public:
        int n;
        Point2D* points;

        Curve2D(int n, Point2D* old_points) : n(n)
        {
            points = new Point2D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = old_points[i];
            }
        }

        explicit Curve2D(int n) : n(n)
        {
            points = new Point2D[n + 1];
        }

        Curve2D() : n(1), points(nullptr) {}

        ~Curve2D()
        {
            delete points;
            points = nullptr;
        }

        Curve2D& operator = (const Curve2D& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point2D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
            }
            return *this;
        }

        [[nodiscard]] Point2D evaluate_de_casteljau(double t) const
        {
            Point2D aux[n + 1];
            double t1 = 1.0 - t;
            for(int i = 0; i <= n; i++) aux[i] = points[i];
            for(int k = 1; k <= n; k++)
                for(int i = 0; i <= n - k; i++)
                    aux[i] = aux[i] * t1 + aux[i + 1] * t;
            return aux[0];
        }

        [[nodiscard]] Point2D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point2D q = points[0];
            if( t <= 0.5)
            {
                u = t / u;
                for(int k = 1; k <= n; k++)
                {
                    h = h * u * (n1 - k);
                    h = h / (k + h);
                    double h1 = 1 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            else
            {
                u = u / t;
                for(int k = 1; k <= n; k++)
                {
                    h = h * (n1 - k);
                    h = h / (k * u + h);
                    double h1 = 1.0 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            return q;
        }

        Point2D evaluate_wozny_chudy_h(double t, double* h) const
        {
            Point2D q = points[0];
            for(int k = 1; k <= n; k++)
            {
                q = q * (1 - h[k]) + points[k] * h[k];
            }
            return q;
        }

        [[nodiscard]] Curve2D derive_lower_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve2D res = Curve2D(n - 1);
            for(int i = 0; i <= n - 1; i++)
            {
                res.points[i] = (points[i + 1] - points[i]) * n;
            }
            return res;
        }

        [[nodiscard]] Curve2D derive_keep_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve2D res = Curve2D(n);
            for(int i = 0; i <= n - 1; i++)
            {
                int i1 = i + 1;
                Point2D delta = points[i1] - points[i];
                res.points[i] = res.points[i] + delta * (n - i);
                res.points[i1] = delta * i1;
            }
            return res;
        }
};

class Curve2DVol
{
    public:
        int n;
        Point2D* points;

        Curve2DVol(int n, Point2D* old_points) : n(n), points(old_points) {}

        ~Curve2DVol()
        {
            points = nullptr;
        }

        Curve2DVol& operator = (const Curve2DVol& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point2D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
            }
            return *this;
        }

        [[nodiscard]] Point2D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point2D q = points[0];
            if (t <= 0.5) {
                u = t / u;
                for (int k = 1; k <= n; k++) {
                    h = h * u * (n1 - k);
                    h = h / (k + h);
                    double h1 = 1 - h;
                    q = q * h1 + points[k] * h;
                }
            } else {
                u = u / t;
                for (int k = 1; k <= n; k++) {
                    h = h * (n1 - k);
                    h = h / (k * u + h);
                    double h1 = 1.0 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            return q;
        }
        Point2D evaluate_wozny_chudy_h(double t, double* h) const
        {
            Point2D q = points[0];
            for(int k = 1; k <= n; k++)
            {
                q = q * (1 - h[k]) + points[k] * h[k];
            }
            return q;
        }
};

class RatCurve2DVol
{
    public:
        int n;
        Point2D* points;
        double* weights;

        RatCurve2DVol(int n, Point2D* old_points, double* old_weights) : n(n), points(old_points), weights(old_weights) {}

        RatCurve2DVol& operator = (const RatCurve2DVol& c)
        = default;

        [[nodiscard]] Point2D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            int n1 = n + 1;
            Point2D Q = points[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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            return Q;
        }
};

class RatCurve2D
{
    public:
        int n;
        Point2D* points;
        double* weights;

        RatCurve2D(int n, const Point2D* old_points,
                   const double* old_weights) : n(n)
        {
            points = new Point2D[n + 1];
            weights = new double[n + 1];

            for(int i = 0; i <= n; i++)
            {
                points[i] = old_points[i];
                weights[i] = old_weights[i];
            }
        }

        RatCurve2D() : n(1), points(nullptr), weights(nullptr) {}

        ~RatCurve2D()
        {
            delete points;
            points = nullptr;
            delete weights;
            weights = nullptr;
        }

        RatCurve2D& operator = (const Curve2D& c)
        {
            delete points;
            n = c.n;
            points = new Point2D[n + 1];
            weights = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
                weights[i] = 1.0;
            }
            return *this;
        }

        RatCurve2D& operator = (const RatCurve2D& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point2D[n + 1];
            weights = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
                weights[i] = c.weights[i];
            }
            return *this;
        }

        [[nodiscard]] Point2D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            int n1 = n + 1;
            Point2D Q = points[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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            return Q;
        }

        void evaluate_wozny_chudy_inter(double t, double* hs, Point2D* Qs) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point2D Q = points[0];
            hs[0] = h;
            Qs[0] = Q;
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                    hs[k] = h;
                    Qs[k] = Q;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                    hs[k] = h;
                    Qs[k] = Q;
                }
            }
        }

        void derive_simple_floater(double t, Point2D* derivatives,
                                   bool skip_second= true) const
        {
            Point2D Q[3][3];
            double ww[3][3];
            unsigned int r;
            Point2D aux[n + 1];
            double w[n + 1];
            double t1 = 1.0 - t;
            double u, v;
            for(int i = 0; i <= n; i++)
            {
                aux[i] = points[i];
                w[i] = weights[i];
            }
            if(n <= 2)
            {
                int j_max = 2 - n;
                for(int j = 0; j <= 2 - j_max; j++)
                {
                    Q[j_max][j] = aux[j];
                    ww[j_max][j] = w[j];
                }
            }
            for(int k = 1; k <= n; k++)
            {
                for(int i = 0; i <= n - k; i++)
                {
                    u = t1 * w[i];
                    v = t * w[i + 1];
                    w[i] = u + v;
                    u /= w[i];
                    v = 1.0 - u;
                    aux[i] = aux[i] * u + aux[i + 1] * v;
                }
                if(k >= n - 2)
                {
                    int j_max = k - n + 2;
                    for(int j = 0; j <= 2 - j_max; j++)
                    {
                        Q[j_max][j] = aux[j];
                        ww[j_max][j] = w[j];
                    }
                }
            }
            derivatives[0] = Q[2][0];
            derivatives[1] = (Q[1][1] - Q[1][0])
                             * (n * ww[1][0] * ww[1][1] / (ww[2][0] * ww[2][0]));
            if(!skip_second)
            {
                double ww203 = ww[2][0] * ww[2][0] * ww[2][0];
                double w1 = n * ww[0][2] / ww203
                            * (2 * n * ww[1][0] * ww[1][0] - (n - 1) * ww[0][0] * ww[2][0] - 2 * ww[1][0] * ww[2][0]);
                double w2 = n * ww[0][0] / ww203
                            * (2 * n * ww[1][1] * ww[1][1] - (n - 1) * ww[0][2] * ww[2][0] - 2 * ww[1][1] * ww[2][0]);
                derivatives[2] = (Q[0][2] - Q[0][1]) * w1 - (Q[0][1] - Q[0][0]) * w2;
            }
        }

        void derive_new_floater(double t, Point2D* derivatives,
                                bool skip_second= true) const
        {
            Point2D Q[3][3];
            double ww[3][3];
            double u, v, t1 = 1.0 - t;

            int r = 1;
            if(skip_second)
            {
                r = 2;
                RatCurve2DVol Qc0 = RatCurve2DVol(n - 1, points, weights);
                RatCurve2DVol Qc1 = RatCurve2DVol(n - 1, points + 1, weights + 1);
                Q[1][0] = Qc0.evaluate_wozny_chudy(t);
                Q[1][1] = Qc1.evaluate_wozny_chudy(t);
                Curve1DVol wc0 = Curve1DVol(n - 1, weights);
                Curve1DVol wc1 = Curve1DVol(n - 1, weights + 1);
                ww[1][0] = wc0.evaluate_wozny_chudy(t);
                ww[1][1] = wc1.evaluate_wozny_chudy(t);
            }
            else
            {
                RatCurve2DVol Qc0 = RatCurve2DVol(n - 2, points, weights);
                RatCurve2DVol Qc1 = RatCurve2DVol(n - 2, points + 1, weights + 1);
                RatCurve2DVol Qc2 = RatCurve2DVol(n - 2, points + 2, weights + 2);
                Q[0][0] = Qc0.evaluate_wozny_chudy(t);
                Q[0][1] = Qc1.evaluate_wozny_chudy(t);
                Q[0][2] = Qc2.evaluate_wozny_chudy(t);
                Curve1DVol wc0 = Curve1DVol(n - 2, weights);
                Curve1DVol wc1 = Curve1DVol(n - 2, weights + 1);
                Curve1DVol wc2 = Curve1DVol(n - 2, weights + 2);
                ww[0][0] = wc0.evaluate_wozny_chudy(t);
                ww[0][1] = wc1.evaluate_wozny_chudy(t);
                ww[0][2] = wc2.evaluate_wozny_chudy(t);
            }
            for(int k = r; k <= 2; k++)
            {
                for(int i = 0; i <= 2 - k; i++)
                {
                    u = t1 * ww[k - 1][i];
                    v = t * ww[k - 1][i + 1];
                    ww[k][i] = u + v;
                    u /= ww[k][i];
                    v = 1.0 - u;
                    Q[k][i] = Q[k - 1][i] * u + Q[k - 1][i + 1] * v;
                }
            }

            derivatives[0] = Q[2][0];

            derivatives[1] = (Q[1][1] - Q[1][0])
                             * (n * ww[1][0] * ww[1][1] / (ww[2][0] * ww[2][0]));
            if(!skip_second)
            {
                double ww203 = ww[2][0] * ww[2][0] * ww[2][0];
                double w1 = n * ww[0][2] / ww203
                            * (2 * n * ww[1][0] * ww[1][0] - (n - 1) * ww[0][0] * ww[2][0] - 2 * ww[1][0] * ww[2][0]);
                double w2 = n * ww[0][0] / ww203
                            * (2 * n * ww[1][1] * ww[1][1] - (n - 1) * ww[0][2] * ww[2][0] - 2 * ww[1][1] * ww[2][0]);
                derivatives[2] = (Q[0][2] - Q[0][1]) * w1 - (Q[0][1] - Q[0][0]) * w2;
            }
        }

        void derive_new_geom(double t, Point2D* derivatives, unsigned int r) const
        {
            Point2D Dk;
            double B[r + 1][n + 1], As[r + 1];
            int bin[r + 1][r + 1];
            evaluate_bernstein(n, t, B[0]);
            bin[0][0] = 1;
            for(int j = 1; j <= r; j++)
            {
                bin[j][0] = 1;
                for(int i = 1; i < j; i++)
                    bin[j][i] = bin[j - 1][i - 1] + bin[j - 1][i];
                bin[j][j] = 1;
            }

            derivatives[0] = evaluate_wozny_chudy(t);

            Curve1D w_curve = Curve1D(n, weights);
            As[0] = w_curve.evaluate_wozny_chudy(t);

            for(int k = 1; k <= r; k++)
            {
                fill_bernstein_derivatives(n, t, B[k - 1], B[k]); // OK

                w_curve.derive_lower_degree_inplace();

                As[k] = w_curve.evaluate_wozny_chudy(t);
                if(As[k] == 0.0)
                {
                    Dk = geomVector(n, weights, points, B[k]);
                    derivatives[k] = Dk;
                }
                else
                {
                    Dk = genAlgPosNeg(n, weights, points, B[k]);
                    derivatives[k] = (Dk - derivatives[0]) * As[k];
                }
                for(int i = 1; i <= k - 1; i++)
                    derivatives[k] = derivatives[k] - derivatives[i] * (bin[k][i] * As[k - i]);
                derivatives[k] = derivatives[k] * (1.0 / As[0]);
            }
        }

        void derive_new_h(double t, Point2D* derivatives, unsigned int r) const
        {
            if(t > 0.5)
            {
                RatCurve2D tmp_curve = RatCurve2D(n, points, weights);
                for(int k = 0; k <= n; k++)
                {
                    tmp_curve.points[k] = points[n - k];
                    tmp_curve.weights[k] = weights[n - k];
                }
                tmp_curve.derive_new_h(1.0 - t, derivatives, r);
            }
            else
            {
                unsigned int bin[r + 1][r + 1];
                bin[0][0] = 1;
                for(int j = 1; j <= r; j++)
                {
                    bin[j][0] = 1;
                    for(int i = 1; i < j; i++)
                        bin[j][i] = bin[j - 1][i - 1] + bin[j - 1][i];
                    bin[j][j] = 1;
                }
                Point2D Qs[r + 1][n + 1];
                double hs[r + 1][n + 1];
                evaluate_wozny_chudy_inter(t, hs[0], Qs[0]);
                derivatives[0] = Qs[0][n];

                for(int k = 1; k <= r; k++)
                {
                    hs[k][0] = 0.0;
                    Qs[k][0] = Point2D(0.0, 0.0);
                    for(int i = 1; i <= n; i++)
                    {

                        hs[k][i] = t * hs[k][i - 1] + k * hs[k - 1][i - 1];
                        hs[k][i] += weights[i - 1] * i * k * hs[k - 1][i]
                                    / (weights[i] * (n - i + 1));
                        for(int j = 1; j <= k; j++)
                        {
                            hs[k][i] -= bin[k][j] * hs[k - j][i]
                                        * (t * hs[j][i - 1] + j * hs[j - 1][i - 1]);
                        }
                        if(t == 0.0)
                        {
                            hs[k][i] *= (weights[i] * (n - i + 1)) / (weights[i - 1] * i);
                        }
                        else hs[k][i] *= hs[0][i] / (hs[0][i - 1] * t);
                        Qs[k][i] = (points[i] - Qs[0][i - 1]) * hs[k][i] + Qs[k][i - 1];
                        for(int j = 0; j <= k - 1; j++)
                        {
                            Qs[k][i] = Qs[k][i] - Qs[k - j][i - 1] * (bin[k][j] * hs[j][i]);
                        }
                    }
                    derivatives[k] = Qs[k][n];
                }
            }
        }
};


class Curve3D
{
    public:
        int n;
        Point3D* points;

        Curve3D(int n, Point3D* old_points) : n(n)
        {
            points = new Point3D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = old_points[i];
            }
        }

        explicit Curve3D(int n) : n(n)
        {
            points = new Point3D[n + 1];
        }

        Curve3D() : n(1), points(nullptr) {}

        ~Curve3D()
        {
            delete points;
            points = nullptr;
        }

        Curve3D& operator = (const Curve3D& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point3D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
            }
            return *this;
        }

        [[nodiscard]] Point3D evaluate_de_casteljau(double t) const
        {
            Point3D aux[n + 1];
            double t1 = 1.0 - t;
            for(int i = 0; i <= n; i++) aux[i] = points[i];
            for(int k = 1; k <= n; k++)
                for(int i = 0; i <= n - k; i++)
                    aux[i] = aux[i] * t1 + aux[i + 1] * t;
            return aux[0];
        }

        [[nodiscard]] Point3D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point3D q = points[0];
            if( t <= 0.5)
            {
                u = t / u;
                for(int k = 1; k <= n; k++)
                {
                    h = h * u * (n1 - k);
                    h = h / (k + h);
                    double h1 = 1 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            else
            {
                u = u / t;
                for(int k = 1; k <= n; k++)
                {
                    h = h * (n1 - k);
                    h = h / (k * u + h);
                    double h1 = 1.0 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            return q;
        }

        Point3D evaluate_wozny_chudy_h(double t, double* h) const
        {
            Point3D q = points[0];
            for(int k = 1; k <= n; k++)
            {
                q = q * (1 - h[k]) + points[k] * h[k];
            }
            return q;
        }

        [[nodiscard]] Curve3D derive_lower_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve3D res = Curve3D(n - 1);
            for(int i = 0; i <= n - 1; i++)
            {
                res.points[i] = (points[i + 1] - points[i]) * n;
            }
            return res;
        }

        [[nodiscard]] Curve3D derive_keep_degree() const
        {
            if (n == 0) throw InvalidDerivation();
            Curve3D res = Curve3D(n);
            for(int i = 0; i <= n - 1; i++)
            {
                int i1 = i + 1;
                Point3D delta = points[i1] - points[i];
                res.points[i] = res.points[i] + delta * (n - i);
                res.points[i1] = delta * i1;
            }
            return res;
        }
};

class Curve3DVol
{
    public:
        int n;
        Point3D* points;

        Curve3DVol(int n, Point3D* old_points) : n(n), points(old_points) {}

        ~Curve3DVol()
        {
            points = nullptr;
        }

        Curve3DVol& operator = (const Curve3DVol& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point3D[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
            }
            return *this;
        }

        [[nodiscard]] Point3D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point3D q = points[0];
            if (t <= 0.5) {
                u = t / u;
                for (int k = 1; k <= n; k++) {
                    h = h * u * (n1 - k);
                    h = h / (k + h);
                    double h1 = 1 - h;
                    q = q * h1 + points[k] * h;
                }
            } else {
                u = u / t;
                for (int k = 1; k <= n; k++) {
                    h = h * (n1 - k);
                    h = h / (k * u + h);
                    double h1 = 1.0 - h;
                    q = q * h1 + points[k] * h;
                }
            }
            return q;
        }
        Point3D evaluate_wozny_chudy_h(double t, double* h) const
        {
            Point3D q = points[0];
            for(int k = 1; k <= n; k++)
            {
                q = q * (1 - h[k]) + points[k] * h[k];
            }
            return q;
        }
};

class RatCurve3DVol
{
    public:
        int n;
        Point3D* points;
        double* weights;

        RatCurve3DVol(int n, Point3D* old_points, double* old_weights) : n(n), points(old_points), weights(old_weights) {}

        RatCurve3DVol& operator = (const RatCurve3DVol& c)
        = default;

        [[nodiscard]] Point3D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            int n1 = n + 1;
            Point3D Q = points[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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            return Q;
        }
};

class RatCurve3D
{
    public:
        int n;
        Point3D* points;
        double* weights;

        RatCurve3D(int n, const Point3D* old_points,
                   const double* old_weights) : n(n)
        {
            points = new Point3D[n + 1];
            weights = new double[n + 1];

            for(int i = 0; i <= n; i++)
            {
                points[i] = old_points[i];
                weights[i] = old_weights[i];
            }
        }

        RatCurve3D() : n(1), points(nullptr), weights(nullptr) {}

        ~RatCurve3D()
        {
            delete points;
            points = nullptr;
            delete weights;
            weights = nullptr;
        }

        RatCurve3D& operator = (const Curve3D& c)
        {
            delete points;
            n = c.n;
            points = new Point3D[n + 1];
            weights = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
                weights[i] = 1.0;
            }
            return *this;
        }

        RatCurve3D& operator = (const RatCurve3D& c)
        {
            if(this == &c) return *this;
            delete points;
            n = c.n;
            points = new Point3D[n + 1];
            weights = new double[n + 1];
            for(int i = 0; i <= n; i++)
            {
                points[i] = c.points[i];
                weights[i] = c.weights[i];
            }
            return *this;
        }

        [[nodiscard]] Point3D evaluate_wozny_chudy(double t) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            int n1 = n + 1;
            Point3D Q = points[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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                }
            }
            return Q;
        }

        void evaluate_wozny_chudy_inter(double t, double* hs, Point3D* Qs) const
        {
            double h = 1.0;
            double h1;
            double u = 1.0 - t;
            unsigned int n1 = n + 1;
            Point3D Q = points[0];
            hs[0] = h;
            Qs[0] = Q;
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                    hs[k] = h;
                    Qs[k] = Q;
                }
            }
            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);
                    h1 = 1.0 - h;
                    Q = Q * h1 + points[k] * h;
                    hs[k] = h;
                    Qs[k] = Q;
                }
            }
        }

        void derive_simple_floater(double t, Point3D* derivatives,
                                   bool skip_second= true) const
        {
            Point3D Q[3][3];
            double ww[3][3];
            unsigned int r;
            Point3D aux[n + 1];
            double w[n + 1];
            double t1 = 1.0 - t;
            double u, v;
            for(int i = 0; i <= n; i++)
            {
                aux[i] = points[i];
                w[i] = weights[i];
            }
            if(n <= 2)
            {
                int j_max = 2 - n;
                for(int j = 0; j <= 2 - j_max; j++)
                {
                    Q[j_max][j] = aux[j];
                    ww[j_max][j] = w[j];
                }
            }
            for(int k = 1; k <= n; k++)
            {
                for(int i = 0; i <= n - k; i++)
                {
                    u = t1 * w[i];
                    v = t * w[i + 1];
                    w[i] = u + v;
                    u /= w[i];
                    v = 1.0 - u;
                    aux[i] = aux[i] * u + aux[i + 1] * v;
                }
                if(k >= n - 2)
                {
                    int j_max = k - n + 2;
                    for(int j = 0; j <= 2 - j_max; j++)
                    {
                        Q[j_max][j] = aux[j];
                        ww[j_max][j] = w[j];
                    }
                }
            }
            derivatives[0] = Q[2][0];
            derivatives[1] = (Q[1][1] - Q[1][0])
                             * (n * ww[1][0] * ww[1][1] / (ww[2][0] * ww[2][0]));
            if(!skip_second)
            {
                double ww203 = ww[2][0] * ww[2][0] * ww[2][0];
                double w1 = n * ww[0][2] / ww203
                            * (2 * n * ww[1][0] * ww[1][0] - (n - 1) * ww[0][0] * ww[2][0] - 2 * ww[1][0] * ww[2][0]);
                double w2 = n * ww[0][0] / ww203
                            * (2 * n * ww[1][1] * ww[1][1] - (n - 1) * ww[0][2] * ww[2][0] - 2 * ww[1][1] * ww[2][0]);
                derivatives[2] = (Q[0][2] - Q[0][1]) * w1 - (Q[0][1] - Q[0][0]) * w2;
            }
        }

        void derive_new_floater(double t, Point3D* derivatives,
                                bool skip_second= true) const
        {
            Point3D Q[3][3];
            double ww[3][3];
            double u, v, t1 = 1.0 - t;

            int r = 1;
            if(skip_second)
            {
                r = 2;
                RatCurve3DVol Qc0 = RatCurve3DVol(n - 1, points, weights);
                RatCurve3DVol Qc1 = RatCurve3DVol(n - 1, points + 1, weights + 1);
                Q[1][0] = Qc0.evaluate_wozny_chudy(t);
                Q[1][1] = Qc1.evaluate_wozny_chudy(t);
                Curve1DVol wc0 = Curve1DVol(n - 1, weights);
                Curve1DVol wc1 = Curve1DVol(n - 1, weights + 1);
                ww[1][0] = wc0.evaluate_wozny_chudy(t);
                ww[1][1] = wc1.evaluate_wozny_chudy(t);
            }
            else
            {
                RatCurve3DVol Qc0 = RatCurve3DVol(n - 2, points, weights);
                RatCurve3DVol Qc1 = RatCurve3DVol(n - 2, points + 1, weights + 1);
                RatCurve3DVol Qc2 = RatCurve3DVol(n - 2, points + 2, weights + 2);
                Q[0][0] = Qc0.evaluate_wozny_chudy(t);
                Q[0][1] = Qc1.evaluate_wozny_chudy(t);
                Q[0][2] = Qc2.evaluate_wozny_chudy(t);
                Curve1DVol wc0 = Curve1DVol(n - 2, weights);
                Curve1DVol wc1 = Curve1DVol(n - 2, weights + 1);
                Curve1DVol wc2 = Curve1DVol(n - 2, weights + 2);
                ww[0][0] = wc0.evaluate_wozny_chudy(t);
                ww[0][1] = wc1.evaluate_wozny_chudy(t);
                ww[0][2] = wc2.evaluate_wozny_chudy(t);
            }
            for(int k = r; k <= 2; k++)
            {
                for(int i = 0; i <= 2 - k; i++)
                {
                    u = t1 * ww[k - 1][i];
                    v = t * ww[k - 1][i + 1];
                    ww[k][i] = u + v;
                    u /= ww[k][i];
                    v = 1.0 - u;
                    Q[k][i] = Q[k - 1][i] * u + Q[k - 1][i + 1] * v;
                }
            }

            derivatives[0] = Q[2][0];

            derivatives[1] = (Q[1][1] - Q[1][0])
                             * (n * ww[1][0] * ww[1][1] / (ww[2][0] * ww[2][0]));
            if(!skip_second)
            {
                double ww203 = ww[2][0] * ww[2][0] * ww[2][0];
                double w1 = n * ww[0][2] / ww203
                            * (2 * n * ww[1][0] * ww[1][0] - (n - 1) * ww[0][0] * ww[2][0] - 2 * ww[1][0] * ww[2][0]);
                double w2 = n * ww[0][0] / ww203
                            * (2 * n * ww[1][1] * ww[1][1] - (n - 1) * ww[0][2] * ww[2][0] - 2 * ww[1][1] * ww[2][0]);
                derivatives[2] = (Q[0][2] - Q[0][1]) * w1 - (Q[0][1] - Q[0][0]) * w2;
            }
        }

        void derive_new_geom(double t, Point3D* derivatives, unsigned int r) const
        {
            Point3D Dk;
            double B[r + 1][n + 1], As[r + 1];
            int bin[r + 1][r + 1];
            evaluate_bernstein(n, t, B[0]);
            bin[0][0] = 1;
            for(int j = 1; j <= r; j++)
            {
                bin[j][0] = 1;
                for(int i = 1; i < j; i++)
                    bin[j][i] = bin[j - 1][i - 1] + bin[j - 1][i];
                bin[j][j] = 1;
            }

            derivatives[0] = evaluate_wozny_chudy(t);

            Curve1D w_curve = Curve1D(n, weights);
            As[0] = w_curve.evaluate_wozny_chudy(t);

            for(int k = 1; k <= r; k++)
            {
                fill_bernstein_derivatives(n, t, B[k - 1], B[k]); // OK

                w_curve.derive_lower_degree_inplace();

                As[k] = w_curve.evaluate_wozny_chudy(t);
                if(As[k] == 0.0)
                {
                    Dk = geomVector(n, weights, points, B[k]);
                    derivatives[k] = Dk;
                }
                else
                {
                    Dk = genAlgPosNeg(n, weights, points, B[k]);
                    derivatives[k] = (Dk - derivatives[0]) * As[k];
                }
                for(int i = 1; i <= k - 1; i++)
                    derivatives[k] = derivatives[k] - derivatives[i] * (bin[k][i] * As[k - i]);
                derivatives[k] = derivatives[k] * (1.0 / As[0]);
            }
        }

        void derive_new_h(double t, Point3D* derivatives, unsigned int r) const
        {
            if(t > 0.5)
            {
                RatCurve3D tmp_curve = RatCurve3D(n, points, weights);
                for(int k = 0; k <= n; k++)
                {
                    tmp_curve.points[k] = points[n - k];
                    tmp_curve.weights[k] = weights[n - k];
                }
                tmp_curve.derive_new_h(1.0 - t, derivatives, r);
            }
            else
            {
                unsigned int bin[r + 1][r + 1];
                bin[0][0] = 1;
                for(int j = 1; j <= r; j++)
                {
                    bin[j][0] = 1;
                    for(int i = 1; i < j; i++)
                        bin[j][i] = bin[j - 1][i - 1] + bin[j - 1][i];
                    bin[j][j] = 1;
                }
                Point3D Qs[r + 1][n + 1];
                double hs[r + 1][n + 1];
                evaluate_wozny_chudy_inter(t, hs[0], Qs[0]);
                derivatives[0] = Qs[0][n];

                for(int k = 1; k <= r; k++)
                {
                    hs[k][0] = 0.0;
                    Qs[k][0] = Point3D(0.0, 0.0, 0.0);
                    for(int i = 1; i <= n; i++)
                    {

                        hs[k][i] = t * hs[k][i - 1] + k * hs[k - 1][i - 1];
                        hs[k][i] += weights[i - 1] * i * k * hs[k - 1][i]
                                    / (weights[i] * (n - i + 1));
                        for(int j = 1; j <= k; j++)
                        {
                            hs[k][i] -= bin[k][j] * hs[k - j][i]
                                        * (t * hs[j][i - 1] + j * hs[j - 1][i - 1]);
                        }
                        if(t == 0.0)
                        {
                            hs[k][i] *= (weights[i] * (n - i + 1)) / (weights[i - 1] * i);
                        }
                        else hs[k][i] *= hs[0][i] / (hs[0][i - 1] * t);
                        Qs[k][i] = (points[i] - Qs[0][i - 1]) * hs[k][i] + Qs[k][i - 1];
                        for(int j = 0; j <= k - 1; j++)
                        {
                            Qs[k][i] = Qs[k][i] - Qs[k - j][i - 1] * (bin[k][j] * hs[j][i]);
                        }
                    }
                    derivatives[k] = Qs[k][n];
                }
            }
        }
};


void evaluate_h(unsigned int n, double t, double* h)
{
    h[0] = 1.0;
    double u = 1.0 - t;
    unsigned int n1 = n + 1;
    if(t <= 0.5)
    {
        u = t / u;
        for(int k = 1; k <= n; k++)
        {
            h[k] = h[k - 1] * u * (n1 - k);
            h[k] = h[k] / (k + h[k]);
        }
    }
    else
    {
        u = u / t;
        for(int k = 1; k <= n; k++)
        {
            h[k] = h[k - 1] * (n1 - k);
            h[k] = h[k] / (k * u + h[k]);
        }
    }
}

void evaluate_derivatives_casteljau(Curve2D& c, unsigned int r, double t,
                                    Point2D* results)
{
    double t1 = 1.0 - t;
    Point2D Q[c.n + 1][c.n + 1];
    for (int k = 0; k <= c.n; k++) Q[0][k] = c.points[k];
    for (int j = 1; j <= c.n; j++)
        for (int k = 0; k <= c.n - j; k++)
            Q[j][k] = Q[j - 1][k] * t1 + Q[j - 1][k + 1] * t;
    double factor = 1.0;
    results[0] = Q[c.n][0];
    for (int j = 1; j <= r; j++)
    {
        factor *= c.n - (j - 1);
        for(int delta_num = j; delta_num > 0; delta_num--)
        {
            for (int k = 0; k <= delta_num - 1; k++)
                Q[c.n - j][k] = Q[c.n - j][k + 1] - Q[c.n - j][k];
        }
        results[j] = Q[c.n - j][0] * factor;
    }
}

void evaluate_derivatives_casteljau(Curve3D& c, unsigned int r, double t,
                                    Point3D* results)
{
    double t1 = 1.0 - t;
    Point3D Q[c.n + 1][c.n + 1];
    for (int k = 0; k <= c.n; k++) Q[0][k] = c.points[k];
    for (int j = 1; j <= c.n; j++)
        for (int k = 0; k <= c.n - j; k++)
            Q[j][k] = Q[j - 1][k] * t1 + Q[j - 1][k + 1] * t;
    double factor = 1.0;
    results[0] = Q[c.n][0];
    for (int j = 1; j <= r; j++)
    {
        factor *= c.n - (j - 1);
        for(int delta_num = j; delta_num > 0; delta_num--)
        {
            for (int k = 0; k <= delta_num - 1; k++)
                Q[c.n - j][k] = Q[c.n - j][k + 1] - Q[c.n - j][k];
        }
        results[j] = Q[c.n - j][0] * factor;
    }
}

void evaluate_derivatives_casteljau_pts(int n, const double* control_points, int r, double t, double* results)
{
    double t1 = 1.0 - t;
    double Q[n + 1][n + 1];
    for (int k = 0; k <= n; k++) Q[0][k] = control_points[k];
    for (int j = 1; j <= n; j++)
        for (int k = 0; k <= n - j; k++)
            Q[j][k] = Q[j - 1][k] * t1 + Q[j - 1][k + 1] * t;
    double factor = 1.0;
    results[0] = Q[n][0];
    for (int j = 1; j <= r; j++)
    {
        factor *= n - (j - 1);
        for(int delta_num = j; delta_num > 0; delta_num--)
        {
            for (int k = 0; k <= delta_num - 1; k++)
                Q[n - j][k] = Q[n - j][k + 1] - Q[n - j][k];
        }
        results[j] = Q[n - j][0] * factor;
    }
}

void evaluate_derivatives_casteljau_pts(int n, Point2D* control_points, int r, double t, Point2D* results)
{
    double t1 = 1.0 - t;
    Point2D Q[n + 1][n + 1];
    for (int k = 0; k <= n; k++) Q[0][k] = control_points[k];
    for (int j = 1; j <= n; j++)
        for (int k = 0; k <= n - j; k++)
            Q[j][k] = Q[j - 1][k] * t1 + Q[j - 1][k + 1] * t;
    double factor = 1.0;
    results[0] = Q[n][0];
    for (int j = 1; j <= r; j++)
    {
        factor *= n - (j - 1);
        for(int delta_num = j; delta_num > 0; delta_num--)
        {
            for (int k = 0; k <= delta_num - 1; k++)
                Q[n - j][k] = Q[n - j][k + 1] - Q[n - j][k];
        }
        results[j] = Q[n - j][0] * factor;
    }
}

void evaluate_derivatives_casteljau_pts(int n, Point3D* control_points, int r, double t, Point3D* results)
{
    double t1 = 1.0 - t;
    Point3D Q[n + 1][n + 1];
    for (int k = 0; k <= n; k++) Q[0][k] = control_points[k];
    for (int j = 1; j <= n; j++)
        for (int k = 0; k <= n - j; k++)
            Q[j][k] = Q[j - 1][k] * t1 + Q[j - 1][k + 1] * t;
    double factor = 1.0;
    results[0] = Q[n][0];
    for (int j = 1; j <= r; j++)
    {
        factor *= n - (j - 1);
        for(int delta_num = j; delta_num > 0; delta_num--)
        {
            for (int k = 0; k <= delta_num - 1; k++)
                Q[n - j][k] = Q[n - j][k + 1] - Q[n - j][k];
        }
        results[j] = Q[n - j][0] * factor;
    }
}

void evaluate_derivatives_wozny_chudy(Curve2D& c, unsigned int r, double t,
                                      Point2D* results)
{
    Curve2D curves[r + 1];
    curves[0] = c;
    for(int k = 1; k <= r; k++)
        curves[k] = curves[k - 1].derive_lower_degree();
    for(int k = 0; k <= r; k++)
        results[k] = curves[k].evaluate_wozny_chudy(t);
}

void evaluate_derivatives_wozny_chudy(Curve3D& c, unsigned int r, double t,
                                      Point3D* results)
{
    Curve3D curves[r + 1];
    curves[0] = c;
    for(int k = 1; k <= r; k++)
        curves[k] = curves[k - 1].derive_lower_degree();
    for(int k = 0; k <= r; k++)
        results[k] = curves[k].evaluate_wozny_chudy(t);
}

double evaluate_wozny_chudy_h_pts(int n, const double* h, const double* points)
{
    double q = points[0];
    for(int k = 1; k <= n; k++) q = q * (1 - h[k]) + points[k] * h[k];
    return q;
}

Point2D evaluate_wozny_chudy_h_pts(int n, double* h, Point2D* points)
{
    Point2D q = points[0];
    for(int k = 1; k <= n; k++) q = q * (1 - h[k]) + points[k] * h[k];
    return q;
}

Point3D evaluate_wozny_chudy_h_pts(int n, double* h, Point3D* points)
{
    Point3D q = points[0];
    for(int k = 1; k <= n; k++) q = q * (1 - h[k]) + points[k] * h[k];
    return q;
}

void derive_polynomial_lower(int n, const double* points, double* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++) result[i] = (points[i + 1] - points[i]) * n;
}

void derive_polynomial_lower(int n, Point2D* points, Point2D* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++) result[i] = (points[i + 1] - points[i]) * n;
}

void derive_polynomial_lower(int n, Point3D* points, Point3D* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++) result[i] = (points[i + 1] - points[i]) * n;
}

void derive_polynomial_keep(int n, const double* points, double* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++)
    {
        int i1 = i + 1;
        double delta = points[i1] - points[i];
        result[i] = result[i] + delta * (n - i);
        result[i1] = delta * i1;
    }
}

void derive_polynomial_keep(int n, Point2D* points, Point2D* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++)
    {
        int i1 = i + 1;
        Point2D delta = points[i1] - points[i];
        result[i] = result[i] + delta * (n - i);
        result[i1] = delta * i1;
    }
}

void derive_polynomial_keep(int n, Point3D* points, Point3D* result)
{
    if (n == 0) throw InvalidDerivation();
    for(int i = 0; i <= n - 1; i++)
    {
        int i1 = i + 1;
        Point3D delta = points[i1] - points[i];
        result[i] = result[i] + delta * (n - i);
        result[i1] = delta * i1;
    }
}

void test_polynomial_many_proper_1d(int n, int r, unsigned int t_num,
                                    unsigned int repetitions, unsigned int curve_count)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_dc = 0.0, elapsed_wc = 0.0, elapsed_wc_h = 0.0;
    double t;
    double pts_c[curve_count][n + 1], pts_wc[curve_count][n + 1], pts_wc_h[curve_count][n + 1];
    double h[n + 1][n + 1];

    auto w = new double[repetitions][100][101];
    for(int v = 0; v < repetitions; v++)
    {
        for (int m = 0; m < curve_count; m++)
        {
            for (int j = 0; j <= n; j++) w[v][m][j] = coordinate(gen);
        }
    }
    auto points_derivative = new double[100][101][101];
    auto points_derivative2 = new double[100][101][101];


    // de Casteljau
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;
            for(int m = 0; m < curve_count; m++)
            {
                evaluate_derivatives_casteljau_pts(n, w[v][m], r, t, pts_c[m]);
            }
        }
    }
    end_t = clock();
    elapsed_dc += (double)(end_t - start_t) / CLOCKS_PER_SEC;
    

    // WC_lower
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++)
                derive_polynomial_lower(n - k + 1, points_derivative[m][k - 1], points_derivative[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            for (int k = n; k >= n - r; k--) evaluate_h(k, t, h[k]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc[m][k] = evaluate_wozny_chudy_h_pts(n - k, h[n - k], points_derivative[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    //WC_keep
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {

        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative2[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++) points_derivative2[m][k][0] = 0.0;
            for(int k = 1; k <= r; k++) derive_polynomial_keep(n, points_derivative2[m][k - 1],
                                                               points_derivative2[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            evaluate_h(n, t, h[n]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc_h[m][k] = evaluate_wozny_chudy_h_pts(n, h[n], points_derivative2[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    std::cout << curve_count << "," << n << "," << r;
    std::cout << "," << elapsed_dc;
    std::cout << "," << elapsed_wc;
    std::cout << "," << elapsed_wc_h;
    std::cout << std::endl;

    delete[] w;
    delete[] points_derivative;
    delete[] points_derivative2;
}

void test_polynomial_many_proper_2d(int n, int r, unsigned int t_num,
                                  unsigned int repetitions, unsigned int curve_count)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_dc = 0.0, elapsed_wc = 0.0, elapsed_wc_h = 0.0;
    double t;
    Point2D pts_c[curve_count][n + 1], pts_wc[curve_count][n + 1], pts_wc_h[curve_count][n + 1];
    double h[n + 1][n + 1];

    auto w = new Point2D[repetitions][100][101];
    for(int v = 0; v < repetitions; v++)
    {
        for (int m = 0; m < curve_count; m++)
        {
            for (int j = 0; j <= n; j++) w[v][m][j] = Point2D(coordinate(gen), coordinate(gen));
        }
    }
    auto points_derivative = new Point2D[100][101][101];
    auto points_derivative2 = new Point2D[100][101][101];


    // de Casteljau
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;
            for(int m = 0; m < curve_count; m++)
            {
                evaluate_derivatives_casteljau_pts(n, w[v][m], r, t, pts_c[m]);
            }
        }
    }
    end_t = clock();
    elapsed_dc += (double)(end_t - start_t) / CLOCKS_PER_SEC;


    // WC_lower
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++)
                derive_polynomial_lower(n - k + 1, points_derivative[m][k - 1], points_derivative[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            for (int k = n; k >= n - r; k--) evaluate_h(k, t, h[k]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc[m][k] = evaluate_wozny_chudy_h_pts(n - k, h[n - k], points_derivative[m][k]);
            }

        }
    }
    end_t = clock();
    elapsed_wc += (double)(end_t - start_t) / CLOCKS_PER_SEC;

     //WC_keep
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {

        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative2[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++) points_derivative2[m][k][0] = Point2D();
            for(int k = 1; k <= r; k++) derive_polynomial_keep(n, points_derivative2[m][k - 1],
                                                               points_derivative2[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            evaluate_h(n, t, h[n]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc_h[m][k] = evaluate_wozny_chudy_h_pts(n, h[n], points_derivative2[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    std::cout << curve_count << "," << n << "," << r;
    std::cout << "," << elapsed_dc;
    std::cout << "," << elapsed_wc;
    std::cout << "," << elapsed_wc_h;
    std::cout << std::endl;

    delete[] w;
    delete[] points_derivative;
    delete[] points_derivative2;
}

void test_polynomial_many_proper_2d_no_dc(int n, int r, unsigned int t_num,
                                          unsigned int repetitions, unsigned int curve_count)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_dc = 0.0, elapsed_wc = 0.0, elapsed_wc_h = 0.0;
    double t;
    Point2D pts_c[curve_count][n + 1], pts_wc[curve_count][n + 1], pts_wc_h[curve_count][n + 1];
    double h[n + 1][n + 1];

    auto w = new Point2D[repetitions][100][101];
    for(int v = 0; v < repetitions; v++)
    {
        for (int m = 0; m < curve_count; m++)
        {
            for (int j = 0; j <= n; j++) w[v][m][j] = Point2D(coordinate(gen), coordinate(gen));
        }
    }
    auto points_derivative = new Point2D[100][101][101];
    auto points_derivative2 = new Point2D[100][101][101];


    // WC_lower
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++)
                derive_polynomial_lower(n - k + 1, points_derivative[m][k - 1], points_derivative[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            for (int k = n; k >= n - r; k--) evaluate_h(k, t, h[k]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc[m][k] = evaluate_wozny_chudy_h_pts(n - k, h[n - k], points_derivative[m][k]);
            }

        }
    }
    end_t = clock();
    elapsed_wc += (double)(end_t - start_t) / CLOCKS_PER_SEC;

     //WC_keep
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {

        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative2[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++) points_derivative2[m][k][0] = Point2D();
            for(int k = 1; k <= r; k++) derive_polynomial_keep(n, points_derivative2[m][k - 1],
                                                               points_derivative2[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            evaluate_h(n, t, h[n]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc_h[m][k] = evaluate_wozny_chudy_h_pts(n, h[n], points_derivative2[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    std::cout << curve_count << "," << n << "," << r;
    std::cout << "," << "n/a";
    std::cout << "," << elapsed_wc;
    std::cout << "," << elapsed_wc_h;
    std::cout << std::endl;

    delete[] w;
    delete[] points_derivative;
    delete[] points_derivative2;
}


void test_polynomial_many_proper_3d(int n, int r, unsigned int t_num,
                                  unsigned int repetitions, unsigned int curve_count)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_dc = 0.0, elapsed_wc = 0.0, elapsed_wc_h = 0.0;
    double t;
    Point3D pts_c[curve_count][n + 1], pts_wc[curve_count][n + 1], pts_wc_h[curve_count][n + 1];
    double h[n + 1][n + 1];

    auto w = new Point3D[repetitions][100][101];
    for(int v = 0; v < repetitions; v++)
    {
        for (int m = 0; m < curve_count; m++)
        {
            for (int j = 0; j <= n; j++) w[v][m][j] = Point3D(coordinate(gen), 
                                                              coordinate(gen),
                                                              coordinate(gen));
        }
    }
    auto points_derivative = new Point3D[100][101][101];
    auto points_derivative2 = new Point3D[100][101][101];


    // de Casteljau
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;
            for(int m = 0; m < curve_count; m++)
            {
                evaluate_derivatives_casteljau_pts(n, w[v][m], r, t, pts_c[m]);
            }
        }
    }
    end_t = clock();
    elapsed_dc += (double)(end_t - start_t) / CLOCKS_PER_SEC;


    // WC_lower
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++)
                derive_polynomial_lower(n - k + 1, points_derivative[m][k - 1], points_derivative[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            for (int k = n; k >= n - r; k--) evaluate_h(k, t, h[k]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc[m][k] = evaluate_wozny_chudy_h_pts(n - k, h[n - k], points_derivative[m][k]);
            }

        }
    }
    end_t = clock();
    elapsed_wc += (double)(end_t - start_t) / CLOCKS_PER_SEC;

     //WC_keep
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {

        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative2[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++) points_derivative2[m][k][0] = Point3D();
            for(int k = 1; k <= r; k++) derive_polynomial_keep(n, points_derivative2[m][k - 1],
                                                               points_derivative2[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            evaluate_h(n, t, h[n]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc_h[m][k] = evaluate_wozny_chudy_h_pts(n, h[n], points_derivative2[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    std::cout << curve_count << "," << n << "," << r;
    std::cout << "," << elapsed_dc;
    std::cout << "," << elapsed_wc;
    std::cout << "," << elapsed_wc_h;
    std::cout << std::endl;

    delete[] w;
    delete[] points_derivative;
    delete[] points_derivative2;
}

void test_polynomial_many_proper_3d_no_dc(int n, int r, unsigned int t_num,
                                  unsigned int repetitions, unsigned int curve_count)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_dc = 0.0, elapsed_wc = 0.0, elapsed_wc_h = 0.0;
    double t;
    Point3D pts_c[curve_count][n + 1], pts_wc[curve_count][n + 1], pts_wc_h[curve_count][n + 1];
    double h[n + 1][n + 1];

    auto w = new Point3D[repetitions][100][101];
    for(int v = 0; v < repetitions; v++)
    {
        for (int m = 0; m < curve_count; m++)
        {
            for (int j = 0; j <= n; j++) w[v][m][j] = Point3D(coordinate(gen), 
                                                              coordinate(gen),
                                                              coordinate(gen));
        }
    }
    auto points_derivative = new Point3D[100][101][101];
    auto points_derivative2 = new Point3D[100][101][101];

    // WC_lower
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {
        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++)
                derive_polynomial_lower(n - k + 1, points_derivative[m][k - 1], points_derivative[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            for (int k = n; k >= n - r; k--) evaluate_h(k, t, h[k]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc[m][k] = evaluate_wozny_chudy_h_pts(n - k, h[n - k], points_derivative[m][k]);
            }

        }
    }
    end_t = clock();
    elapsed_wc += (double)(end_t - start_t) / CLOCKS_PER_SEC;

     //WC_keep
    start_t = clock();
    for(int v = 0; v < repetitions; v++)
    {

        for(int m = 0; m < curve_count; m++)
        {
            for(int k = 0; k <= n; k++) points_derivative2[m][0][k] = w[v][m][k];
            for(int k = 1; k <= r; k++) points_derivative2[m][k][0] = Point3D();
            for(int k = 1; k <= r; k++) derive_polynomial_keep(n, points_derivative2[m][k - 1],
                                                               points_derivative2[m][k]);
        }

        for(int i = 0; i <= t_num; i++)
        {
            t = (1. * i) / t_num;

            evaluate_h(n, t, h[n]);

            for(int m = 0; m < curve_count; m++)
            {
                for (int k = 0; k <= r; k++)
                    pts_wc_h[m][k] = evaluate_wozny_chudy_h_pts(n, h[n], points_derivative2[m][k]);
            }
        }
    }
    end_t = clock();
    elapsed_wc_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;

    std::cout << curve_count << "," << n << "," << r;
    std::cout << "," << "n/a";
    std::cout << "," << elapsed_wc;
    std::cout << "," << elapsed_wc_h;
    std::cout << std::endl;

    delete[] w;
    delete[] points_derivative;
    delete[] points_derivative2;
}

double good_digits(double value, double correct)
{
    if(value == correct) return 53.0 * std::log10(2.0);
    return -std::log10( std::abs((value - correct)/ correct));

}

void test_polynomial_num(int n, unsigned int r, unsigned int t_num,
                         unsigned int repetitions)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    double t;
    Point2D w[n + 1];
    Point2D pts_c[n + 1], pts_wc[n + 1], pts_wc_h[n + 1];
    double digits_min[n + 1], digits_max[n + 1], digits_sum[n + 1];
    for(int k = 0; k <= r; k++)
    {
        digits_min[k] = 53.0 * std::log10(2.0);
        digits_max[k] = 0.0;
        digits_sum[k] = 0.0;
    }
    Curve2D curves[repetitions];
    for(int i = 0; i <= t_num; i++)
    {
        t = (1. * i) / t_num;

        double h[n + 1];
        evaluate_h(n, t, h);
        for(int v = 0; v < repetitions; v++) {
            for (int j = 0; j <= n; j++)
                w[j] = Point2D(coordinate(gen), coordinate(gen));
            Curve2D c = Curve2D(n, w);

            evaluate_derivatives_casteljau(c, r, t, pts_c);
            evaluate_derivatives_wozny_chudy(c, r, t, pts_wc);
            for (int k = 0; k <= r; k++) {
                double digits = (good_digits(pts_c[k].x, pts_wc[k].x) + good_digits(pts_c[k].y, pts_wc[k].y)) / 2.0;
                digits_min[k] = std::min(digits, digits_min[k]);
                digits_max[k] = std::max(digits, digits_max[k]);
                digits_sum[k] += digits;
            }
        }
    }
    for(int k = 0; k <= r; k++)
    {
        std::cout << n << "," << r << "," << k << ",";
        std::cout << digits_min[k] << ",";
        std::cout << digits_sum[k] / (repetitions * t_num) << ",";
        std::cout << digits_max[k] << std::endl;
    }
}

void test_rational_2d(int n, unsigned int r, unsigned int t_num,
                unsigned int repetitions)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_nf = 0.0, elapsed_sf = 0.0, elapsed_g = 0.0, elapsed_h = 0.0;
    double t;
    double weights[n + 1];
    Point2D w[n + 1];
    Point2D pts_sf[n + 1], pts_nf[n + 1], pts_h[n + 1], pts_g[n + 1];
    RatCurve2D curves[repetitions];
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= n; i++)
        {
            w[i] = Point2D(coordinate(gen), coordinate(gen));
            weights[i] = weight(gen);
        }
        curves[v] = RatCurve2D(n, w, weights);
    }
    for(int i = 0; i <= t_num; i++)
    {
        t = (1. * i) / t_num;
        
        for(int v = 0; v < repetitions; v++)
        {
            if(r == 2)
            {
                start_t = clock();
                curves[v].derive_simple_floater(t, pts_sf, false);
                end_t = clock();
                elapsed_sf += (double)(end_t - start_t) / CLOCKS_PER_SEC;

                start_t = clock();
                curves[v].derive_new_floater(t, pts_nf, false);
                end_t = clock();
                elapsed_nf += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            }
            else if(r == 1)
            {
                start_t = clock();
                curves[v].derive_simple_floater(t, pts_sf, true);
                end_t = clock();
                elapsed_sf += (double)(end_t - start_t) / CLOCKS_PER_SEC;

                start_t = clock();
                curves[v].derive_new_floater(t, pts_nf, true);
                end_t = clock();
                elapsed_nf += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            }
            
            start_t = clock();
            curves[v].derive_new_geom(t, pts_g, r);
            end_t = clock();
            elapsed_g += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            

            start_t = clock();
            curves[v].derive_new_h(t, pts_h, r);
            end_t = clock();
            elapsed_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;
        }

    }
    std::cout << n << "," << r << ",";
    if(r <= 2)
    {
        std::cout << elapsed_sf << ",";
        std::cout << elapsed_nf << ",";
    }
    std::cout << elapsed_g << ",";
    std::cout << elapsed_h << std::endl;
}

void test_rational_3d(int n, unsigned int r, unsigned int t_num,
                unsigned int repetitions)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    clock_t start_t, end_t;
    double elapsed_nf = 0.0, elapsed_sf = 0.0, elapsed_g = 0.0, elapsed_h = 0.0;
    double t;
    double weights[n + 1];
    Point3D w[n + 1];
    Point3D pts_sf[n + 1], pts_nf[n + 1], pts_h[n + 1], pts_g[n + 1];
    RatCurve3D curves[repetitions];
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= n; i++)
        {
            w[i] = Point3D(coordinate(gen), coordinate(gen), coordinate(gen));
            weights[i] = weight(gen);
        }
        curves[v] = RatCurve3D(n, w, weights);
    }
    for(int i = 0; i <= t_num; i++)
    {
        t = (1. * i) / t_num;
        
        for(int v = 0; v < repetitions; v++)
        {
            if(r == 2)
            {
                start_t = clock();
                curves[v].derive_simple_floater(t, pts_sf, false);
                end_t = clock();
                elapsed_sf += (double)(end_t - start_t) / CLOCKS_PER_SEC;

                start_t = clock();
                curves[v].derive_new_floater(t, pts_nf, false);
                end_t = clock();
                elapsed_nf += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            }
            else if(r == 1)
            {
                start_t = clock();
                curves[v].derive_simple_floater(t, pts_sf, true);
                end_t = clock();
                elapsed_sf += (double)(end_t - start_t) / CLOCKS_PER_SEC;

                start_t = clock();
                curves[v].derive_new_floater(t, pts_nf, true);
                end_t = clock();
                elapsed_nf += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            }
            
            start_t = clock();
            curves[v].derive_new_geom(t, pts_g, r);
            end_t = clock();
            elapsed_g += (double)(end_t - start_t) / CLOCKS_PER_SEC;
            

            start_t = clock();
            curves[v].derive_new_h(t, pts_h, r);
            end_t = clock();
            elapsed_h += (double)(end_t - start_t) / CLOCKS_PER_SEC;
        }

    }
    std::cout << n << "," << r << ",";
    if(r <= 2)
    {
        std::cout << elapsed_sf << ",";
        std::cout << elapsed_nf << ",";
    }
    std::cout << elapsed_g << ",";
    std::cout << elapsed_h << std::endl;
}

void test_rational_num(int n, unsigned int r, unsigned int t_num,
                         unsigned int repetitions)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> coordinate(-1.0, 1.0);
    std::uniform_real_distribution<double> weight(0.01, 2.0);
    double t;
    double weights[n + 1];
    Point2D w[n + 1];
    Point2D pts_sf[n + 1], pts_nf[n + 1], pts_h[n + 1], pts_g[n + 1];
    RatCurve2D curves[repetitions];
    for(int v = 0; v < repetitions; v++)
    {
        for(int i = 0; i <= n; i++)
        {
            w[i] = Point2D(coordinate(gen), coordinate(gen));
            weights[i] = weight(gen);
        }
        curves[v] = RatCurve2D(n, w, weights);
    }

    double digits_min_nf[n + 1], digits_max_nf[n + 1], digits_sum_nf[n + 1];
    double digits_min_g[n + 1], digits_max_g[n + 1], digits_sum_g[n + 1];
    double digits_min_h[n + 1], digits_max_h[n + 1], digits_sum_h[n + 1];
    for(int k = 0; k <= r; k++)
    {
        digits_min_nf[k] = 53.0 * std::log10(2.0);
        digits_max_nf[k] = 0.0;
        digits_sum_nf[k] = 0.0;
        digits_min_g[k] = 53.0 * std::log10(2.0);
        digits_max_g[k] = 0.0;
        digits_sum_g[k] = 0.0;
        digits_min_h[k] = 53.0 * std::log10(2.0);
        digits_max_h[k] = 0.0;
        digits_sum_h[k] = 0.0;
    }
    for(int i = 0; i <= t_num; i++)
    {
        t = (1. * i) / t_num;

        for(int v = 0; v < repetitions; v++)
        {
            if(r == 2)
            {
                curves[v].derive_simple_floater(t, pts_sf, false);
                curves[v].derive_new_floater(t, pts_nf, false);
            }
            else if(r == 1)
            {
                curves[v].derive_simple_floater(t, pts_sf, true);
                curves[v].derive_new_floater(t, pts_nf, true);
            }

            curves[v].derive_new_geom(t, pts_g, r);
            curves[v].derive_new_h(t, pts_h, r);
            if(r <= 2)
            {
                for (int k = 0; k <= r; k++) {
                    double digits = (good_digits(pts_nf[k].x, pts_sf[k].x) + good_digits(pts_nf[k].y, pts_sf[k].y)) / 2.0;
                    digits_min_nf[k] = std::min(digits, digits_min_nf[k]);
                    digits_max_nf[k] = std::max(digits, digits_max_nf[k]);
                    digits_sum_nf[k] += digits;

                    digits = (good_digits(pts_g[k].x, pts_sf[k].x) + good_digits(pts_g[k].y, pts_sf[k].y)) / 2.0;
                    digits_min_g[k] = std::min(digits, digits_min_g[k]);
                    digits_max_g[k] = std::max(digits, digits_max_g[k]);
                    digits_sum_g[k] += digits;

                    digits = (good_digits(pts_h[k].x, pts_sf[k].x) + good_digits(pts_h[k].y, pts_sf[k].y)) / 2.0;
                    digits_min_h[k] = std::min(digits, digits_min_h[k]);
                    digits_max_h[k] = std::max(digits, digits_max_h[k]);
                    digits_sum_h[k] += digits;
                }
            }
            else
            {
                for (int k = 0; k <= r; k++) {
                    double digits = (good_digits(pts_h[k].x, pts_g[k].x) + good_digits(pts_h[k].y, pts_g[k].y)) / 2.0;
                    digits_min_h[k] = std::min(digits, digits_min_h[k]);
                    digits_max_h[k] = std::max(digits, digits_max_h[k]);
                    digits_sum_h[k] += digits;
                }
            }
        }

    }
    for(int k = 0; k <= r; k++)
    {
        if(r <= 2)
        {
            std::cout << n << "," << r << "," << k << ",";
            std::cout << "NF,";
            std::cout << digits_min_nf[k] << ",";
            std::cout << digits_sum_nf[k] / (repetitions * t_num) << ",";
            std::cout << digits_max_nf[k] << std::endl;

            std::cout << n << "," << r << "," << k << ",";
            std::cout << "G,";
            std::cout << digits_min_g[k] << ",";
            std::cout << digits_sum_g[k] / (repetitions * t_num) << ",";
            std::cout << digits_max_g[k] << std::endl;

            std::cout << n << "," << r << "," << k << ",";
            std::cout << "H,";
            std::cout << digits_min_h[k] << ",";
            std::cout << digits_sum_h[k] / (repetitions * t_num) << ",";
            std::cout << digits_max_h[k] << std::endl;
        }
        else
        {
            std::cout << n << "," << r << "," << k << ",";
            std::cout << digits_min_h[k] << ",";
            std::cout << digits_sum_h[k] / (repetitions * t_num) << ",";
            std::cout << digits_max_h[k] << std::endl;
        }
    }
}

int main()
{
    int n_candidates[8] = {2, 3, 4, 5,
                           10, 20, 30, 50};
    int r_candidates[3] = {1, 2, 3};
    int m_candidates[3] = {2,5, 10};
    int repetitions = 1000; // for 1000 reps, it takes 3.5 hrs
    int t_num = 500;

    /*
    std::cout << "n,r,k,min,mean,max" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r) test_polynomial_num(n, r, t_num, repetitions);
        }
    }

    std::cout << "Table 2.1" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
     for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r) test_polynomial_many_proper_2d(n, r, t_num, repetitions, 1);
        }
        if(n > 3) test_polynomial_many_proper_2d(n, n, t_num, repetitions, 1);
    }

    std::cout << "Table 2.2" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
    for(int n : n_candidates)
    {
        if(n >= 20)
        {
            for(int r : r_candidates)
            {
                if(n >= r) test_polynomial_many_proper_1d(n, r, t_num, repetitions, 1);
            }
            if(n > 3) test_polynomial_many_proper_1d(n, n, t_num, repetitions, 1);
        }
    }

    std::cout << "Table 2.3" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
    for(int n : {2, 3, 5, 10, 50})
    {
        for(int r : {2, n})
        {
            for(int m : {2, 5, 10})
            {
                test_polynomial_many_proper_2d(n, r, t_num, repetitions, m);
            }
        }
    }

    std::cout << "n,r,simpleFloater,newFloater,newGeom,newH" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r <= 2) test_rational_2d(n, r, t_num, repetitions);
        }
    }

    std::cout << "n,r,k,method,min,mean,max" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r <= 2) test_rational_num(n, r, t_num, repetitions);
        }
    }

    std::cout << "n,r,newGeom,newH" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r > 2) test_rational_2d(n, r, t_num, repetitions);
        }
    }

    std::cout << "n,r,k,diff" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r > 2) test_rational_num(n, r, t_num, repetitions);
        }
    }
    */

    /*
    std::cout << "Table 2.1 (3D)" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
     for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r) test_polynomial_many_proper_3d(n, r, t_num, repetitions, 1);
        }
        if(n > 3) test_polynomial_many_proper_3d(n, n, t_num, repetitions, 1);
    }
    */
    /*
    std::cout << "Table 3.1 (3D)" << std::endl;
    std::cout << "n,r,simpleFloater,newFloater,newGeom,newH" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r <= 2) test_rational_3d(n, r, t_num, repetitions);
        }
    }
    */
    /*
    std::cout << "Table 3.2 (3D)" << std::endl;
    std::cout << "n,r,newGeom,newH" << std::endl;
    for(int n : n_candidates)
    {
        for(int r : n_candidates)
        {
            if(n >= r && r > 2) test_rational_3d(n, r, t_num, repetitions);
        }
    }
    */
    int n_big[3] = {100, 200, 300};

    /*
    std::cout << "Table 2.1 (2D) extended" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
     for(int n : n_big)
    {
        for(int r : r_candidates)
        {
            if(n >= r) test_polynomial_many_proper_2d(n, r, t_num, repetitions, 1);
        }
        if(n > 3)
        {
            if(n == 300) test_polynomial_many_proper_2d_no_dc(n, n, t_num, repetitions, 1);
            else test_polynomial_many_proper_2d(n, n, t_num, repetitions, 1);
        }
    }
    std::cout << "Table 2.1 (3D) extended" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
     for(int n : n_big)
    {
        for(int r : r_candidates)
        {
            if(n >= r) test_polynomial_many_proper_3d(n, r, t_num, repetitions, 1);
        }
        if(n > 3)
        {
            if(n == 300) test_polynomial_many_proper_3d_no_dc(n, n, t_num, repetitions, 1);
            else test_polynomial_many_proper_3d(n, n, t_num, repetitions, 1);
        }
    }
    */
    
    /*
    std::cout << "Table 3.1 (2D) extended" << std::endl;
    std::cout << "n,r,simpleFloater,newFloater,newGeom,newH" << std::endl;
    for(int n : n_big)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r <= 2) test_rational_2d(n, r, t_num, repetitions);
        }
    }
    std::cout << "Table 3.1 (3D) extended" << std::endl;
    std::cout << "n,r,simpleFloater,newFloater,newGeom,newH" << std::endl;
    for(int n : n_big)
    {
        for(int r : r_candidates)
        {
            if(n >= r && r <= 2) test_rational_3d(n, r, t_num, repetitions);
        }
    }
    */

    int n_big2[2] = {100, 200};

    /*
    std::cout << "Table 3.2 (2D) extended" << std::endl;
    std::cout << "n,r,newGeom,newH" << std::endl;
    for(int n : n_big2)
    {
        for(int r : n_candidates)
        {
            if(n >= r && r > 2) test_rational_2d(n, r, t_num, repetitions);
        }
        if(n > 3) test_rational_2d(n, n, t_num, repetitions);
    }
    */
    /*
    std::cout << "Table 3.2 (3D) extended" << std::endl;
    std::cout << "n,r,newGeom,newH" << std::endl;
    for(int n : n_big2)
    {
        for(int r : n_candidates)
        {
            if(n >= r && r > 2) test_rational_3d(n, r, t_num, repetitions);
        }
        if(n > 3) test_rational_3d(n, n, t_num, repetitions);
    }
    */
    
    std::cout << "Table 2.2 extended" << std::endl;
    std::cout << "m,n,r,deCasteljau,WoznyChudy,WoznyChudyH" << std::endl;
    for(int n : n_big)
    {
        if(n >= 20)
        {
            for(int r : r_candidates)
            {
                if(n >= r) test_polynomial_many_proper_1d(n, r, t_num, repetitions, 1);
            }
            if(n > 3) test_polynomial_many_proper_1d(n, n, t_num, repetitions, 1);
        }
    }

    return 0;
}