diff --git a/criticals.c b/criticals.c
index 877b0cf..11233b7 100644
--- a/criticals.c
+++ b/criticals.c
@@ -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;
diff --git a/fun-gen.c b/fun-gen.c
index 77bad49..ddd5260 100644
--- a/fun-gen.c
+++ b/fun-gen.c
@@ -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);