#include <iostream>
#include <iomanip>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <limits>
#include <stdexcept>

// POLYNOMIAL ALGORITHMS
//
// A polynomial p(x) = p[0] + p[1]x + ... + p[d]x^d is stored
// in a dynamically allocated array of d + 1 coefficients.
// The zero polynomial is represented by degree 0 and {0.0}.
//
// Each operation producing a polynomial allocates a NEW array.
// The caller owns that array and must release it with delete[].
// Output pointers must not alias any input array.
//
// This exercise covers differentiation, addition, multiplication,
// copying, scalar multiplication, Horner evaluation, synthetic
// division, and numerical root finding.

int main();

void diff(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q
);
void multiply(
    double       const p[],
    unsigned int const deg_p,
    double       const q[],
    unsigned int const deg_q,
    double            *&r,
    unsigned int      &deg_r
);
void add(
    double       const p[],
    unsigned int const deg_p,
    double       const q[],
    unsigned int const deg_q,
    double            *&r,
    unsigned int      &deg_r
);
void copy(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q
);
double eval(
    double       const p[],
    unsigned int const deg_p,
    double       const x
);
void scalar_multiply(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q,
    double       const c
);
double root_divide(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q,
    double       const c
);
double newton(
    double       const p[],
    unsigned int const deg_p
);
void real_roots(
    double       const p[],
    unsigned int const deg_p
);
void print(
    double       const p[],
    unsigned int const deg_p
);

// Remove exactly zero leading coefficients. Floating-point
// calculations can produce tiny nonzero coefficients; do not
// silently discard those without an application-specific tolerance.
void normalize_degree(
    double       const p[],
    unsigned int      &degree
) {
    assert( p != nullptr );
    while ( degree > 0 && p[degree] == 0.0 ) {
        --degree;
    }
}

int main() {
    double const p1[]{ 0.0 };
    double const p2[]{ 3.2 };
    double const p3[]{ 3.2, 4.7, -1.2 };
    double const p4[]{ -7.2, 0.0, 5.2, 0.2 };
    double const p5[]{ -8.3, 4.7, -5.2, -0.2 };
    double const p6[]{ 7.2, 0.0, -5.2, -0.2 };

    double const *polynomials[]{ p1, p2, p3, p4, p5, p6 };
    unsigned int const degrees[]{ 0, 0, 2, 3, 3, 3 };
    unsigned int const count{ 6 };

    for ( unsigned int i{ 0 }; i < count; ++i ) {
        std::cout << "p" << i + 1 << "(x) = ";
        print( polynomials[i], degrees[i] );
    }

    double *result{ nullptr };
    unsigned int degree{ 0 };

    std::cout << "\nDifferentiation\n";
    for ( unsigned int i{ 0 }; i < count; ++i ) {
        diff( polynomials[i], degrees[i], result, degree );
        print( result, degree );
        delete[] result;
        result = nullptr;
    }

    std::cout << "\nAddition\n";
    for ( unsigned int i{ 0 }; i < count; ++i ) {
        for ( unsigned int j{ 0 }; j < count; ++j ) {
            add( polynomials[i], degrees[i],
                 polynomials[j], degrees[j], result, degree );
            print( result, degree );
            delete[] result;
            result = nullptr;
        }
        std::cout << std::endl;
    }

    std::cout << "\nMultiplication\n";
    for ( unsigned int i{ 0 }; i < count; ++i ) {
        for ( unsigned int j{ 0 }; j < count; ++j ) {
            multiply( polynomials[i], degrees[i],
                      polynomials[j], degrees[j], result, degree );
            print( result, degree );
            delete[] result;
            result = nullptr;
        }
        std::cout << std::endl;
    }

    std::cout << "\nSynthetic division by x - 3\n";
    double const remainder{ root_divide( p6, 3, result, degree, 3.0 ) };
    std::cout << "Remainder: " << remainder << "\nQuotient: ";
    print( result, degree );
    assert( std::abs( remainder - eval( p6, 3, 3.0 ) ) < 1e-10 );
    delete[] result;
    result = nullptr;

    std::cout << "\nNewton's method\n";
    try {
        double const root{ newton( p6, 3 ) };
        std::cout << "Approximate root: " << std::setprecision( 15 )
                  << root << "\n";
        double const rem{ root_divide( p6, 3, result, degree, root ) };
        std::cout << "Remainder: " << rem << "\nQuotient: ";
        print( result, degree );
        delete[] result;
        result = nullptr;
    } catch ( std::runtime_error const &error ) {
        std::cout << error.what() << std::endl;
    }

    std::cout << "\nReal roots found by repeated Newton iterations\n";
    real_roots( p6, 3 );

    // Additional tests for cancellation, constant polynomials,
    // multiplication accumulation, and synthetic division.
    add( p4, 3, p6, 3, result, degree );
    assert( degree == 0 && result[0] == 0.0 );
    delete[] result;
    result = nullptr;

    multiply( p3, 2, p3, 2, result, degree );
    assert( degree == 4 );
    assert( std::abs( eval( result, degree, 2.0 )
                      - eval( p3, 2, 2.0 )*eval( p3, 2, 2.0 ) ) < 1e-10 );
    delete[] result;
    result = nullptr;

    diff( p2, 0, result, degree );
    assert( degree == 0 && result[0] == 0.0 );
    delete[] result;
    result = nullptr;

    scalar_multiply( p3, 2, result, degree, 0.0 );
    assert( degree == 0 && result[0] == 0.0 );
    delete[] result;
    result = nullptr;

    assert( eval( p2, 0, 123.0 ) == 3.2 );
    return 0;
}

