A 1400 years old approximation to the sine function by Mahabhaskariya of Bhaskara I

Question

The approximation $$\sin(x) \simeq \frac{16 (\pi -x) x}{5 \pi ^2-4 (\pi -x) x}\qquad (0\leq x\leq\pi)$$ was proposed by Mahabhaskariya of Bhaskara I, a seventh-century Indian mathematician.

I wondered how much this could be improved using our computers and so I tried (very immodestly) to see if we could do better using $$\sin(x) \simeq \frac{a (\pi -x) x}{5 \pi ^2-b (\pi -x) x}$$ I so computed $$\Phi(a,b)=\int_0^{\pi} \left(\sin (x)-\frac{a (\pi -x) x}{5 \pi ^2-b (\pi -x)x}\right)^2 dx$$ the analytical expression of which not being added to the post. Settings the derivatives equal to $0$ and solving for $a$ and $b$, I arrived to $a=15.9815,b=4.03344$ so close to the original approximation !

What is interesting is to compare the values of $\Phi$ : $2.98 \times 10^{-6}$ only decreased to $2.17 \times 10^{-6}$. Then, no improvement and loss of attractive coefficients.

Now, since this is a matter of etiquette on this site, I ask a simple question:

with all the tools and machines we have in our hands, could any of our community propose something as simple (or almost) for basic trigonometric functions ?

In the discussions, I mentioned one I made (it is probable that I reinvented the wheel) in the same spirit $$\cos(x) \simeq\frac{\pi ^2-4x^2}{\pi ^2+x^2}\qquad (-\frac \pi 2 \leq x\leq\frac \pi 2)$$ which is amazing too !

Answer 1

While you're at it, try also $\cos\bigg(\dfrac\pi2x\bigg)\simeq\Big(1-x^2\Big)^\tfrac65$ and $\cos\bigg(\dfrac\pi2x\bigg)\simeq\Big(1-x^2\Big)^\tfrac76$.

But since the numerical evaluation of fractional powers is significantly more time-

consuming in terms of CPU, we can substantially improve this by using the binomial

series for $\Big(1-x^2\Big)^\tfrac15$, and experimentally adjusting the coefficient, finally arriving at

$\color{seagreen}{\cos\bigg(\dfrac\pi2x\bigg)\simeq\Big(1-x^2\Big)\bigg(1-\dfrac{x^2}{4.5}\bigg)}$, which yields an absolute error of about $\pm1$ ‰

Answer 2

This might be more explicable if you observe that it is the same thing as

$$ \csc(x) \approx -\frac{1}{4} + \frac{5 \pi}{16} \left( \frac{1}{x} + \frac{1}{\pi - x} \right) $$

The two summands in the parentheses are obvious if you want to get the poles of $\csc$ correct. If you wanted a good approximation of $\csc$ near the poles, then the coefficient out front should be $1$. But since we're approximating $\sin$, it's okay to get that wrong because anything near zero is near zero.

The extreme point is at $\csc(\pi/2) = 1$; in the approximation, this would become

$$ -\frac{1}{4} + \frac{5 \pi}{16} \left( \frac{2}{\pi} + \frac{2}{\pi} \right) = -\frac{1}{4} + \frac{5}{4} = 1$$

and so we see the appearance of the remaining copy of $\pi$ is to cancel out the other two $\pi$'s. All that's left is to tune the factor $\frac{5}{16}$ to something appropriate, and adjust the $-\frac{1}{4}$ to compensate. I'm not sure where the choice of $\frac{5}{16}$ comes from, although it's quite plausible it ought to be near $\frac{1}{\pi}$; maybe it was chosen just to be a small fraction whose denominator was divisible by $4$, so as to cancel the $4$ in $\frac{4}{\pi}$.

---

As a bit of an aside, my comment about the poles suggests an infinite sum for $\csc(x)$ that I hadn't seen before:

$$ \csc(x) = \sum_{n=-\infty}^{+\infty} (-1)^n \frac{1}{x - \pi n} $$

Answer 3

One simple way to derive this is to come up with a parabola approximation. Just getting the roots correct we have

$$f(x)=x(\pi-x)$$

Then, we need to scale it (to get the heights correct). And we are gonna do that by dividing by another parabola $p(x)$

$$f(x)=\frac{x(\pi-x)}{p(x)}$$

Let's fix this at three points (thus defining a parabola). Easy rational points would be when $\sin$ is $1/2$ or $1$. So we fix it at $x=\pi/6,\pi/2,5\pi/6$.

We want $$f(\pi/6)=f(5\pi/6)=1/2=\frac{5\pi^2/36}{p(\pi/6)}=\frac{5\pi^2/36}{p(5\pi/6)}$$ And we conclude that $p(\pi/6)=p(5\pi/6)=5\pi^2/18$

We do the same at $x=\pi/2$ to conclude that $p(\pi/2)=\pi^2/4$.

The only parabola through those points is

$$p(x)=\frac{1}{16}(5\pi^2-4x(\pi-x))$$

And thus we have the original approximation.

In the spirit of answering the question: This method could be applied for most trig functions on some small symmetric bound.

