For quick-and-dirty approximations, I'm fond of the musical logarithms (writeup by Sanjoy Mahajan of work due to I. J. Good, who credits his father). Essentially this boils down to a mathematical fact:
- $2^{10} \approx 10^3$, and taking 120th roots, $2^{1/12} \approx 10^{1/40}$.
and a "musical" fact:
- many rational numbers with small numerator and denominator can be approximated as powers of $2^{1/12}$.
I call this a "musical" fact because $2^{1/12}$ is the frequency ratio corresponding to an (equal-tempered) semitone.
For example: $3/2$ is the frequency ratio corresponding to the musical interval of a perfect fifth, which is seven semitones; thus $3/2 \approx 2^{7/12} \approx 10^{7/40}$ and so $\log_{10} 3/2 \approx 7/40$.
For the Newton-Raphson case, if the number were rounded to $8$ decimal places between iterations, then if you get the same number twice in a row, every successive number has to be the same after that. Since Newton-Raphson is quadratically convergent you normally get about twice the significant figures on every iteration after a pretty good approximation has been achieved. Of course you can make up examples where no convergence is reached or even force the algorithm to cycle between a limit set of values.
If you are tracking the values by hand you can see when you have over half the significant figures you want and perform the last iteration with extra precision to create a better probability of getting the last digit right. But you are right in thinking that it's a big problem to try to get the last digit right every time in floating point computations. In general you probably need about twice the significant figures of your final output to get the last digit right every time.
Let's look at an example with $\cosh x$ where we will run out all the $8$-digit numbers between $0$ and $1$ and then see how many might require $15$ digits to get the $8^{\text{th}}$ digit right. Out program searches for outputs with digits $9:15$ having values between $4999999$ and $5000001$. Since this is a range of $2$ out of $1\times10^7$, we might expect about $20$ hits for a program testing $1\times10^8$ inputs, and we are not far off in that estimate.
program round
use ISO_FORTRAN_ENV, only:wp=>REAL128,qp=>REAL128,wi=>INT64
implicit none
real(wp) x,y,z
real(qp) qx, qy
integer(wi) i
do i = 0,10_wi**8
x = 1.0e-8_wp*i
y = cosh(x)
z = y*1.0e8_wp+0.5_wp
if(abs(z-nint(z))<1.0e-7) then
qx = 1.0e-8_qp*i
qy = cosh(qx)
write(*,'(f10.8,1x,f22.20)') x,qy
end if
end do
end program round
Output:
0.00010000 1.00000000500000000417
0.00030000 1.00000004500000033750
0.05844226 1.00170823500000040322
0.07380594 1.00272489500000039410
0.44746231 1.10179282500000099007
0.45315675 1.10444463500000093431
0.47303029 1.11398059500000076845
0.47962980 1.11724437499999901204
0.49332468 1.12417258499999983893
0.49725888 1.12620181499999997563
0.51736259 1.13684395499999912340
0.59708506 1.18361444500000091792
0.60322956 1.18752751500000094494
0.64184055 1.21314873500000023866
0.72946242 1.27806675499999917355
0.75252442 1.29676328499999998383
0.90813720 1.44148689499999975793
0.93055850 1.46512925500000051185
The output for $0.00010000$ was expected from the Taylor series for $\cosh x$, but check out how close we got with $0.75252442$: we would have had to calculate the $17^{\text{th}}$ digit to round the $8^{\text{th}}$ digit correctly.
Best Answer
To add emphasis, the simple continued fraction of a rational number is finite. It can be found using computer operations that exclusively use integers. For your purpose, the most important part is that, if your rational number is $r,$ and we have consecutive "convergents" $p/q$ and $s/t,$ then $$ \left| r - \frac{p}{q} \right| \leq \frac{1}{qt} < \frac{1}{q^2}. $$ This means that, for $q > 10^9,$ we get $\left| r - \frac{p}{q} \right| < 10^{-18}$
https://en.wikipedia.org/wiki/Continued_fraction#Finite_continued_fractions
Alright, gp-pari has continued fractions. I wrote my own version in C++ with GMP, all I really wanted was the extended GCD, that is, given integers $a,b$ with $\gcd(a,b) = 1,$ find integers $x,y$ so $ax+by = 1.$
Where was I: to do many such fractions, all you need to do is implement the Euclidean Algorithm, with the extra step that, every time you calculate an intermediate partial quotient command "m / n" which is $$ \left\lfloor \frac{m}{n} \right\rfloor $$ you update the "convergent."