fractal/fun-gen.c

302 lines
5.7 KiB
C
Raw Normal View History

2020-03-07 23:32:57 +01:00
/*
* 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
2020-03-07 23:32:57 +01:00
#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;
}