/*
    polynom_test.cpp version 0.11               Pieter Suurmond, jan 2 18, 2020.
    Test our 'polynomial library'.
*/

#include <iostream>
#include <iomanip>
#include <cstdlib>          // For random().
#include <complex>
#include <vector>

using namespace std;

#include "polynom.hpp"      // Needs namespace std, <complex> and <vector>.


static double rr() // Floating point random between -10 (incl.) and +10 (excl.).
{
    return (double)((random() & 32767) - 16384) / 1638.40;
}


static long test(int degree, long num_polynomials)
{
    vector<complex<double> >    coeffs(degree + 1);         // INPUT.
    vector<complex<double> >    roots(degree);              // OUT/INPUT.
    vector<complex<double> >    coeffs_back(degree + 1);    // OUTPUT.
    
    long failed_laguerre  = 0L; // Count number of Laguerre failures.
    long failed_coeffdiff = 0L; // Count number of coefficient differences.

    cout << "Factoring and re-expanding " << num_polynomials
         << " polynomials of degree " << degree << "\n";
    for (long test = 0; test < num_polynomials; test++)
        {
        int n, e;
        //----------------------------- INVENT SOME RANDOM COEFFICIENTS:
        for (n = 0; n <= degree; n++)
            coeffs[n] = complex<double>(rr(), rr());
        //----------------------------- FIND COMPLEX POLYNOMIAL ROOTS:
        e = factor(coeffs,           // Input.
                   roots);           // Output.
        if (e)
            failed_laguerre++;       // Or synth.div.remainder unacceptable.
        else
            {
            //------------------------- CONVERT BACK TO COEFFICIENTS AND TEST:
            e = expand(roots,        // Input.
                       coeffs_back); // Output.
            if (e)
                { cout << "ERROR: expand()=" << e << "!\n"; break; }
            // Scale coeffs. Coeff of highest degree cannot be zero.
            int             n_roots = roots.size();
            complex<double> f = coeffs[n_roots] / coeffs_back[n_roots];
            for (n = 0; n <= n_roots; n++)
                coeffs_back[n] *= f;
            // Finally, compare with the original coefficients:
            for (n = 0; (!e) && (n <= n_roots); n++)
                {
                double d = abs(coeffs_back[n] - coeffs[n]);
                if (d > 1e-4)           // Differences above 0.0001 are counted.
                    failed_coeffdiff++;
                }
            }
        }
    if (failed_laguerre)
        cout << " " << failed_laguerre << " laguerre failed to converge.\n";
    if (failed_coeffdiff)
        cout << " " << failed_coeffdiff << " had differing coefficients.\n";
    long failed = failed_coeffdiff + failed_laguerre;
    cout << " " << 100.0 * (num_polynomials - failed)
                         / double(num_polynomials) << "% successful.\n";
    return failed;
}


int main()
{
    cout << "\nTesting 'POLYNOM LIBRARY' with a million complex polynomials.\n";
    cout << setiosflags(ios::fixed) << setprecision(3);

    try {
        long failed = 0L;

        for (int degree = 1; degree <= 40; degree++)
            failed += test(degree, 25000);

        cout << "Done. In total, " << failed << " polynomial-conversions failed.\n";
        }
    catch (const std::bad_alloc& exception)
        {
        cout << "Out of memory!\n";
        }
    
    return 0;
}