#include "Base/Debug/Error.hh"
#include <MathAlgo/Laguerre.hh>
#include <Base/Common/CompareFloatingPoint.hh>
#include <iostream>
#include <algorithm>
#ifdef WIN32
#include <functional>
#endif
using namespace BIAS;
using namespace std;
double frand(double min = -10.0 , double max = 10.0)
{
double tmin = min;
double tmax = max;
const unsigned int res = 10000000;
unsigned int i=rand()%res;
double scale = ((double) i) / ((double) res);
double diff = tmax - tmin;
return(tmin+scale*diff);
}
unsigned int GetNumberOfBinaryOnes(unsigned int val,
unsigned int mask)
{
unsigned int bitmask = mask;
unsigned int temp = val;
unsigned int count = 0;
while((bitmask)!=0)
{
if((temp&1)==1)
{
count++;
}
temp = temp>>1;
bitmask = bitmask>>1;
}
return count;
}
void GetIndexListFromBinaryZeros(unsigned int val,
unsigned int mask,
std::vector<unsigned int>& indices)
{
indices.clear();
unsigned int bitmask = mask;
unsigned int temp = val;
unsigned int count = 0;
while((bitmask)!=0)
{
if((temp&1)==0)
{
indices.push_back(count);
}
temp = temp>>1;
bitmask = bitmask>>1;
count++;
}
}
void CreateRandomPolynomial(unsigned int degree,
vector< complex<double> >& coeffs,
vector< complex<double> >& roots)
{
BIASASSERT(degree>0);
roots.clear();
coeffs.clear();
std::vector<unsigned int> tempIDs;
unsigned int coeffidx ;
for(unsigned int r=0; r<degree; r++)
{
roots.push_back(frand());
}
coeffs.resize(degree+1,complex<double>(0,0));
unsigned int maxiter = 1<<(degree);
unsigned int bitmask = maxiter - 1;
for(unsigned int c=0; c<maxiter; c++)
{
GetIndexListFromBinaryZeros(c,bitmask,tempIDs);
coeffidx = GetNumberOfBinaryOnes(c,bitmask);
complex<double> someroots(1,0);
for(unsigned int r=0;r<tempIDs.size();r++)
{
someroots = someroots * (-roots[tempIDs[r]]);
}
BIASASSERT(coeffidx<coeffs.size());
coeffs[coeffidx]= coeffs[coeffidx] + someroots;
}
}
struct complex_comp :
public binary_function< complex<double>, complex<double>, bool>
{
bool operator()(complex<double> x, complex< double> y)
{
if(x.real()<=y.real())
{
if(x.real()==y.real())
{
return(x.imag()<y.imag());
}
return(true);
}
return(false);
}
};
double MeanError(vector< complex<double> >& v1,
vector< complex<double> >& v2,
double& min, double& max)
{
BIASASSERT(v1.size()==v2.size())
sort(v1.begin(),v1.end(),complex_comp());
sort(v2.begin(),v2.end(),complex_comp());
double error =0 ;
double tmax=0, tmin=10000000;
for(unsigned int roo=0;roo<v1.size(); roo++)
{
double n = norm(v2[roo]-v1[roo]);
if(n>tmax)tmax= n;
if(n<tmin)tmin= n;
error+= n;
}
max=tmax;
min=tmin;
return (error/v1.size());
}
bool Equal(
const complex<double>& a,
const complex<double>& b)
{
return((
Equal(a.real(),b.real(),10e-12))&&
((
Equal(a.imag(),b.imag(),10e-12))));
}
int main(int argc, char *argv[])
{
vector< complex<double> > realcoeff;
vector< complex<double> > expectedroots;
vector< complex<double> > roots;
unsigned int deg = 8;
double max = -10e6;
double min = 10e6;
double meanmeanerror=0;
double meanmin = 10000000;
double meanmax = 0;
double tmin;
double tmax;
for(unsigned run=0; run<100; run++)
{
CreateRandomPolynomial(deg,realcoeff,expectedroots);
l.
Solve(realcoeff, roots);
double err;
if(run==0)
err=MeanError(expectedroots,roots,min,max);
else
{
err=MeanError(expectedroots,roots,tmin,tmax);
if(min>tmin)min=tmin;
if(max<tmax)max=tmax;
}
cout<<err<<" min: "<< min<<" max: "<< max<<endl;
if(err<meanmin)meanmin=err;
if(err>meanmax)meanmax=err;
meanmeanerror +=err;
}
meanmeanerror /= 100.0;
cout<<"-----------------------------------------------------------"<<endl;
cout<<"min error (no mean):"<< min <<endl;
cout<<"max error (no mean):"<< max <<endl;
cout<<"mean min error:"<< meanmin <<endl;
cout<<"mean max error:"<< meanmax <<endl;
cout<<"mean mean error:"<< meanmeanerror <<endl;
return 0;
}