[R] About the multiprecision computing package in R

Richard M. Heiberger rmh @end|ng |rom temp|e@edu
Sun Mar 15 00:55:11 CET 2020


Rmpfr product and inner product are not doubling the nominal precision.
Small precision inner products are stored as 53-bit numeric.

library(Rmpfr)

 ## double with precBits=53
(5/17) * (5/17) ## 0.08650519 ## base 10
formatHex( (5/17) * (5/17) ) ##  +0x1.6253443526171p-4

## various precisions
fs5 <- 5*mpfr(1, precBits=5)/17
fs6 <- 5*mpfr(1, precBits=6)/17
fs7 <- 5*mpfr(1, precBits=7)/17
fs8 <- 5*mpfr(1, precBits=8)/17
fs9 <- 5*mpfr(1, precBits=9)/17

## product
formatHex(fs5*fs5) ## +0x1.7p-4

## All these are smaller
formatHex(fs6*fs6)
formatHex(fs7*fs7)
formatHex(fs8*fs8)
formatHex(fs9*fs9)

## and all round DOWN to +0x1.6p-4
formatHex(mpfr(  fs6*fs6 , precBits=5))
formatHex(mpfr(  fs7*fs7 , precBits=5))
formatHex(mpfr(  fs8*fs8 , precBits=5))
formatHex(mpfr(  fs9*fs9 , precBits=5))


## inner product, for low precision, is stored as 53-bit numeric
rounded to 6 bits
formatHex(fs5 %*% fs5) ## +0x1.6900000000000p-4
class(fs5 %*% fs5)
getPrec(fs5 %*% fs5)
formatHex((fs5 %*% fs5)*(fs5 %*% fs5)) ## 12 bits ## +0x1.fd11000000000p-8

## All these are smaller
formatHex(fs6 %*% fs6) ## same
formatHex(fs7 %*% fs7)
formatHex(fs8 %*% fs8)
formatHex(fs9 %*% fs9)

## and all round DOWN to +0x1.6p-4
formatHex(mpfr(  fs6 %*% fs6 , precBits=5)) ## same
formatHex(mpfr(  fs7 %*% fs7 , precBits=5))
formatHex(mpfr(  fs8 %*% fs8 , precBits=5))
formatHex(mpfr(  fs9 %*% fs9 , precBits=5))



## inner product, for high precision, stays at the specified precision
fs65 <- 5*mpfr(1, precBits=65)/17
formatHex(fs65 * fs65)   ## +0x1.62534435261707f8p-4

formatHex(fs65 %*% fs65) ## +0x1.62534435261707f8p-4
class(fs65 %*% fs65)
getPrec(fs65 %*% fs65)

fs130 <- 5*mpfr(1, precBits=130)/17
formatHex(fs130 * fs130)   ## +0x1.62534435261707f8e9dacbbcad9e8f800p-4
formatHex(mpfr( fs130 * fs130 , precBits=65))   ## +0x1.62534435261707f9p-4
## not the same as 65-bit inner product

