/* vim: set sts=4 sw=4 noet : */

/*
 * Copyright (c) 2007, Fernando J. Pereda <ferdy@gentoo.org>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 *     * Redistributions of source code must retain the above copyright notice,
 *       this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of the program nor the names of its
 *       contributors may be used to endorse or promote products derived from
 *       this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * Compile with:
 *
 *     gcc -o pi_mandel_gmp pi_mandel_gmp.c -lgmp
 *
 * And run with:
 *
 *     ./pi_mandel_gmp
 *
 * Output is something like:
 *
 *  a      || k         || aprox. pi
 * --------++-----------++-----------
 *  1      || 3         || 3.0000000
 *  0.1    || 33        || 3.3000000
 *  0.01   || 315       || 3.1500000
 *  0.001  || 3143      || 3.1430000
 *  0.0001 || 31417     || 3.1417000
 *  1e-05  || 314160    || 3.1416000
 *  1e-06  || 3141593   || 3.1415930
 *  1e-07  || 31415927  || 3.1415927
 */

#include <stdlib.h>
#include <stdio.h>
#include <gmp.h>

static unsigned mandelbrot_it(mpf_t cX, mpf_t cY)
{
	mpf_t x;
	mpf_t y;

	mpf_init_set(x, cX);
	mpf_init_set(y, cY);

	mpf_t x2;
	mpf_t y2;

	mpf_init_set(x2, x);
	mpf_mul(x2, x, x);
	mpf_init_set(y2, y);
	mpf_mul(y2, y, y);

	mpf_t x2y2;
	mpf_init_set(x2y2, x2);
	mpf_add(x2y2, x2, y2);

	unsigned it = 0;
	for (;;) {
		it++;
		if (mpf_cmp_ui(x2y2, 4) > 0)
			break;
		mpf_mul_ui(y, y, 2);
		mpf_mul(y, y, x);
		mpf_add(y, y, cY);

		mpf_sub(x, x2, y2);
		mpf_add(x, x, cX);

		mpf_mul(x2, x, x);
		mpf_mul(y2, y, y);

		mpf_add(x2y2, x2, y2);
	}

	mpf_clear(x);
	mpf_clear(y);
	mpf_clear(x2);
	mpf_clear(y2);
	mpf_clear(x2y2);

	return it;
}

int main(int argc, char *argv[])
{

	mpf_t cX;
	mpf_init_set_d(cX, -0.75L);

	mpf_t cY;

	printf(" %-6s || %-9s || %s\n", "a", "k", "aprox. pi");
	printf("--------++-----------++-----------\n");

	double a;
	for (a = 1L; a > 0.0000001L; a /= 10L) {
		mpf_init_set_d(cY, a);

		unsigned k = mandelbrot_it(cX, cY);
		mpf_mul_ui(cY, cY, k);
		gmp_printf(" %-6g || %-9d || %.7Ff\n", a, k, cY);

		mpf_clear(cY);
	}

	mpf_clear(cX);

	return EXIT_SUCCESS;
}