Answer 4

I wondered how much this could be improved using our computers

Here is an overview of how much Bhāskara's formula can be improved, given as abolute and relative errors of the respective formulas, and as %-values relative to Bhāskara's approximation.

Approximation of $\cos x$ for $|x|\leqslant\frac\pi2$	Absolute Error	Relative Error
Bhāskara I's $(\pi^2-4x^2)/(\pi^2+x^2)$	0.00163 (100%)	0.019 (100%)
Best degree [2/2] rational function	0.00098 (60%)	0.0047 (25%)
Best degree 4 polynomial	0.0006 (37%)	0.0027 (15%)

The values and functions are worked out in the remainder.

Absolute Error

If we want to minimize the absolute error, it's simpler to start with approximating $\cos$ for $|x|\leqslant \pi/2$. We start by approximating

$$ f(x) = \cos(\sqrt x) \quad \text{ with }\quad x\in[0,\pi^2/4] \tag 1$$ and then use $$ \cos x = f(x^2) \tag2 $$ That way, we get the symmetry of $\cos$ and $f$ can be approximated by a rational function $g$ of degree [1/1]. Doing the numerics (Remez algorithm), one finds

$$ g(x) = \frac{p(x)}{q(x)} = \frac{9.8473383 - 3.9958133 x}{9.8377958 + x} \tag 3 $$

where the function is normalized so that coefficient $q_1=1$. The coefficients of the numerator and denominator polynomial only loosely resemble the values of $\cos \approx(\pi^2-4x^2)/(\pi^2+x^2)$ worked out in the other answer.

The absolte error is given by

$$ \max_{x\in[0,\pi^2/4]} |g(x)-f(x)| = \max_{|x|\leqslant \pi/2} |g(x^2)-\cos(x)| < 0.00097 \tag 4$$ which is an accuracy of around 3 decimals or almost exactly 10 bits.

Plot: Absolute error of $g(x^2)$ against $\cos(x)$ for $|x|\leqslant\pi/2$.

The optimality of $g(x)$ (or of $g(x^2)$ for that matter) follows from the alternant property of the absolute error: If $\Delta=\max|g-f|$ the absolute error oscillates between $+\Delta$ and $-\Delta$ which is nicely demonstrated by the plot above.

In contrast, the non-optimality of the classic

$$ \frac{\pi^2-4x^2}{\pi^2+x^2} \tag 5$$

is immediately clear from the plot below because the absolue error does not oscillate nicely with constant "amplitude".

Plot: Absolute error of $(π²−4x²)/(π²+x²)$ against $\cos(x)$ for $|x|\leqslant\pi/2$.

Conclusion

• The classic $(\pi^2-4x^2)/(\pi^2+x^2)$ has an absolute error of ca. $1.63\cdot 10^{-3}$, whereas the optimal $$ g(x^2) = \frac{9.8473383 - 3.9958133 x^2}{9.8377958 + x^2} \tag6 $$ has an absolute error of less than $10^{-3}$.

• $g(x^2)$ is the best approximation in terms of absolute error amongst all rational functions of degree [2/2] or less.

• Desmos plot of approximations and their absolute errors

Relative Error

For optimizing the relative error, the situation is similar: The classic formula is sub-optimal, and we can numerically approch the best possible rational function of degree [2/2].

As $\cos$ has zeros at $\pm\pi/2$, we have to divide them out so that the relative error is well-behaved. Using Remez, we find:

$$ \frac{\cos x}{(\frac\pi2-x)(\frac\pi2+x)} \approx r(x)=\frac{3.676844}{9.030178 + x^2} \tag7 $$ with a relative error of $0.0047$ over the target interval as shown by the following plot.

Hence, the best approximation in terms of relative error amongst all rational functions of degree [2/2] or less is:

$$ r(x)\left(\frac{\pi^2}4-x^2\right) = \frac{3.676844}{9.030178 + x^2}\left(\frac{\pi^2}4-x^2\right) \tag8 $$

The relative error of the classic approximaton is around $0.019$, which is 4 times the relative error of $(8)$.

Desmos plot of approximations and their relative errors

Improving the Results by using Polynomials of Degree 4

The best approximation in terms of absolute error and degree [2/2] is $g(x^2)$ given by $(6)$ with an absolute error around below 0.001.

Evaluating $g$ costs: 1 division, 2 multiplications and 2 additions.

The question is whether this result can be improved with the same or less compute. This is actually the case:

$$ 0.99940323 + x^2 (-0.49558085 + 0.036791683 x^2) \tag 9 $$

costs 3 multiplications and 2 additions. As divisions are usually not cheaper than multiplications, this is an improvement compute-wise, and the absolut error is below 0.0006 which is just 60% of the absolute error of $g$ and just 37% of the absolute error of Bhāskara's $(5)$.

The situation with relative error can also be improved:

$$ (0.4042199 - 0.03515697 x^2)\left(\frac{\pi^2}4 - x^2\right) \tag{10} $$ has a relative error of less than 0.0027 over the target interval, which is only 58% of the error of $(8)$, and only 15% of the relative error of Bhāskara's $(5)$.

