fractal/fun-gen.c

/*
 * fractal
 * Written in 2020 by Lucas
 * CC0 1.0 Universal/Public domain - No rights reserved
 *
 * To the extent possible under law, the author(s) have dedicated all
 * copyright and related and neighboring rights to this software to the
 * public domain worldwide. This software is distributed without any
 * warranty. You should have received a copy of the CC0 Public Domain
 * Dedication along with this software. If not, see
 * <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
/* gepopt visibility for GLIBC */
#if !defined(__OpenBSD__)
#define _DEFAULT_SOURCE
#endif

#include <errno.h>
#include <inttypes.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))

static char *argv0;

struct poly {
	unsigned int degree;
	double *coefs;
	const char *expr;
};

static void
poly_init(struct poly *p, unsigned int degree, const char *expr)
{
	if ((p->coefs = malloc((degree + 1) * sizeof(double))) == NULL) {
		fprintf(stderr, "out of memory\n");
		exit(1);
	}

	p->degree = degree;
	p->expr = expr;
}

static void
poly_derivate(struct poly p, struct poly *dp, const char *expr)
{
	unsigned int i;

	poly_init(dp, p.degree == 0 ? 0 : p.degree - 1, expr);
	for (i = p.degree; i > 0; i--)
		dp->coefs[i - 1] = p.coefs[i] * i;
}

static int
poly_cmp(const void *vp, const void *vq)
{
	const struct poly *p = vp, *q = vq;
	return p->degree - q->degree;
}

static void
poly_free(struct poly *p)
{
	free(p->coefs);
}

static void
print_fun_header(const char *fun_name)
{
	printf("static long double complex\n");
	printf("%s(long double complex z)\n", fun_name);
	printf("{\n");
}

static void
print_fun_body(size_t n, struct poly *polys)
{
	size_t i, j;
	unsigned int d;

	printf("\tlong double complex a");
	for (i = 0; i < n; i++)
		printf(", %s = 0.0L", polys[i].expr);
	printf(";\n\n");

	d = 0;
	for (i = 0; i < n; i++) {
		for (; d <= polys[i].degree; d++) {
			if (d == 0)
				printf("\ta = 1.0L;\n");
			else
				printf("\ta *= z;\n");
			for (j = n; j > i; j--)
				if (polys[j - 1].coefs[d] != 0.0L)
					printf("\t%s += %.20f * a;\n",
					    polys[j - 1].expr,
					    polys[j - 1].coefs[d]);
			printf("\n");
		}
	}
}

static void
print_fun_footer(void)
{
	printf("}\n");
}

static void
print_funs(struct poly p, struct poly q)
{
	struct poly polys[6];
	struct poly dp, ddp, dq, ddq;

	polys[0] = p;
	polys[1] = q;
	qsort(polys, 2, sizeof(struct poly), &poly_cmp);

	print_fun_header("f");
	print_fun_body(2, polys);
	printf("\treturn p / q;\n");
	print_fun_footer();

	poly_derivate(p, &dp, "dp");
	poly_derivate(q, &dq, "dq");
	polys[0] = p;
	polys[1] = q;
	polys[2] = dp;
	polys[3] = dq;
	qsort(polys, 4, sizeof(struct poly), &poly_cmp);

	printf("\n");
	print_fun_header("df");
	print_fun_body(4, polys);
	printf("\treturn (dp * q - p * dq) / (q * q);\n");
	print_fun_footer();

	poly_derivate(dp, &ddp, "ddp");
	poly_derivate(dq, &ddq, "ddq");
	polys[0] = p;
	polys[1] = q;
	polys[2] = dp;
	polys[3] = dq;
	polys[4] = ddp;
	polys[5] = ddq;
	qsort(polys, 6, sizeof(struct poly), &poly_cmp);

	printf("\n");
	print_fun_header("ddf");
	print_fun_body(6, polys);
	printf("\treturn ("
	    "q * q * ddp "
	    "- q * (2.0L * dp * dq + p * ddq) "
	    "+ 2.0L * p * dq * dq"
	    ") / (q * q * q);\n");
	/*printf("\treturn (ddp "
	    "- 2.0 * ((dp * q - p * dq) / (q * q)) * dq "
	    "- (p / q) * ddq) / q;\n");*/
	print_fun_footer();

	poly_free(&dp);
	poly_free(&dq);
	poly_free(&ddp);
	poly_free(&ddq);
}

static unsigned int
parse_uint(const char *s)
{
	char *end;
	uintmax_t v;

	errno = 0;
	v = strtoumax(s, &end, 10);
	if (s == end || *end != '\0' || errno != 0 || v > UINT_MAX) {
		fprintf(stderr, "%s: invalid value: %s\n", argv0, s);
		exit(1);
	}

	return (unsigned int)v;
}

static double
parse_double(const char *s)
{
	char *end;
	double d;

	errno = 0;
	d = strtod(s, &end);
	if (s == end || *end != '\0' || errno != 0 || !isfinite(d)) {
		fprintf(stderr, "%s: invalid value: %s\n", argv0, s);
		exit(1);
	}

	return d;
}

static void
usage(void)
{
	fprintf(stderr, "Usage: %s p-degree q-degree p-coefs q-coefs\n"
	    "    %s -n p-degree p-coefs\n", argv0, argv0);
	exit(1);
}

static void
do_poly_div(int argc, char *argv[])
{
	struct poly p, q;
	unsigned int i, j, p_degree, q_degree;

	if (argc < 2)
		usage();

	p_degree = parse_uint(argv[0]);
	q_degree = parse_uint(argv[1]);
	if (argc != p_degree + q_degree + 2 + 2)
		usage();

	poly_init(&p, p_degree, "p");
	poly_init(&q, q_degree, "q");

	for (j = 0, i = 1 + 1 + p.degree + 1 + q.degree;
	    i > 1 + 1 + p.degree;
	    i--, j++)
		q.coefs[j] = parse_double(argv[i]);
	for (j = 0;
	    i > 1;
	    i--, j++)
		p.coefs[j] = parse_double(argv[i]);

	print_funs(p, q);
	poly_free(&p);
	poly_free(&q);
}

static void
do_newton(int argc, char *argv[])
{
	struct poly f, df, p;
	unsigned int i, j, f_degree;

	if (argc < 1)
		usage();

	f_degree = parse_uint(argv[0]);
	if (argc != f_degree + 1 + 1)
		usage();

	poly_init(&f, f_degree, "f");
	for (j = 0, i = 1 + f.degree;
	    i > 0;
	    i--, j++)
		f.coefs[j] = parse_double(argv[i]);

	poly_derivate(f, &df, "q");
	poly_init(&p, f_degree, "p");
	for (i = 0; i <= f.degree; i++)
		p.coefs[i] = (i == 0 ? 0.0 : df.coefs[i - 1]) - f.coefs[i];

	print_funs(p, df);
	poly_free(&f);
	poly_free(&df);
	poly_free(&p);
}

int
main(int argc, char *argv[])
{
	struct poly p, q;
	char *end;
	uintmax_t v;
	unsigned int i, j, p_degree, q_degree;
	int ch, nflag;

	argv0 = argv[0];
	nflag = 0;
	while ((ch = getopt(argc, argv, "n")) != -1)
		switch (ch) {
		case 'n':
			nflag = 1;
			break;
		default:
			usage();
		}
	argc -= optind;
	argv += optind;

	if (nflag)
		do_newton(argc, argv);
	else
		do_poly_div(argc, argv);

	return 0;
}