#include <Base/Common/BIASpragma.hh>
#include <MathAlgo/PolynomialSolve.hh>
#include <Base/Math/Random.hh>
#include <algorithm>
using namespace BIAS;
using namespace std;
bool debug=false;
double maxerr=1e-8;
namespace std{
bool operator<(complex<double>& left, complex<double>& right)
{
return (left.real()<right.real());
}
bool operator<(const complex<double>& left, const complex<double>& right)
{
return (left.real()<right.real());
}
}
void solout(vector<POLYNOMIALSOLVE_TYPE>& sol)
{
sort(sol.begin(), sol.end());
int s=(int)sol.size();
cout.precision(20);
for (int i=0; i<s; i++)
if (debug) cout << sol[i] << " ";
if (debug) cout << endl;
}
void solout(vector<complex<POLYNOMIALSOLVE_TYPE> >& sol)
{
sort(sol.begin(), sol.end());
int s=(int)sol.size();
cout.precision(20);
for (int i=0; i<s; i++)
if (debug) cout << sol[i] << " ";
if (debug) cout << endl;
}
void errout(vector<POLYNOMIALSOLVE_TYPE>& sol,
vector<POLYNOMIALSOLVE_TYPE>& gtsol)
{
sort(sol.begin(), sol.end());
int s=(int)sol.size();
cout.precision(20);
if (debug) cout << "error: ";
for (int i=0; i<s; i++){
if (debug) cout << fabs(sol[i]-gtsol[i]) << " ";
if (fabs(sol[i]-gtsol[i])>maxerr) {
cerr << "\n\n\nerror bigger than "<<maxerr<<"\n" <<endl;
debug=true;
cout << "ground truth: ";
solout(gtsol);
cout << "computed: ";
solout(sol);
cout << "error: ";
for (int i=0; i<s; i++) cout << fabs(sol[i]-gtsol[i]) << " ";
cout << endl;
cerr << "\n\n\nthis happens if the roots of the polynomial are close together\n\n\n"
<<endl;
BIASASSERT(false);
}
}
if (debug) cout << endl << endl;
}
void errout(vector<complex<POLYNOMIALSOLVE_TYPE> >& sol,
vector<POLYNOMIALSOLVE_TYPE>& gtsol)
{
sort(sol.begin(), sol.end());
int s=(int)sol.size();
cout.precision(20);
if (debug) cout << "error: ";
for (int i=0; i<s; i++){
if (debug) cout << fabs((sol[i]).real()-gtsol[i]) << " ";
if (fabs((sol[i]).real()-gtsol[i])>maxerr) {
cerr << "\n\n\nerror bigger than "<<maxerr<<"\n" <<endl;
debug=true;
cout << "ground truth: ";
solout(gtsol);
cout << "computed: ";
solout(sol);
cout << "error: ";
for (int i=0; i<s; i++) cout << fabs((sol[i]).real()-gtsol[i]) << " ";
cout << endl;
cerr << "\n\n\nthis happens if the roots of the polynomial are close together\n\n\n"
<<endl;
BIASASSERT(false);
}
}
if (debug) cout << endl << endl;
}
int main()
{
int num=9000;
int maxdeg=8;
double maxsol=5e1;
int deg, nums;
vector<POLYNOMIALSOLVE_TYPE> coeff, sol, gtsol;
vector<complex<POLYNOMIALSOLVE_TYPE> > csol;
for (int n=0; n<num; n++){
if (debug){
cout << "-----------------"<< setw(3) << n <<"------------------"<<endl;
cout << "using a polynom of degree "<<deg<<endl;
} else {
if (n%10==0){
cout << n << flush;
} else {
cout << "." << flush;
}
}
gtsol.resize(deg);
for (int d=0; d<deg; d++){
if (debug)
cout << d << "th root: "<<gtsol[d]<<endl;
}
sort(gtsol.begin(), gtsol.end());
if (debug){
for (int i=0; i<deg+1; i++){
cout << i << "th coefficient: "<< coeff[i]<<endl;
}
}
if (debug) {
cout << "polynomial at ground truth roots: "<<endl;
for (int i=0; i<deg; i++){
<< " ";
cout << endl;
}
}
nums=ps.
Solve(coeff, sol);
if (debug) cout <<"Solve (real) found "<<nums<<" solutions:\n";
if (nums>0) { solout(sol); errout(sol, gtsol); }
nums=ps.
Solve(coeff, csol);
if (debug) cout <<"Solve (complex) found "<<nums<<" solutions:\n";
if (nums>0) { solout(csol); errout(csol, gtsol); }
if (debug)cout<<"Numeric (real) found "<<nums<<" solutions:\n";
if (nums>0) { solout(sol); errout(sol, gtsol); }
if (debug) cout<<"Numeric (complex) found "<<nums<<" solutions:\n";
if (nums>0) { solout(csol); errout(csol, gtsol); }
if (deg<=4){
if (debug) cout << "Analytic (real) found "<<nums<<" solutions:\n";
if (nums>0) { solout(sol); errout(sol, gtsol); }
}
if (deg<=3){
if (debug) cout << "Analytic (complex) found "<<nums<<" solutions:\n";
if (nums>0) { solout(csol); errout(csol, gtsol); }
}
}
return 0;
}