Both $(8)$ and $(10)$ have the alternante property, which means that these polynomials are the best ones amongst all polynomials of degree 4 or less.

• Desmos plot of the errors of $(8)$ and $(10)$

Answer 5

I would like to add an addendum to @BeaumontTaz's wonderful derivation of the approximation, because I think I have a very simple answer to user @kaka's comment "What is the intuition behind scaling of $x(\pi−x)$ by dividing with $p(x)$, and not by some multiplicative function ?". The key prerequisite to this approximation is I think the knowledge of the five nicest/simplest values of the sine function on $[0,\pi]$, namely $\sin(0) = 0, \sin(\frac \pi 6) = \frac 12, \sin(\frac \pi 2) = 1, \sin (\frac{5\pi}6) = \frac 12, \sin (\pi) = 0$. It is no surprise that these values were known to Bhaskara I in the 7th century, as all the values are rational (none of that square-rooty business that apocryphally got someone killed), and can be easily derived through the 30-60-90 triangle.

Anyways, as BeaumontTaz did, the first order of business is to approximate using the zeroes, so $f(x)= x (\pi-x)$. Then obviously, we want $f(\frac \pi 2) = \sin (\frac \pi 2) = 1$, so the easiest thing we could do is scale by $\frac 1{f(\pi/2)}$ and get the approximation $$a_1(x) = \frac{x(\pi-x)}{f(\pi/2)}.$$ This is not too good of an approximation, so we hope to do better. Lucky for us, we still have two pieces of information left! We want our approximation to agree with $\sin x$ at $x = \frac \pi6, \frac{5\pi}6$ as well. So at $x=\frac \pi6$, we would need to scale by $\frac{\sin(\pi/6)}{f(\pi/6)} = \frac{1}{2f(\pi/6)}$, and similarly at $x= \frac{5\pi}6$ we would scale by $\frac{\sin(5\pi/6)}{f(5\pi/6)} = \frac{1}{2f(5\pi/6)}$. In summary, at $x = \frac \pi6, \frac \pi2, \frac{5\pi}6$ respectively, we would like to divide by $2f(\pi/6), f(\pi/2), 2f(5\pi/6)$ respectively, i.e. we have a "scaling function" $s(x)$ that we divide by that we have pinned down at three values, and an approximation $a(x) = \frac{f(x)}{s(x)}$. What's the easiest way to extend a function defined at three points? Well, there is exactly one parabola passing through those three points, and parabolas are very easy to deal with, so let us simply take $s(x)$ to be the parabola passing those three points.

So we see that the reason it is "more natural" to divide by the scaling instead of multiplying like @kaka suggested was because in this problem the first idea we had of scaling was in the form of a division ("priming our mind for division"), and all the scaling factors looked like $\frac{1}{\text{something}}$, and it is easier to deal with/look at interpolation at points $(\frac \pi 6, 2f(\frac \pi 6)), (\frac \pi 2, f(\frac \pi2)), (\frac{5\pi}6, 2f(\frac{5\pi}6))$ than interpolation at points $(\frac \pi 6, \frac{1}{2f(\frac \pi 6)}), (\frac \pi 2, \frac{1}{f(\frac \pi2)}), (\frac{5\pi}6, \frac1{2f(\frac{5\pi}6)})$. Said another way, when we want to "renormalize" $f(\frac \pi 6)$ to $\frac 12$, I think people naturally think "divide by $f(\frac \pi 6)$ to normalize to 1, and then divide by 2/multiply by $\frac 12$ to get to $\frac 12$", i.e. the idea of "normalization to 1" is "more naturally" thought of as a division rather than a multiplication by reciprocal.

I do think aesthetics/simplicity play some role in formulas like these (they were written out in verses in Sanskrit), and the formula for the parabola passing through $(\frac \pi 6, 2f(\frac \pi 6)), (\frac \pi 2, f(\frac \pi2)), (\frac{5\pi}6, 2f(\frac{5\pi}6))$ is much nicer (especially if we write things in terms of degrees instead of radians) than the formula for the parabola passing through $(\frac \pi 6, \frac{1}{2f(\frac \pi 6)}), (\frac \pi 2, \frac{1}{f(\frac \pi2)}), (\frac{5\pi}6, \frac1{2f(\frac{5\pi}6)})$ (taking reciprocals in the latter case messes up things quite a lot). Indeed explicitly writing out the interpolation w.r.t. the first and second set of points in degrees, we have $s_1(x) = \frac{1}{4}x^{2}-45x+10125$ and $s_2(x)=-\frac{1}{291600000}x^{2}+\frac{1}{1620000}x+\frac{31}{324000}$; you can visualize them here: https://www.desmos.com/calculator/aiqcjbbcn2.

I hope this answer gave a little bit more details to the "naturality" of dividing by scaling factors (and the "naturality" of the entire construction/approximation --- basically I argue that the formula in Claude's question is the simplest and most reasonable approximation given exactly the five "trivial" values of $\sin x$), and didn't just repeat what people already knew.