On Sat, Mar 14, 2020 at 4:51 PM J C Nash <profjcnash using gmail.com> wrote:
>
> The issue is to avoid the storage and operational penalty. 100 x 100 matrix in 100 decimals vs 100 x 100 matrix in 50
> decimals for many operations like copy, scale, etc. But accumulation of inner products, you want to avoid digit loss,
> e.g., A and B are long vectors -- say 100000 long, with a few "large" elements but most smaller. However 99000 "small"
> numbers could add up and ...
>
> As with many such issues, it doesn't affect the vast majority of occurrences. Indeed few of my own. IEEE arithmetic
> usually does the job (I was one of the members of the original committee). This particular
> concern arose when trying to tease out a signal from a great deal of noise from a spacecraft sensor.
>
> JN
>
>
> On 2020-03-14 4:28 p.m., Jeff Newmiller wrote:
> > Not sure I understand the concern. IEEE 754 double precision floating point was invented to allow for avoiding loss of precision when manipulating single precision floating point numbers... but then C just ignores single precision and you are expected to know that the precision of your answers may only be relied on at single precision for certain operations.
> >
> > Wouldn't you just set Rmpfr precision to double your actual desired precision and move on? Though I suppose you might consider  more than doubling the desired precision to deal with exponentiation [1].
> >
> > [1] https://en.m.wikipedia.org/wiki/Extended_precision#Working_range
> >
> > On March 14, 2020 1:10:19 PM PDT, J C Nash <profjcnash using gmail.com> wrote:
> >> Rmpfr does "support" matrix algebra, but I have been trying for some
> >> time to determine if it computes "double" precision (i.e., double the
> >> set level of precision) inner products. I suspect that it does NOT,
> >> which is unfortunate. However, I would be happy to be wrong about
> >> this.
> >>
> >> JN
> >>
> >> On 2020-03-14 3:41 p.m., Bert Gunter wrote:
> >>> Read its documentation yourself and unless you have good reason not
> >> to,
> >>> always cc the list (which I have done here).
> >>>
> >>>
> >>> Bert Gunter
> >>>
> >>> "The trouble with having an open mind is that people keep coming
> >> along and
> >>> sticking things into it."
> >>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
> >>>
> >>>
> >>> On Sat, Mar 14, 2020 at 12:28 PM 林伟璐 <13917987541 using 163.com> wrote:
> >>>
> >>>> Thanks. Does it support matrix algebra?
> >>>>
> >>>>
> >>>>
> >>>> 林伟璐
> >>>> 邮箱:13917987541 using 163.com
> >>>>
> >>>>
> >> <https://maas.mail.163.com/dashi-web-extend/html/proSignature.html?iconUrl=http%3A%2F%2Fmail-online.nosdn.127.net%2F17c23bc1722125aa6261c33736f525c5.jpg&name=%E6%9E%97%E4%BC%9F%E7%92%90&uid=13917987541%40163.com&ftlId=1&items=%5B%22%E9%82%AE%E7%AE%B1%EF%BC%9A13917987541%40163.com%22%5D>
> >>>>
> >>>> 签名由 网易邮箱大师 <https://mail.163.com/dashi/dlpro.html?from=mail88> 定制
> >>>>
> >>>> On 03/15/2020 03:18, Bert Gunter <bgunter.4567 using gmail.com> wrote:
> >>>> Use google instead, as I recommended. If that's impossible in China,
> >> you
> >>>> should state that and what you tried (Baidu) in your query.
> >>>>
> >>>> You'll get the Rmpfr package.
> >>>>
> >>>> Bert Gunter
> >>>>
> >>>> "The trouble with having an open mind is that people keep coming
> >> along and
> >>>> sticking things into it."
> >>>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
> >>>>
> >>>>
> >>>> On Sat, Mar 14, 2020 at 11:48 AM 林伟璐 <13917987541 using 163.com> wrote:
> >>>>
> >>>>> Thanks, I tried Baidu but get nothing
> >>>>>
> >>>>>
> >>>>>
> >>>>> 林伟璐
> >>>>> 邮箱:13917987541 using 163.com
> >>>>>
> >>>>>
> >> <https://maas.mail.163.com/dashi-web-extend/html/proSignature.html?iconUrl=http%3A%2F%2Fmail-online.nosdn.127.net%2F17c23bc1722125aa6261c33736f525c5.jpg&name=%E6%9E%97%E4%BC%9F%E7%92%90&uid=13917987541%40163.com&ftlId=1&items=%5B%22%E9%82%AE%E7%AE%B1%EF%BC%9A13917987541%40163.com%22%5D>
> >>>>>
> >>>>> 签名由 网易邮箱大师 <https://mail.163.com/dashi/dlpro.html?from=mail88> 定制
> >>>>>
> >>>>> On 03/15/2020 02:44, Bert Gunter <bgunter.4567 using gmail.com> wrote:
> >>>>> Here's a novel idea:
> >>>>> Do a google search on "multiprecision computing package R" for an
> >> answer.
> >>>>>
> >>>>> Bert Gunter
> >>>>>
> >>>>> "The trouble with having an open mind is that people keep coming
> >> along
> >>>>> and sticking things into it."
> >>>>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
> >>>>>
> >>>>>
> >>>>> On Sat, Mar 14, 2020 at 10:36 AM 林伟璐 <13917987541 using 163.com> wrote:
> >>>>>
> >>>>>> Dear all
> >>>>>>
> >>>>>>
> >>>>>> I need a multiprecision computing package in R, if anyone in the
> >> list
> >>>>>> knows, please let me known...
> >>>>>>
> >>>>>>
> >>>>>> Many thanks
> >>>>>>
> >>>>>>
> >>>>>> Weilu Lin
> >>>>>>         [[alternative HTML version deleted]]
> >>>>>>
> >>>>>> ______________________________________________
> >>>>>> R-help using 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 using 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 using 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 using 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.



More information about the R-help mailing list