Lucas/fractal
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

[fun-gen] add Halley method

Browse Source
This commit is contained in:
Lucas 2020-03-16 01:54:20 +00:00
parent 5078ef1578
commit 022e2cf2e1
2 changed files with 82 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
criticals.c
View File

@ -48,7 +48,7 @@ criticals_add(long double complex z)
if (idx >= MAX_CRITICALS)
return;
for (i = 0; i < idx; i++)
if (cabsl(z - criticals[i]) < 1.0e-9)
if (cabsl(z - criticals[i]) < 1.0e-8)
return;
criticals[idx++] = z;

87
fun-gen.c
View File

@ -46,13 +46,26 @@ poly_init(struct poly *p, unsigned int degree, const char *expr)
static void
poly_derivate(struct poly p, struct poly *dp, const char *expr)
{
unsigned int i;
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 void
poly_multiply(struct poly p, struct poly q, struct poly *f, const char *expr)
{
int i, j;
poly_init(f, p.degree + q.degree, expr);
for (i = 0; i <= f->degree; i++)
f->coefs[i] = 0.0;
for (i = 0; i <= p.degree; i++)
for (j = 0; j <= q.degree; j++)
f->coefs[j + i] += p.coefs[i] * q.coefs[j];
}
static int
poly_cmp(const void *vp, const void *vq)
{
@ -113,6 +126,17 @@ print_funs(struct poly p, struct poly q)
{
struct poly polys[6];
struct poly dp, ddp, dq, ddq;
int i, j;
for (i = 0; i <= MIN(p.degree, q.degree)
&& p.coefs[i] == 0.0 && q.coefs[i] == 0.0; i++)
;
for (j = i; j <= p.degree; j++)
p.coefs[j - i] = p.coefs[j];
p.degree -= i;
for (j = i; j <= q.degree; j++)
q.coefs[j - i] = q.coefs[j];
q.degree -= i;
polys[0] = p;
polys[1] = q;
@ -169,7 +193,8 @@ usage(void)
const char *p = xgetprogname();
fprintf(stderr, "Usage:\n"
"\t%s p-degree q-degree p-coefs q-coefs\n"
"\t%s -n p-degree p-coefs\n", p, p);
"\t%s -h p-degree p-coefs\n"
"\t%s -n p-degree p-coefs\n", p, p, p);
exit(1);
}
@ -204,6 +229,51 @@ do_poly_div(int argc, char *argv[])
poly_free(&q);
}
static void
do_halley(int argc, char *argv[])
{
struct poly f, df, ddf, p, q, s, t, u;
int i, j, f_degree;
if (argc < 1)
usage();
f_degree = parse_integer(argv[0], 0, MAX_DEGREE);
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, "df");
poly_derivate(df, &ddf, "ddf");
poly_multiply(f, df, &s, "f_df");
poly_multiply(df, df, &t, "df_df");
poly_multiply(f, ddf, &u, "f_ddf");
poly_init(&p, MAX(MAX(t.degree, u.degree) + 1, s.degree), "p");
poly_init(&q, MAX(t.degree, u.degree), "q");
for (i = 0; i <= q.degree; i++)
q.coefs[i] = 2.0 * t.coefs[i] - u.coefs[i];
for (i = 0; i <= p.degree; i++)
p.coefs[i] = (i == 0 ? 0.0 : q.coefs[i - 1])
- 2.0L * (i <= s.degree ? s.coefs[i] : 0.0);
print_funs(p, q);
poly_free(&f);
poly_free(&df);
poly_free(&ddf);
poly_free(&s);
poly_free(&t);
poly_free(&u);
poly_free(&p);
poly_free(&q);
}
static void
do_newton(int argc, char *argv[])
{
@ -237,11 +307,14 @@ do_newton(int argc, char *argv[])
int
main(int argc, char *argv[])
{
int ch, nflag;
int ch, hflag, nflag;
nflag = 0;
while ((ch = getopt(argc, argv, "n")) != -1)
hflag = nflag = 0;
while ((ch = getopt(argc, argv, "hn")) != -1)
switch (ch) {
case 'h':
hflag = 1;
break;
case 'n':
nflag = 1;
break;
@ -251,7 +324,9 @@ main(int argc, char *argv[])
argc -= optind;
argv += optind;
if (nflag)
if (hflag)
do_halley(argc, argv);
else if (nflag)
do_newton(argc, argv);
else
do_poly_div(argc, argv);