Writing Heron's formula that is "safer in a giant triangle"

When Implementing Heron's Formula, We Use Kahan's Algorithm "Without Question."

S = \sqrt{s(s-a)(s-b)(s-c)} \quad \left( s = \frac{a+b+c}{2} \right)

This is an all too familiar formula. A monumental achievement discovered by the genius engineer of ancient Greece, Archimedes, and preserved and transmitted by Heron.

And we have "believed" that using this great formula in the form organized without a single flaw by William Kahan leaves no room for exceptions.

But is that really true? We will pull improved and degraded algorithms across different eras, compare their execution results, and contrast their methods.

And allow me to introduce a function written in "my style."

The Most Fundamental Idea

First, when looking at the history of algorithms, I consider who first "properly elucidated" them, regardless of nationality.

For example, the person who proved the existence of 9th-order orthogonal Latin squares and pairs of 10th-order orthogonal Latin squares is Choi Seok-jeong. And the earliest known statement of the Pythagorean theorem is found in the Baudhayana Sulba Sutra.

However, this has nothing to do with cultural relativism or pluralism.

Heron's formula was discovered by Archimedes and preserved by Heron.

So we should pay homage to both while admiring this implementation.

double heron_european(double a, double b, double c) {
    double s = (a + b + c) * 0.5;
    return sqrt(s * (s - a) * (s - b) * (s - c));
}

Aside from how numerically stable it is, the parentheses grouping the numbers are regular and elegant.

Archimedes, just how far ahead were you looking?

In the Distant Future, Medieval China Produces This Formulation

Qin Jiushao documented the following approach. For those familiar with mathematical history, Qin Jiushao was quite an accomplished mathematician.

double heron_asian(double a, double b, double c) {
    double a2 = a * a;
    double b2 = b * b;
    double c2 = c * c;
    double t = (a2 + c2 - b2) * 0.5;
    return 0.5 * sqrt(a2 * c2 - t * t);
}

It is slow. Substantially so.

Does this mean ink was wasted without value? Certainly not.

Each squared term unfolds step-by-step. One can follow along and appreciate visually why Heron's formula takes this form.

However, we do not use an abacus or counting rods. If computed via computer, it is slower than the rudimentary Archimedean formulation, leaving little reason to consider it as-is.

Perhaps by applying factoring and borrowing its core intuition, it can be improved.

Kahan's Method

double heron_kahan(double a, double b, double c) {
    double term1 = a + (b + c);
    double term2 = c - (a - b);
    double term3 = c + (a - b);
    double term4 = a + (b - c);

    double radicand = term1 * term2 * term3 * term4;
    return 0.25 * sqrt(fmax(0.0, radicand));
}

This well-known method feels right at home somewhere inside a low-level C program.

Each variable follows a termN pattern, decoupled so operations can execute independently in parallel.

Moreover, it proceeds using only additions and subtractions, ultimately multiplying all four together and dividing by 4.

A clean design philosophy that inherits classical Greek geometry.

My Style

My approach differs. It blends elements from both East and West. I tend to think a practitioner's cultural background influences algorithmic thinking. From that perspective, this might feel quite Korean.

double heron_mystyle(double a, double b, double c) {
    double inv_h = 1.0 / c;
    double h2 = c * c;

    a *= inv_h;
    b *= inv_h;
    c = 1.0;

    double a2 = a*a;
    double b2 = b*b;
    // 1.0 * 1.0 == 1.0
    // c2 is not necessary.
    // c2 == c

    double t = 0.25 * (a2 + c - b2);
    // c == 1.0
    double half_ac = a * 0.5;

    double l = half_ac - t;
    double r = half_ac + t;
    double rad = l * r;

    double area = sqrt(rad <= 0.0 ? 0.0 : rad);

    area *= h2;
    return area;
}

To prevent values from exploding, we scale down the triangle so the base side becomes 1.

Then, taking inspiration from Qin Jiushao's method through factoring, we factor (ac)^2 - t^2 = (ac - t)(ac + t) and provide boundary defense.

---

Starting from the original equation:

S = \frac{1}{2} \left[ a^2 c^2 - \left(\frac{a^2 + c^2 - b^2}{2}\right)^2 \right]^{1/2}

S = \left[ \frac{1}{4} a^2 c^2 - \left(\frac{a^2 + c^2 - b^2}{4}\right)^2 \right]^{1/2}

h = c

