Initial

1 files changed, 420 insertions, 0 deletions

diff --git a/roots.cpp b/roots.cpp
new file mode 100644
index 0000000..617eea3
--- /dev/null
+++ b/roots.cpp
@@ -0,0 +1,420 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <tuple>
+#include <functional>
+#include <cstring>
+#include <cmath>
+#include <cstdint>
+#include <ctime>
+#include <cassert>
+#include "lodepng.h"
+
+using namespace std;
+
+using i64 = int64_t;
+
+const double EPSILON=1e-8;
+
+
+i64 ipow(i64 b,i64 e){
+	assert(e>=0);
+	if(e==0)return 1;
+	if(e==1)return b;
+	i64 x=1;
+	while(e!=0){
+		if(e&1)x*=b;
+		b*=b;
+		e>>=1;
+	}
+	return x;
+}
+
+int ilog(i64 n){
+	assert(n!=0);
+	int l=0;
+	while(n>>=1)l++;
+	return l;
+}
+
+
+class Polynomial;
+ostream& operator<<(ostream &os,const Polynomial &poly);
+
+class Polynomial : public vector<double> {
+public:
+	Polynomial() = default;
+	Polynomial(size_t sz):vector<double>(sz){}
+	Polynomial(initializer_list<double> init):vector<double>(init){}
+
+	Polynomial d() const {
+		Polynomial res(size()-1);
+		for(size_t i=1;i<size();i++){
+			res[i-1]=operator[](i)*i;
+		}
+		return res;
+	}
+
+	double ev(double x) const {
+		if(size()==0)return 0;
+		if(size()==1)return operator[](0);
+		assert(operator[](size()-1)!=0);
+		double res=operator[](0);
+		double xn=x;
+		for(size_t i=1;i<size();i++){
+			res+=operator[](i)*xn;
+			xn*=x;
+		}
+		return res;
+	}
+
+	bool odd() const {
+		return size()%2==0;
+	}
+	bool even() const {
+		return size()%2==1;
+	}
+
+	int degree() const {
+		assert(size()>0);
+		return size()-1;
+	}
+
+	// Returns quotient, *this becomes remainder
+	Polynomial div(const Polynomial &p){
+		// cerr<<*this<<" / "<<p<<" = ";
+		assert(p.size()>0);
+		Polynomial quo(size()-p.size()+1);
+		for(ssize_t i=size()-p.size();i>=0;i--){
+			double q=operator[](i+p.size()-1)/p.back();
+			quo[i]=q;
+			for(size_t j=0;j<p.size();j++){
+				operator[](i+j)-=q*p[j];
+			}
+		}
+		while(size()>1&&abs(back())<EPSILON)pop_back();
+		// cerr<<quo<<"  r. "<<*this<<endl;
+		return quo;
+	}
+};
+
+ostream& operator<<(ostream &os,const Polynomial &poly){
+	if(poly.degree()==0)return os<<poly[0];
+	bool first=true;
+	for(size_t i=0;i<poly.size();i++){
+		if(poly[i]==0)continue;
+		if(first)first=false;
+		else os<<" + ";
+		if(i==0)os<<poly[i];
+		else {
+			if(poly[i]==1)os<<"X";
+			else if(poly[i]==-1)os<<"-X";
+			else os<<poly[i]<<" X";
+			if(i!=1)os<<'^'<<i;
+		}
+	}
+	return os;
+}
+
+struct Com{
+	double r,i;
+};
+
+ostream& operator<<(ostream &os,const Com &com){
+	return os<<'('<<com.r<<','<<com.i<<')';
+}
+
+namespace FR {
+	double bisect(const Polynomial &poly,int niters){
+		assert(poly.odd());
+		double L=-10,R=10;
+		if(poly.ev(L)==0)return L;
+		if(poly.ev(R)==0)return R;
+		while((poly.ev(L)>0)==(poly.ev(R)>0)){
+			L-=20;
+			R+=20;
+		}
+		double M=(L+R)/2;
+		for(int i=0;i<niters;i++){
+			M=(L+R)/2;
+			if(poly.ev(M)==0)break;
+			if((poly.ev(M)>0)==(poly.ev(L)>0))L=M;
+			else R=M;
+		}
+		return M;
+	}
+
+	double newton(const Polynomial &poly){
+		assert(poly.odd());
+		Polynomial deriv=poly.d();
+		double x=0,xn;
+		do {
+			double d=deriv.ev(x);
+			while(d==0){
+				x+=0.01;
+				d=deriv.ev(x);
+			}
+			xn=x-poly.ev(x)/d;
+		} while(xn-x>EPSILON);
+		return xn;
+	}
+
+	pair<double,double> bairstow(const Polynomial &poly){
+		assert(poly.even()&&poly.degree()>=2);
+		const size_t n=poly.degree();
+		double u=poly[n-1]/poly[n],v=poly[n-2]/poly[n];
+		if(abs(u)<EPSILON)u=1;
+		if(abs(v)<EPSILON)v=1;
+		for(int iter=0;iter<30;iter++){
+			double barr[n+1],farr[n+1];
+			barr[n]=barr[n-1]=0;
+			farr[n]=farr[n-1]=0;
+			for(ssize_t i=n-2;i>=0;i--){
+				barr[i]=poly[i+2]-u*barr[i+1]-v*barr[i+2];
+				farr[i]=barr[i+2]-u*farr[i+1]-v*farr[i+2];
+			}
+			double c=poly[1]-u*barr[0]-v*barr[1];
+			double d=poly[0]-v*barr[0];
+			double g=barr[1]-u*farr[0]-v*farr[1];
+			double h=barr[0]-v*farr[1];
+
+			double discr=1/(v*g*g+h*(h-u*g));
+			// cerr<<"u="<<u<<" v="<<v<<" g="<<g<<" h="<<h<<" c="<<c<<" d="<<d<<" discr="<<discr<<endl;
+			if(abs(c)<EPSILON&&abs(d)<EPSILON){
+				break;
+			}
+
+			v-=discr*(-g*v*c+(g*u-h)*d);
+			u-=discr*(-h*c+g*d);
+		}
+		// cerr<<"bairstow("<<poly<<") = {"<<u<<','<<v<<'}'<<endl;
+		return {u,v};
+	}
+
+	pair<Com,Com> quadratic(double u,double v){
+		double D=u*u-4*v;
+		if(D>=0){
+			double sqrtD2=sqrt(D)/2;
+			return {{-u/2+sqrtD2,0},{-u/2-sqrtD2,0}};
+		} else {
+			double sqrtD2=sqrt(-D)/2;
+			return {{-u/2,sqrtD2},{-u/2,-sqrtD2}};
+		}
+	}
+
+	vector<Com> allroots(const Polynomial &poly){
+		assert(poly.degree()>=1);
+		if(poly.degree()==1){
+			return {{-(double)poly[0]/poly[1],0}};
+		}
+		if(poly.odd()){
+			double rt=bisect(poly,30);
+			vector<Com> r=allroots(Polynomial(poly).div({-rt,1}));
+			r.push_back({rt,0});
+			return r;
+		}
+		pair<double,double> quadr=bairstow(poly);
+		vector<Com> r=poly.degree()==2
+			? vector<Com>()
+			: allroots(Polynomial(poly).div({quadr.second,quadr.first,1}));
+		pair<Com,Com> rts=quadratic(quadr.first,quadr.second);
+		r.push_back(rts.first);
+		r.push_back(rts.second);
+		return r;
+	}
+}
+
+void check(const Com &c){
+	assert(!::isnan(c.r)&&!::isnan(c.i));
+}
+void check(const vector<Com> &v){
+	for(const Com &c : v)check(c);
+}
+
+template <Polynomial(*F)(i64 n)>
+class Generator{
+	i64 n=0;
+	bool isEnd=false;
+	Polynomial current;
+
+	Generator(bool isEnd):isEnd(isEnd){}
+
+	static Generator endGen;
+
+public:
+	Generator(){
+		operator++();
+	}
+
+	Generator begin() const {
+		return Generator(*this);
+	}
+
+	const Generator& end(){
+		return endGen;
+	}
+
+	bool operator!=(const Generator &g){
+		return isEnd!=g.isEnd;
+	}
+
+	Generator& operator++(){
+		if(!isEnd){
+			current=F(n++);
+			if(current.size()==0)isEnd=true;
+		}
+		return *this;
+	}
+
+	const Polynomial& operator*(){
+		return current;
+	}
+};
+
+template <Polynomial(*F)(i64)>
+Generator<F> Generator<F>::endGen = Generator<F>(true);
+
+namespace Gens {
+	template <int N,int D>
+	Polynomial christensen(i64 n){
+		int d=1;
+		while(d<=D){
+			i64 p=ipow(2*N+1,d);
+			if(n<p)break;
+			n-=p;
+			d++;
+		}
+		if(d>D)return {};
+		Polynomial poly(d+1);
+		for(int dd=0;dd<d;dd++){
+			poly[dd]=n%(2*N+1)-N;
+			n/=2*N+1;
+		}
+		poly[d]=1;
+		return poly;
+	}
+
+	template <int D>
+	Polynomial derbyshire(i64 n){
+		n+=2;
+		const int d=ilog(n);
+		if(d>D)return {};
+		Polynomial poly(d+1);
+		for(int dd=0;dd<d;dd++){
+			poly[dd]=n&(1<<dd)?-1:1;
+		}
+		poly[d]=1;
+		return poly;
+	}
+}
+
+
+void writeData(const char *fname,int *tally,int width,int height){
+	ofstream f(fname);
+	assert(f);
+	f<<width<<' '<<height<<'\n';
+	for(int i=0;i<width*height;i++){
+		f<<tally[i]<<'\n';
+	}
+}
+
+tuple<int*,int,int> readData(const char *fname){
+	ifstream f(fname);
+	assert(f);
+	int width,height;
+	f>>width>>height;
+	assert(width>0&&width<10000&&height>0&&height<10000);
+	int *tally=new int[width*height];
+	for(int i=0;i<width*height;i++){
+		f>>tally[i];
+	}
+	return make_tuple(tally,width,height);
+}
+
+
+inline double interpCurve(double x){
+	// return atan(30*x);
+	return pow(x,1.0/2);
+	// return sqrt(log(18*pow(x,0.7)+1));
+	// return x;
+}
+
+inline double interp(double x,double clip){
+	return min(interpCurve(x)/interpCurve(clip),1.0);
+}
+
+void buildImage(unsigned char *image,const int *tally,int width,int height){
+	int maxval=0;
+	for(int i=0;i<width*height;i++){
+		maxval=max(maxval,(int)tally[i]);
+	}
+	for(int i=0;i<width*height;i++){
+		// image[3*i+0]=interp((double)min(tally[i],300)/300,0.4)*255;
+		image[3*i+0]=interp((double)min(tally[i],40)/40,0.4)*255;
+		// image[3*i+0]=interp((double)tally[i],0.4)*255;
+		image[3*i+1]=0;
+		image[3*i+2]=0;
+	}
+}
+
+void writePlot(int *tally,int width,int height){
+	unsigned char *image=new unsigned char[width*height*3];
+	buildImage(image,tally,width,height);
+	int ret=lodepng::encode("out.png",image,width,height,LCT_RGB);
+	if(ret){
+		cerr<<lodepng_error_text(ret)<<endl;
+		exit(1);
+	}
+	delete[] image;
+}
+
+
+int main(int argc,char **argv){
+
+	int *tally;
+	int width,height;
+
+	if(argc==1){
+		// const double minx=-3,maxx=3,miny=-3,maxy=3;
+		const double minx=0,maxx=1.5,miny=0,maxy=1.5;
+		width=2000;
+		height=width/(maxx-minx)*(maxy-miny);
+
+		tally=new int[width*height]();
+
+		int prevdegree=0;
+		clock_t clk=clock();
+
+		// for(const Polynomial &poly : Generator<Gens::christensen<4,5>>()){
+		for(const Polynomial &poly : Generator<Gens::derbyshire<24>>()){
+			// cerr<<"poly = "<<poly<<endl;
+			int degree=poly.degree();
+			if(degree!=prevdegree){
+				if(prevdegree!=0){
+					clock_t c2=clock();
+					cerr<<'('<<(double)(c2-clk)/CLOCKS_PER_SEC<<"sec)"<<endl;
+					clk=c2;
+				}
+				prevdegree=degree;
+				cerr<<"d="<<degree<<' ';
+			}
+			vector<Com> r=FR::allroots(poly);
+			for(const Com &c : r){
+				int x=(c.r-minx)/(maxx-minx)*width;
+				int y=(c.i-miny)/(maxy-miny)*height;
+				if(x<0||x>=width||y<0||y>=height)continue;
+				tally[width*y+x]++;
+			}
+		}
+
+		writeData("out.txt",tally,width,height);
+	} else if(argc==2){
+		tie(tally,width,height)=readData(argv[1]);
+	} else {
+		cerr<<"Usage: "<<argv[0]<<" [out.txt]"<<endl;
+		return 1;
+	}
+
+	writePlot(tally,width,height);
+	delete[] tally;
+}