aberth.cc

#include <complex>
#include <cstdlib>
#include <iostream>
#include <fstream>

typedef double number_t;

typedef std::complex<number_t> complex_t;
void evaluate_derivative(complex_t& y, complex_t& yp, const complex_t* ap, const int n, const complex_t& x) {
	y = ap[n - 1];
	yp = 0;

	for(int i = n - 2; i >= 0; i--) {
		y = y * x + ap[i];
		yp = yp * x + ap[i + 1] * complex_t(i + 1, 0);
	}
}


complex_t synthetic_division(complex_t* bp, const complex_t* ap, const int n, const complex_t& x) {
	complex_t acc = ap[n - 1];

	for(int i = n - 2; i >= 0; i--) {
		bp[i] = acc;
		acc = acc * x + ap[i];
	}

	return acc;
}

int aberth(complex_t* bp, const complex_t* app, const int nn) {
	int n = nn;
	const int max_iterations = 10;
	const int max_retries = 10;
	const number_t TOL = 2.22e-16;
	int iterations = 0;
	int retries = 0;
	complex_t wp[n];
	complex_t rp[n];
	int r_count = 0;
	complex_t ap[n];

	number_t tol = 0;

	for(int i = 0; i <= n; i++) {
		ap[i] = app[i];
		tol += norm(ap[i]);
	}

	tol *= TOL / n;
	tol = TOL;

	for(int i = 0; i < n; i++) {
		bp[i] = std::polar(1.0, M_PI * i * 2 / n);
	}

	while(r_count < nn && retries < max_retries) {
		/*if(n==2) {
			complex_t disc_part = std::sqrt(ap[1] * ap[1] - complex_t(4,0)*ap[0] * ap[2])/(complex_t(2,0)*ap[2]);
			complex_t b_part = -ap[1] / (complex_t(2,0)*ap[2]);
			rp[r_count] = b_part + disc_part;
			rp[r_count + 1] = b_part - disc_part;
			r_count+=2;
			break;
		}*/



		for(int k = 0; k < n; k++) {
			complex_t zk;
			complex_t zpk;

			evaluate_derivative(zk, zpk, ap, n + 1, bp[k]);

			if(norm(zpk) == 0) {
				continue;
			}

			complex_t t = zk / zpk;

			complex_t s = 0;

			for(int j = 0; j < n; j++) {
				if(j != k) {
					s += complex_t(1, 0) / (bp[k] - bp[j]);
				}
			}

			complex_t u = complex_t(1, 0) - t * s;

			if(norm(u) == 0) {
				wp[k] = 0;
				continue;
			}

			wp[k] = -t / u;
		}

		for(int k = 0; k < n; k++) {
			bp[k] += wp[k];

			/*complex_t zk;
			complex_t zpk;

			evaluate_derivative(zk, zpk, ap, n+1, bp[k]);
			std::cerr << r_count << " :: " << abs(zk) << std::endl;*/
		}

		complex_t cp[n + 1];

		for(int k = 0; k < n; k++) {

			complex_t y = synthetic_division(cp, ap, n + 1, bp[k]);

			if(norm(y) < tol) {


				//std::cerr << tol << std::endl;
				rp[r_count] = bp[k];
				r_count++;
				/*for(int j=0;j<r_count;j++) {
				complex_t a,b;
				evaluate_derivative(a,b,app, nn+1, rp[j]);
				std::cerr << abs(a) << std::endl;
				}*/

				for(int j = k + 1; j < n; j++) {
					bp[j - 1] = bp[j];
				}

				k--;

				n--;

				for(int j = 0; j <= n; j++) {
					ap[j] = cp[j];
				}

				tol = 0;

				for(int i = 0; i <= n; i++) {
					tol += norm(ap[i]);
				}

				tol = TOL / n;
				tol = TOL;


			}
		}

		iterations++;

		if(iterations >= max_iterations) {
			for(int i = 0; i < n; i++) {
				bp[i] = complex_t(rand(), rand()) * complex_t(2.0 / RAND_MAX, 0) - complex_t(1, 1);
			}

			iterations = 0;
			retries++;
		}
	}

	for(int i = 0; i < r_count; i++) {
		/*complex_t zk, zpk;
		do {
		    evaluate_derivative(zk,zpk, app, n+1, rp[i]);
		std::cerr <<i << "/"<<r_count << " "<< abs(zk) << std::endl;
		    if(norm(zpk)>0) {
		        rp[i] -= zk/zpk;
		    }
		} while(abs(zk) > tol);*/

		bp[i] = rp[i];
	}

	return r_count;
}

int partition(int* ap, const int n) {
	{
		int index = rand() % n;
		int t = ap[n - 1];
		ap[n - 1] = ap[index];
		ap[index] = t;
	}
	int index = 0;

	for(int i = 0; i < n - 1; i++) {
		if(ap[i] <= ap[n - 1]) {
			int t = ap[i];
			ap[i] = ap[index];
			ap[index] = t;
			index++;
		}
	}

	{
		int t = ap[n - 1];
		ap[n - 1] = ap[index];
		ap[index] = t;
	}

	return index;
}
void quicksort(int* ap, const int n) {
	if(n <= 1) {
		return;
	}

	int pivot = partition(ap, n);
	//std::cerr << pivot << std::endl;
	int pivot0 = pivot;

	while(pivot0 > 0 && ap[pivot0] == ap[pivot]) {
		pivot0--;
	}

	quicksort(ap, pivot0);
	quicksort(&ap[pivot + 1], n - pivot - 1);
}


