"Are there any good indefinite (or much higher) precision packages"
A simple search on "arbitrary precision arithmetic in R" would have immediately gotten you to the Rmpfr package. See also: https://cran.r-project.org/web/packages/Ryacas/vignettes/arbitrary-precision.html -- Bert On Fri, Aug 18, 2023 at 4:58 PM <avi.e.gr...@gmail.com> wrote: > This discussion is sooo familiar. > > If you want indefinite precision arithmetic, feel free to use a language > and data type that supports it. > > Otherwise, only do calculations that fit in a safe zone. > > This is not just about this scenario. Floating point can work well when > adding (or subtracting) two numbers of about the same size. But if one > number is .123456789... and another is the same except raised to the -45th > power of ten, then adding them effectively throws away the second number. > > This is a well-known problem for any finite binary representation of > numbers. In the example given, yes, the smaller the number is, the worse > the behavior in this case tends to be. > > There are many solutions and some are fairly expensive in terms of > computation time and sometimes memory usage. > > Are there any good indefinite (or much higher) precision packages out > there that would not only support the data type needed but also properly be > used and passed along to the functions used to do complex calculations? No, > I am not asking for indefinite precision complex numbers, but generally > that would be a tuple of such numbers. > > > -----Original Message----- > From: R-help <r-help-boun...@r-project.org> On Behalf Of Bert Gunter > Sent: Friday, August 18, 2023 7:06 PM > To: Leonard Mada <leo.m...@syonic.eu> > Cc: R-help Mailing List <r-help@r-project.org>; Martin Maechler < > maech...@stat.math.ethz.ch> > Subject: Re: [R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2 > > "The ugly thing is that the error only gets worse as x decreases. The > value neither drops to 0, nor does it blow up to infinity; but it gets > worse in a continuous manner." > > If I understand you correctly, this is wrong: > > > x <- 2^(-20) ## considerably less then 1e-4 !! > > y <- 1 - x^2/2; > > 1/(1 - y) - 2/x^2 > [1] 0 > > It's all about the accuracy of the binary approximation of floating point > numbers (and their arithmetic) > > Cheers, > Bert > > > On Fri, Aug 18, 2023 at 3:25 PM Leonard Mada via R-help < > r-help@r-project.org> wrote: > > > I have added some clarifications below. > > > > On 8/18/2023 10:20 PM, Leonard Mada wrote: > > > [...] > > > After more careful thinking, I believe that it is a limitation due to > > > floating points: > > > [...] > > > > > > The problem really stems from the representation of 1 - x^2/2 as shown > > > below: > > > x = 1E-4 > > > print(1 - x^2/2, digits=20) > > > print(0.999999995, digits=20) # fails > > > # 0.99999999500000003039 > > > > The floating point representation of 1 - x^2/2 is the real culprit: > > # 0.99999999500000003039 > > > > The 3039 at the end is really an error due to the floating point > > representation. However, this error blows up when inverting the value: > > x = 1E-4; > > y = 1 - x^2/2; > > 1/(1 - y) - 2/x^2 > > # 1.215494 > > # should be 1/(x^2/2) - 2/x^2 = 0 > > > > > > The ugly thing is that the error only gets worse as x decreases. The > > value neither drops to 0, nor does it blow up to infinity; but it gets > > worse in a continuous manner. At least the reason has become now clear. > > > > > > > > > > Maybe some functions of type cos1p and cos1n would be handy for such > > > computations (to replace the manual series expansion): > > > cos1p(x) = 1 + cos(x) > > > cos1n(x) = 1 - cos(x) > > > Though, I do not have yet the big picture. > > > > > > > Sincerely, > > > > > > Leonard > > > > > > > > > > > On 8/17/2023 1:57 PM, Martin Maechler wrote: > > >>>>>>> Leonard Mada > > >>>>>>> on Wed, 16 Aug 2023 20:50:52 +0300 writes: > > >> > Dear Iris, > > >> > Dear Martin, > > >> > > >> > Thank you very much for your replies. I add a few comments. > > >> > > >> > 1.) Correct formula > > >> > The formula in the Subject Title was correct. A small glitch > > >> swept into > > >> > the last formula: > > >> > - 1/(cos(x) - 1) - 2/x^2 > > >> > or > > >> > 1/(1 - cos(x)) - 2/x^2 # as in the subject title; > > >> > > >> > 2.) log1p > > >> > Actually, the log-part behaves much better. And when it fails, > > >> it fails > > >> > completely (which is easy to spot!). > > >> > > >> > x = 1E-6 > > >> > log(x) -log(1 - cos(x))/2 > > >> > # 0.3465291 > > >> > > >> > x = 1E-8 > > >> > log(x) -log(1 - cos(x))/2 > > >> > # Inf > > >> > log(x) - log1p(- cos(x))/2 > > >> > # Inf => fails as well! > > >> > # although using only log1p(cos(x)) seems to do the trick; > > >> > log1p(cos(x)); log(2)/2; > > >> > > >> > 3.) 1/(1 - cos(x)) - 2/x^2 > > >> > It is possible to convert the formula to one which is > > >> numerically more > > >> > stable. It is also possible to compute it manually, but it > > >> involves much > > >> > more work and is also error prone: > > >> > > >> > (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x))) > > >> > And applying L'Hospital: > > >> > (2*x - 2*sin(x)) / (2*x * (1 - cos(x)) + x^2*sin(x)) > > >> > # and a 2nd & 3rd & 4th time > > >> > 1/6 > > >> > > >> > The big problem was that I did not expect it to fail for x = > > >> 1E-4. I > > >> > thought it is more robust and works maybe until 1E-5. > > >> > x = 1E-5 > > >> > 2/x^2 - 2E+10 > > >> > # -3.814697e-06 > > >> > > >> > This is the reason why I believe that there is room for > > >> improvement. > > >> > > >> > Sincerely, > > >> > Leonard > > >> > > >> Thank you, Leonard. > > >> Yes, I agree that it is amazing how much your formula suffers from > > >> (a generalization of) "cancellation" --- leading you to think > > >> there was a problem with cos() or log() or .. in R. > > >> But really R uses the system builtin libmath library, and the > > >> problem is really the inherent instability of your formula. > > >> > > >> Indeed your first approximation was not really much more stable: > > >> > > >> ## 3.) 1/(1 - cos(x)) - 2/x^2 > > >> ## It is possible to convert the formula to one which is numerically > > >> more > > >> ## stable. It is also possible to compute it manually, but it > > >> involves much > > >> ## more work and is also error prone: > > >> ## (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x))) > > >> ## MM: but actually, that approximation does not seem better (close > > >> to the breakdown region): > > >> f1 <- \(x) 1/(1 - cos(x)) - 2/x^2 > > >> f2 <- \(x) (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x))) > > >> curve(f1, 1e-8, 1e-1, log="xy" n=2^10) > > >> curve(f2, add = TRUE, col=2, n=2^10) > > >> ## Zoom in: > > >> curve(f1, 1e-4, 1e-1, log="xy",n=2^9) > > >> curve(f2, add = TRUE, col=2, n=2^9) > > >> ## Zoom in much more in y-direction: > > >> yl <- 1/6 + c(-5, 20)/100000 > > >> curve(f1, 1e-4, 1e-1, log="x", ylim=yl, n=2^9) > > >> abline(h = 1/6, lty=3, col="gray") > > >> curve(f2, add = TRUE, n=2^9, col=adjustcolor(2, 1/2)) > > >> > > >> Now, you can use the Rmpfr package (interface to the GNU MPFR > > >> multiple-precision C library) to find out more : > > >> > > >> if(!requireNamespace("Rmpfr")) install.packages("Rmpfr") > > >> M <- function(x, precBits=128) Rmpfr::mpfr(x, precBits) > > >> > > >> (xM <- M(1e-8))# yes, only ~ 16 dig accurate > > >> ## 1.000000000000000020922560830128472675327e-8 > > >> M(10, 128)^-8 # would of course be more accurate, > > >> ## but we want the calculation for the double precision number 1e-8 > > >> > > >> ## Now you can draw "the truth" into the above plots: > > >> curve(f1, 1e-4, 1e-1, log="xy",n=2^9) > > >> curve(f2, add = TRUE, col=2, n=2^9) > > >> ## correct: > > >> curve(f1(M(x, 256)), add = TRUE, col=4, lwd=2, n=2^9) > > >> abline(h = 1/6, lty=3, col="gray") > > >> > > >> But, indeed we take note how much it is the formula instability: > > >> Also MPFR needs a lot of extra bits precision before it gets to > > >> the correct numbers: > > >> > > >> xM <- c(M(1e-8, 80), M(1e-8, 96), M(1e-8, 112), > > >> M(1e-8, 128), M(1e-8, 180), M(1e-8, 256)) > > >> ## to and round back to 70 bits for display: > > >> R <- \(x) Rmpfr::roundMpfr(x, 70) > > >> R(f1(xM)) > > >> R(f2(xM)) > > >> ## [1] 0 0 > > >> 0.15407439555097885670915 > > >> ## [4] 0.16666746653133802175779 0.16666666666666666749979 > > >> 0.16666666666666666750001 > > >> > > >> ## 1. f1() is even worse than f2() {here at x=1e-8} > > >> ## 2. Indeed, even 96 bits precision is *not* sufficient at all, ... > > >> ## which is amazing to me as well !! > > >> > > >> Best regards, > > >> Martin > > > > ______________________________________________ > > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