void diff(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q
) {
    assert( p != nullptr );
    assert( deg_p == 0 || p[deg_p] != 0.0 );
    deg_q = (deg_p == 0) ? 0 : deg_p - 1;
    q = new double[deg_q + 1]{};
    for ( unsigned int k{ 1 }; k <= deg_p; ++k ) {
        q[k - 1] = k*p[k];
    }
    normalize_degree( q, deg_q );
}

void multiply(
    double       const p[],
    unsigned int const deg_p,
    double       const q[],
    unsigned int const deg_q,
    double            *&r,
    unsigned int      &deg_r
) {
    assert( p != nullptr && q != nullptr );
    assert( deg_p == 0 || p[deg_p] != 0.0 );
    assert( deg_q == 0 || q[deg_q] != 0.0 );
    assert( deg_p <= std::numeric_limits<unsigned int>::max() - deg_q );
    if ( (deg_p == 0 && p[0] == 0.0)
      || (deg_q == 0 && q[0] == 0.0) ) {
        deg_r = 0;
        r = new double[1]{ 0.0 };
        return;
    }
    deg_r = deg_p + deg_q;
    r = new double[deg_r + 1]{};
    for ( unsigned int i{ 0 }; i <= deg_p; ++i ) {
        for ( unsigned int j{ 0 }; j <= deg_q; ++j ) {
            r[i + j] += p[i]*q[j];
        }
    }
    normalize_degree( r, deg_r );
}

void add(
    double       const p[],
    unsigned int const deg_p,
    double       const q[],
    unsigned int const deg_q,
    double            *&r,
    unsigned int      &deg_r
) {
    assert( p != nullptr && q != nullptr );
    assert( deg_p == 0 || p[deg_p] != 0.0 );
    assert( deg_q == 0 || q[deg_q] != 0.0 );
    deg_r = (deg_p > deg_q) ? deg_p : deg_q;
    r = new double[deg_r + 1]{};
    for ( unsigned int k{ 0 }; k <= deg_p; ++k ) {
        r[k] += p[k];
    }
    for ( unsigned int k{ 0 }; k <= deg_q; ++k ) {
        r[k] += q[k];
    }
    normalize_degree( r, deg_r );
}

void copy(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q
) {
    assert( p != nullptr );
    deg_q = deg_p;
    q = new double[deg_q + 1]{};
    for ( unsigned int k{ 0 }; k <= deg_q; ++k ) {
        q[k] = p[k];
    }
}

double eval(
    double       const p[],
    unsigned int const deg_p,
    double       const x
) {
    assert( p != nullptr );
    double result{ p[deg_p] };
    for ( unsigned int k{ deg_p }; k > 0; --k ) {
        result = result*x + p[k - 1];
    }
    return result;
}

void scalar_multiply(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q,
    double       const c
) {
    assert( p != nullptr );
    if ( c == 0.0 ) {
        deg_q = 0;
        q = new double[1]{ 0.0 };
        return;
    }
    deg_q = deg_p;
    q = new double[deg_q + 1]{};
    for ( unsigned int k{ 0 }; k <= deg_q; ++k ) {
        q[k] = c*p[k];
    }
    normalize_degree( q, deg_q );
}

// Synthetic division: p(x) = (x - c)q(x) + remainder.
// For a constant polynomial, the quotient is zero.
double root_divide(
    double       const p[],
    unsigned int const deg_p,
    double            *&q,
    unsigned int      &deg_q,
    double       const c
) {
    assert( p != nullptr );
    if ( deg_p == 0 ) {
        deg_q = 0;
        q = new double[1]{ 0.0 };
        return p[0];
    }
    deg_q = deg_p - 1;
    q = new double[deg_q + 1]{};
    q[deg_q] = p[deg_p];
    for ( unsigned int k{ deg_q }; k > 0; --k ) {
        q[k - 1] = p[k] + c*q[k];
    }
    double const remainder{ p[0] + c*q[0] };
    normalize_degree( q, deg_q );
    return remainder;
}