int binary_search(const int* ap, const int n, const int k) {
	int high = n - 1;
	int low = 0;

	while(low < high) {
		int mid = low + (high - low) / 2;

		int y = ap[mid];

		if(y < k) {
			low = mid + 1;
		} else if(k < y) {
			high = mid - 1;
		} else {
			return mid;
		}
	}

	return high;
}
int main() {
	const int n = 2;
	const int m = 24;
	complex_t cc[n];


	cc[0] = complex_t(-1, 0);
	cc[1] = complex_t(1, 0);
	/*for(int i=0; i<n; i++) {
	    //cc[i] = std::polar(1.0,static_cast<number_t>(i)*M_PI*2/n);
	    cc[i] = complex_t(pow(2,i));
	}*/


	number_t w = 2;

	const int screen_size = 1 << 12;
	int* sp;

	sp = new int[screen_size * screen_size];

	for(int i = 0; i < screen_size * screen_size; i++) {
		sp[i] = 0;
	}

	int max_iterations = 2 << m;

	for(;;) {

		#pragma omp parallel for

		for(int iterations = 0; iterations < max_iterations; iterations++) {
			complex_t bp[m];
			complex_t cp[m + 1];
			int idx[m + 1];

#if 1

			for(int i = 0; i <= m; i++) {
				idx[i] = (iterations >> i) & 1;
			}

#else

			for(int i = 0; i < m; i++) {
				idx[i] = rand() % n;
			}

			if(cc[0] != complex_t(0)) {
				idx[0] = rand() % n;
			} else {
				idx[0] = (rand() % (n - 1)) + 1;
			}

			if(cc[0] != complex_t(0)) {
				idx[m] = rand() % n;
			} else {
				idx[m] = (rand() % (n - 1)) + 1;
			}

#endif

			bool reverse_larger = false;

			for(int i = m; i >= 0; i--) {
				if(idx[i] < idx[m - i]) {
					reverse_larger = true;
					break;
				} else if(idx[i] > idx[m - 1]) {
					break;
				}
			}

			if(reverse_larger) {
				continue;
			}

			for(int i = 0; i <= m; i++) {
				cp[i] = cc[idx[i]];
			}

			for(int i = 0; i <= m; i++) {
				cp[i] /= cp[m];
			}

			int bn = aberth(bp, cp, m);
			/*for(int i=0;i<bn;i++) {
				complex_t y, yp;

				evaluate_derivative(y,yp, cp, m+1, bp[i]);


				std::cerr <<i << "/" << m << " : " << bp[i] << " " <<abs(y) <<std::endl;
			}*/



			for(int i = 0; i < bn * 2; i++) {
				complex_t t = bp[i];

				if(i >= bn) {
					t = complex_t(1, 0) / t;
				}

				complex_t z = (t + complex_t(w, w)) * complex_t(screen_size / (2 * w), 0);

				int x, y;
				x = z.real();
				y = z.imag();

				//std::cerr << x << "," << y << std::endl;
				if(0 <= x && x < screen_size && 0 <= y && y < screen_size) {
					#pragma omp atomic
					sp[y * screen_size + x]++;
				}
			}

			if(iterations & 0xff == 0) {
				std::cerr << iterations << std::endl;
			}
		}

		{
			std::ofstream out;
			int* tp = new int[screen_size * screen_size];

			for(int i = 0; i < screen_size * screen_size; i++) {
				tp[i] = sp[i];
			}

			quicksort(tp, screen_size * screen_size);

			int uniques = 1;

			for(int i = 1; i < screen_size * screen_size; i++) {
				if(tp[uniques - 1] != tp[i]) {
					tp[uniques] = tp[i];
					uniques++;
				}
			}

			out.open("out.pnm");
			out << "P2" << std::endl;
			out << screen_size << " " << screen_size << std::endl;
			int max = 0;
			number_t avg = 0;

			for(int i = 0; i < screen_size * screen_size; i++) {
				if(max < sp[i]) {
					max = sp[i];
					avg += static_cast<number_t>(sp[i]) / (screen_size * screen_size);
				}
			}

			std::cerr << "max: " << max << std::endl;
			std::cerr << "uniques: " << uniques << std::endl;
			const int graymax = 65535;
			out << graymax << std::endl;

			for(int i = 0; i < screen_size * screen_size; i++) {
				//int c = graymax.0*pow(static_cast<number_t>(sp[i])/avg,1.0/2.0);
				int c = graymax * (static_cast<number_t>(binary_search(tp, uniques, sp[i])) / (uniques - 1));
				c = c < graymax ? c : graymax;
				out << c << std::endl;
				//std::cerr << c << std::endl;
			}

			out.close();
			return 0;
		}
	}

	delete [] sp;
}