[R] Rmpfr question

[Martin Maechler]

>>>>> Martin Maechler <maechler at lynne.ethz.ch>
>>>>>     on Wed, 25 Sep 2013 11:08:07 +0200 writes:

>>>>> Michel  <michelgomez at free.fr>
>>>>>     on Tue, 24 Sep 2013 12:22:10 +0200 writes:

    >> Hello, Thanks for your answer The file does not contains
    >> numbers in high precision but all the calculation applied
    >> to these data will be

    >> In attachment a text file containing some lines And her
    >> few values:
    >> 11111.0054014388224326026488598,-68239.4114845811819662912967033,10053.1824584878233990181684021
    >> 3.05363891777017837760380136736,-1.40443175474415104684797195311,1.36766960877453817890022497172

    > (two lines, typically split by mail writers/readers)

    > Aha!  But these actually *ARE* of  'high precision', i.e.
    > you cannot store them as usual double precision numbers in full accuracy.

    > So, I've prepared the following -- 100% reproducible -- script
    > to show you how to get such numbers into an Rmpfr matrix :


    > ## MM: Reproducible example
    > set.seed(17); x <- mpfr(matrix(rnorm(28), 7, 4), precBits=128)
    > x <- x^2
    > mfile <- tempfile()
    > write.table(array(format(x), dim=dim(x)), file = mfile,
    > row.names=FALSE, col.names=FALSE, sep=",")
    > ## to check:
    > writeLines(readLines(mfile) [1:2])


    > ## Now, let's assume mfile contains the "high precision" matrix we want to get as
    > ## mpfrMatrix :
    > ## 1) Assume you know the number of columns, then this is fastest :
    > m.ncol <- 4
    > chmat <- matrix(scan(mfile, "", sep=","), ncol = m.ncol )

and -- oops! -  the above was missing the ominous 'byrow = TRUE' 
i.e. the above needs to be

    chmat <- matrix(scan(mfile, "", sep=","), ncol = m.ncol, byrow = TRUE)

... which makes the detour via read.table() even more attractive ...

    > ## 2) Nothing is known, ... ok, why not make the detour via read.table():
    > chmat <- as.matrix(read.table(mfile, colClasses="character", sep=","))

    > ## in both cases:

    > require(Rmpfr)# the package

    > M <- mpfr(chmat)#-> determines precision *from* the input, here  133 .. 143 bits

    > ## or set the precision yourself, high enough:
    > M <- mpfr(chmat, precBits = 144)

    > ## and e.g. this works:

    > crossprod(M) ## == M'M

    > -------

    > Hoping this helps,
    > Martin