// Try Newton's method from several starting points. A failure
// means only that these attempts did not find a root; it does
// NOT prove that the polynomial has no real roots.
double newton(
    double       const p[],
    unsigned int const deg_p
) {
    assert( p != nullptr );
    if ( deg_p == 0 ) {
        throw std::runtime_error{ "A constant polynomial has no isolated root" };
    }
    double *dp{ nullptr };
    unsigned int deg_dp{ 0 };
    diff( p, deg_p, dp, deg_dp );
    // A small fixed set of initial guesses makes demonstrations
    // reproducible without requiring a random-number generator.
    double const starts[]{ 0.0, 1.0, -1.0, 2.0, -2.0,
                           4.0, -4.0, 8.0, -8.0,
                           16.0, -16.0, 32.0, -32.0 };
    for ( double const start : starts ) {
        double x{ start };
        for ( unsigned int iteration{ 0 }; iteration < 1000; ++iteration ) {
            double const fx{ eval( p, deg_p, x ) };
            double const dfx{ eval( dp, deg_dp, x ) };
            if ( !std::isfinite( fx ) || !std::isfinite( dfx ) ) {
                break;
            }
            // Test the residual relative to a simple coefficient scale.
            double scale{ 1.0 };
            double power{ 1.0 };
            for ( unsigned int k{ 0 }; k <= deg_p; ++k ) {
                scale += std::abs( p[k] )*power;
                power *= std::abs( x );
            }
            if ( std::abs( fx ) <= 1e-12*scale ) {
                delete[] dp;
                return x;
            }
            if ( dfx == 0.0 ) {
                break;
            }
            double const next{ x - fx/dfx };
            if ( !std::isfinite( next ) || next == x ) {
                break;
            }
            x = next;
        }
    }
    delete[] dp;
    throw std::runtime_error{ "Newton's method did not find a real root" };
}

// Find roots by Newton's method and repeatedly divide by (x - root).
// Numerical deflation introduces rounding errors, and Newton's
// method can miss roots. The printed roots are candidates, not a
// guarantee that every real root has been found.
void real_roots(
    double       const p[],
    unsigned int const deg_p
) {
    assert( p != nullptr );
    if ( deg_p == 0 ) {
        std::cout << (p[0] == 0.0 ? "Every real number is a root"
                                      : "No roots (nonzero constant)")
                  << std::endl;
        return;
    }
    double *current{ nullptr };
    unsigned int degree{ 0 };
    copy( p, deg_p, current, degree );
    bool first{ true };
    while ( degree > 0 ) {
        double root{ 0.0 };
        try {
            root = newton( current, degree );
        } catch ( std::runtime_error const & ) {
            break;
        }
        double *quotient{ nullptr };
        unsigned int quotient_degree{ 0 };
        double const remainder{
            root_divide( current, degree, quotient, quotient_degree, root )
        };
        double scale{ 1.0 };
        double power{ 1.0 };
        for ( unsigned int k{ 0 }; k <= degree; ++k ) {
            scale += std::abs( current[k] )*power;
            power *= std::abs( root );
        }
        if ( !std::isfinite( remainder )
          || std::abs( remainder ) > 1e-9*scale ) {
            delete[] quotient;
            break;
        }
        if ( !first ) {
            std::cout << ", ";
        }
        first = false;
        std::cout << std::setprecision( 15 ) << root;
        delete[] current;
        current = quotient;
        degree = quotient_degree;
    }
    if ( first ) {
        std::cout << "No real roots found";
    }
    if ( degree > 0 ) {
        std::cout << " (remaining factor of degree " << degree << ")";
    }
    std::cout << std::endl;
    delete[] current;
}

void print(
    double       const p[],
    unsigned int const deg_p
) {
    assert( p != nullptr );
    bool first{ true };
    for ( unsigned int k{ deg_p + 1 }; k > 0; --k ) {
        unsigned int const exponent{ k - 1 };
        double const coefficient{ p[exponent] };
        if ( coefficient == 0.0 ) {
            continue;
        }
        if ( !first ) {
            std::cout << (coefficient < 0.0 ? " - " : " + ");
        } else if ( coefficient < 0.0 ) {
            std::cout << "-";
        }
        std::cout << std::abs( coefficient );
        if ( exponent > 0 ) {
            std::cout << "x";
            if ( exponent > 1 ) {
                std::cout << "^" << exponent;
            }
        }
        first = false;
    }
    if ( first ) {
        std::cout << "0";
    }
    std::cout << std::endl;
}