\tilde{a} = \frac{a}{c}, \quad \tilde{b} = \frac{b}{c}, \quad \tilde{c} = 1

S = h^2 \tilde{S}

\tilde{S} = \left[ \frac{1}{4}\tilde{a}^2(1)^2 - \left(\frac{\tilde{a}^2 + 1^2 - \tilde{b}^2}{4}\right)^2 \right]^{1/2}

\tilde{S} = \left[ \left(\frac{1}{2}\tilde{a}\right)^2 - \left(\frac{\tilde{a}^2 + 1 - \tilde{b}^2}{4}\right)^2 \right]^{1/2}

\text{half\_ac} = \frac{1}{2}\tilde{a}

t = \frac{1}{4}(\tilde{a}^2 + 1 - \tilde{b}^2)

\tilde{R} = (\text{half\_ac})^2 - t^2

\tilde{R} = (\text{half\_ac} - t)(\text{half\_ac} + t)

l = \text{half\_ac} - t

r = \text{half\_ac} + t

\tilde{R} = l \cdot r

\tilde{S} = [\max(0, \, l \cdot r)]^{1/2}

S = \tilde{S} \cdot h^2

In Kahan's method, there was no step to shrink the base side to 1 and scale it back up. While the formulation shifted somewhat, it preserves a measure of the readability found in Qin Jiushao's Chinese derivation.

Now, let us examine whether it genuinely defends against failure.

yjlee@elegant:~$ ./main5
1.0e100 1.0e100 1.0e100
==Classical Asian==
Result: -nan
121 ns elapsed
==Classical European==
Result: inf
30 ns elapsed
==My Style==
Result: 43301270189221932480949137768903387893892541093125184441720041813872084532521811805551529931486414663480764685208819956717969838856881289676610107485659692114536313096674369156159907634415957176221696.000000
40 ns elapsed
Baseline(Kahan's Algorithm)
Result: inf
30 ns elapsed

Area	Classical Asian	Classical European	My Style	Kahan
.	-nan	inf	\approx 4.330127 \times 10^{199}	inf

It naturally safeguards against extraordinarily large triangles.

However, if this protection costs an excessive runtime overhead, this implementation should not be deployed in production programs.

Elapsed Time	Classical Asian	Classical European	My Style	Baseline (Kahan)
20 ns	0 (0.0%)	17 (56.7%)	3 (10.0%)	14 (46.7%)
30~31 ns	0 (0.0%)	10 (33.3%)	24 (80.0%)	12 (40.0%)
40 ns	0 (0.0%)	2 (6.7%)	3 (10.0%)	2 (6.7%)
50~51 ns	6 (20.0%)	0 (0.0%)	0 (0.0%)	2 (6.7%)
60 ns	15 (50.0%)	1 (3.3%)	0 (0.0%)	0 (0.0%)
70~71 ns	7 (23.3%)	0 (0.0%)	0 (0.0%)	0 (0.0%)
80 ns	1 (3.3%)	0 (0.0%)	0 (0.0%)	0 (0.0%)
90 ns	1 (3.3%)	0 (0.0%)	0 (0.0%)	0 (0.0%)
Total	30	30	30	30
Average	62.8 ns	25.7 ns	30.1 ns	27.4 ns

It demands a sacrifice of roughly 10 nanoseconds.

Still, this implementation rarely exceeds 40 ns. It is predictable and fast enough.

If you must compute the area of extremely large triangles, it carries real consideration value.

Granted, This Might Not Immediately Go into Production, But...

I know this cannot be immediately plugged into production systems.

Yet when fed 10^{100}, an ostensibly flawless modern algorithm exposed a vulnerability in an unexpected place.

Perhaps we were merely mistaking a formula specialized for needle-like triangles as universally infallible?

Come to think of it, this can also be handled by branching with an if statement between "My Style," which cannot handle such special cases, and Kahan's method, which handles them well.

And once again, we are reminded:

Is something flawless everywhere simply because it is renowned and frequently used?

Even if outside the scope of daily production work, this experiment demonstrated one point clearly:

Even a formula devised by a brilliant pioneer was not an all-piercing spear capable of surmounting the final 0.001% edge case.

P.S

double heron_hybrid(double a, double b, double c) {
    if (c > 1e75) {
        return heron_mystyle(a, b, c);
    }
    return heron_kahan(c, b, a);
}

Wouldn't this just suffice? Let us remain adaptable, like water. As the ancient saying goes: "The highest goodness is like water" (上善若水).