In a computer application, I need to solve trillions of times an equation which can be reduced to $$f(x)=\sin(x)-a x=0$$ Newton methods (quadratic and higher orders) are used for the solution. Parameter $a$ is random and the solution $x$ looked for is the one between $0$ and $\pi$ if it exists. The retained solution is immediate ($x=0$) if $a \geq 1$ or $a \leq 0$.
For the other cases ($0 \lt a \lt 1$), since I need to save as many iterations as I can, I focused on how to establish a good and unexpensive approximation of the solution. After some empirical experiments, what I found is that writing $$\sin(x) \simeq \frac{4}{\pi^2} x(\pi-x)$$ is a quite good approximation if $0 \leq a \leq 0.7$ leading to $$x \simeq \pi -\frac{\pi ^2 a}{4}$$
For the remaining interval, using Pade approximation $$\sin(x) \simeq \frac{x-\frac{7 x^3}{60}}{1+\frac{x^2}{20}}$$ which leads to $$x \simeq \frac{2 \sqrt{15} \sqrt{1-a}}{\sqrt{3 a+7}}$$ seems interesting.
I wonder if this could be improved, the goal being a simple explicit expression for the estimate of the solution.
Any idea and/or suggestion would really be welcome.
Added later
In comments and answers, Bhoot suggested me to look at the approximation $$\sin(x) \simeq \frac{16 (\pi -x) x}{5 \pi ^2-4 (\pi -x) x}$$ proposed by Mahabhaskariya of Bhaskara I, a seventh-century Indian mathematician. This spendid approximation leads to $$x \simeq \frac{2 \sqrt{-\pi ^2 a^2+2 \pi a+4}+\pi a-4}{2 a}$$ which is effectively very good except very close to $a=1$. This already makes a significant improvement.
Temporary improvement
Extremely impressed by the quality of the approximation proposed by Mahabhaskariya of Bhaskara I (more than 1400 years ago), I tried to undertand why it was very good except in the vicinity of $a=1$. I suspected that the derivative could be in error at the end points. Effectively, this formula gives a slope equal to $\frac{16}{5\pi} \simeq 1.01859$ instead of $1$. On the other side, the area under the curve is given by $$A=\pi \left(-4+\pi +\tan ^{-1}\left(\frac{3116}{237}\right)\right) \simeq 1.99955$$ So, modestly, I built an approximation which is $$\sin(x) \simeq \frac{\pi x(\pi-x)}{\pi^2+(\pi-4)x(\pi-x)}$$ which allows to match exactly the function and derivative values at $x=0,\frac{\pi}{2},\pi$; with respect to accuracy, it is not as good as the original showing a maximum error of $0.0052$ instead of $0.0016$. The area under the curve is then given by $$A=\frac{\pi \left(-4 \pi +\pi ^2+4 \sqrt{(4-\pi ) \pi } \tan ^{-1}\left(\sqrt{\frac{4}{\pi }-1}\right)\right)}{(\pi -4)^2} \simeq 1.99161$$ quite significantly worse than the original.
The estimate of the solution is given by $$x \simeq \frac{\pi \left((\pi -4) a-\sqrt{(\pi -4) a (\pi a-2)+1}+1\right)}{2 (\pi -4) a}$$ which makes Newton scheme converging in less than two iterations for the whole range.
Added after Christian Blatter's answer
I used what has been kindly proposed by Christian Blatter in his answer and set
$$\tilde f^2(a):={p(a)\over q(a)},\qquad p(a):=c_0+c_1 a+ c_2 a^2,\quad q(a):=1+d_1a +d_2 a^2\ $$ Using nonlinear regression, I adjusted the five involved parameters in order to minimize $$SSQ=\sum_{i=1}^n \Big(\sin(\tilde f(a_i))-a_i \tilde f(a_i)\Big)^2$$ The values of the $a_i$ were generated using $1000$ equally spaced values of the $x_i$ between $0$ and $\pi$. I have not been able to compute formally $$\int_0^1 \Big(\sin(\tilde f(a_i))-a_i \tilde f(a_i)\Big)^2 da$$
Starting with the coefficients given in Christian Blatter's answer, the initial $SSQ=5.968\times 10^{-5}$ which is already very good. I arrived to $SSQ= 3.800\times 10^{-6}$. The corresponding parameters are $$c_0=9.86774920$$ $$c_1=4.91765690$$ $$c_2=-14.77935381$$ $$d_1=2.48744104$$ $$d_2=0.63396306$$ For these values, the largest error is $0.000295$ and the average error is $0.000052$ which is incredibly good. As a result, a single Newton iteration is basically required for the desired accuracy. In the following plot the function $\tilde f$ is denoted $g$:
I would like to thank all people who contribute to this work. You have been extremely helpful.
Added later
Continuing working the problem, I set $$\tilde f^2(a):={p(a)\over q(a)},\qquad p(a):=\sum_{i=0}^n c_i a^i,\quad q(a):=1+\sum_{i=1}^n d_i a^i\ $$ and played with $n$. The first result is that moving to cubic polynomials changes a lot the result : $SSQ=4.9609379\times 10^{-10}$, maximum error $= 0.000004$, average error $< 0.000001$. This means that one single Newton iteration is required for a high accuracy. Mowing to fourth oder gives the solution without any Newton iteration.