1 -*- org -*--> Emacs [Tab] key + [Org] menu; C-c C-o follows links 2* Very Short Term 3** TODO 61) Q: Why is quantile() so slow [hence summary() very slow!]? A: because sort() |--> rank() is so slow (in C code!) 4 l2x <- seqMpfr(mpfr(4, 1024), 513, by=1/16) # 8000 numbers of 1024 bits 5 system.time(ql2x <- quantile(l2x, names=FALSE)) # user: 10.735 (nb-mm5, Nov.2020) 6*** quantile() -> sort(*, partial=.) -> xtfrm.default() -> rank(l2x) is so slow 7*** Partial SOLUTION (not yet implemented): use is.unsorted(.) which is *fast* FALSE for sorted vetors (as 'l2x' above) 8** TODO 55b) possibly more documentation on .mpfr* functions, e.g. .getSign(), at least *internally* (roxygen) 9** TODO 31) Valgrind problems + leaks: Brian's e-mail 2014-06-19; ~/R/Pkgs/Rmpfr.Rcheck_valgrind/ 10** TODO 11) format() method for "mpfr", "mpfrArray" (and hence "mpfrMatrix") which nicely 11 and correctly *jointly* formats (for "mpfr") and aligns columns ! 12 Then, formatDec() would be unnecessary. drop0trailing is not really sensible there. 13** TODO 19) outer() now works always ? {as rep() now S3 dispatches 14 ---> need systematic checks *AND* docu changes! 15** TODO 17b) see 'Ops' in R/Arith.R , Rmpfr:::.Math.codes, and 16 design a "test all Ops" with all combinations of "mpfr", "numeric","logical" 17 (and possibly more). 18 19 20* Short or Mid Term 21*** TODO Split this section into Short | Mid (?) 22** TODO 60) Should have *exact* as.bigq.mpfr(), i.e, "mpfr" --> "bigq". (R's "bigq" is C 'mpq') 23** TODO 62) integrateR(): option 'all.sums=TRUE' --> R/integrate-Romberg.R 24*** Inside the GMP library, have 25# -- Function: void mpq_set_f (mpq_t ROP, const mpf_t OP) 26# Set ROP to the value of OP. There is no rounding, this conversion is exact. 27*** MPFR documents 'mpf2mpfr.h' after which you can compile any mpf_ program.. 28** DONE 63) R interface to mpfr's functions `mpfr_get_ld_2exp()` and `mpfr_frexp()` 29 compatibly to DPQ's `ldexp(f, E)` and `frexp(x)`: -> frexpMpfr() and ldexpMpfr(). 30** TODO 53) plogis() {and dlogis, qlogis} are "easy": do use <R>/src/nmath/[dpq]logis.c, as they 31 already use all the numerical tricks including R_Log1_Exp(x) .. : 32 R_Log1_Exp(x) := ((x) > -M_LN2 ? log(-expm1(x)) : log1p(-exp(x))) 33** TODO 54) zapsmall() would be nice, base::zapsmall() fails 34** DONE 59) *exact* dhyper(), phyper().. -> in package 'DPQmpfr' dhyperQ() etc exact via 'gmp' exact rationals 35** TODO 35) tests/bit-repr.R : Bits() fails to work with 2^k 36 37** TODO 37) mpfrXport()/*Import() should work for arrays. Test Windows-portability: 38 - as save() seems not portable; see ~/R/MM/Pkg-ex/Rmpfr/save-load-ex.R 39** TODO 2) Now have working "mpfrMatrix", dim(.) <- ..; t(), %*%, crossprod()... 40 - <mpfr> %*% <mpfr> should work the same as with numeric vectors 41 - <mpfr> %*% t(<mpfr>) ditto 42 43 Note that matrix multiplication seems too slow --> ./Savicky-matrix-mult_tst.R 44 45 - <mpfr>[i] & <mpfrMatrix>[i] work 46 but <mpfrMatrix> [i,j] not yet 47 48 --> want things to work like 49 which(<mpfrMatrix> == ., arr.ind = TRUE) - ok 50 51 [No longer sure if this is true :] 52 For this, we must ensure that the methods are used, instead of the 53 .Primitive base functions : 54 One way: --> see ~/R/MM/NUMERICS/bessel-large-x.R 55 -------------------------------- 56 ## really interesting is bI(x., nu) {for "mpfr" argument}: 57 ## it uses outer(), but that needs to dispatch on, e.g. "^", 58 ## i.e., not only look at "base" 59 environment(outer) <- as.environment("package:Rmpfr") 60 environment(dim) <- as.environment("package:Rmpfr") 61 environment(dimnames) <- as.environment("package:Rmpfr") 62 environment(`dim<-`) <- as.environment("package:Rmpfr") 63 environment(`dimnames<-`) <- as.environment("package:Rmpfr") 64 environment(which) <- as.environment("package:Rmpfr") 65 66** TODO 5) have seqMpfr(), but would like seq() methods, but that seems currently impossible because of 67 a "design infelicity" in base::seq.default --- ???? E-mail to R-core ?? 68 --> R/mpfr.R 69** TODO 6) It is "wrong" that the class "Mnumber" also extends "character", "list"; 70 but it's not clear we can find better definitions, see R/AllClasses.R 71 72** TODO 7) Add tests for hypot() & atan2() to tests/special-fun-ex.R 73 74** TODO 8) round(x, .) & signif(x, .) currently return "mpfr" numbers of the same precision. 75 That *looks* ugly. 76 Potentially add a swith 'keepPrec = FALSE' -- i.e. by default *reduce* 77 precision to "match" 'digits' argument. 78 79** TODO 16) psigamma(x, n) {and digamma(), trigamma() aliases} 80 --> experiments in ~/R/MM/Pkg-ex/Rmpfr/psigamma.R ) 81 Note that since, MPFR 3.0.0, there is a digamma(); .. which we now interface to 82 83** TODO 18) ifelse() fails ... maybe I should mask it 84 {or "fix it" via assign*() in base ?? -- they will love that!} 85 or provide ifelse2() -- a fast simplified version ? 86 87 88** TODO 24) Bernoulli(): we use builtin zeta(); alternatively, use *exact* 89 rationals from 'gmp', using "bigq" (and "bigz") -- and R code from ~/R/Pkgs/copula/R/special-func.R 90** TODO 26) (?) Revert the decision to *not* care about rounding mode in Ops/function, 91 and rather expose that as in mpfr(), e.g., in R/Math.R roundMpfr 92 -- see MPFR_RNDN in src/utils.c and others; --> use src/convert.c R_rnd2MP() 93 -- and all the SEXP functions get an additional SEXP rnd_mode 94 argument, the same as SEXP d2mpfr1() in src/convert.c has already. 95** TODO 29) Our sum() should use system mpfr_sum() : 96 mpfr_sum (mpfr_t ROP, mpfr_ptr const TAB[], unsigned long int N, mpfr_rnd_t RND) 97** TODO 32) Use ./LUP.R to compute the lu() decomposition of an mpfrMatrix 98 ---> solve() and use this for determinant() for larger (n >= 4 ?) dimensions! 99** TODO 50) For *complex* arithmetic, build interface to the "MPC" library 100 ---> http://www.multiprecision.org/mpc -- which is LGPL and itself 101 builds on MPFR and GMP. 102 Ubuntu now has 'libmpc-dev' (!) {but there's no '*-doc' package yet; 103 on nb-mm, I've installed from source 104 --> Info comes along} 105 One idea: Since the names are so much in parallel, try 106 to take each src/*.c file and mechanically s/mpfr/mpc/ producing a 107 mpc version of that; 108 then "the same" for the R interface code. 109 110** TODO 51) Incomplete gamma, i.e. pgamma(), is being added to MPFR. -> do in Rmpfr! 111 ~/F/IN-lists--2016-08 : 112 From: paul zimmermann <Paul.Zimmermann@inria.fr> 113 Subject: [MPFR] incomplete Gamma function 114 Date: Mon, 18 Jan 2016 09:54:51 +0100 115 116 zimmerma@tomate:~/mpfr/tests$ ./tgamma_inc 60 30 100 117 1.3868299023788801161747839921242e80 118 119*** Will be in MPFR 3.2.x (Oct.2016: current is 3.1.5) 120*** MPFR-devel(svn): https://gforge.inria.fr/scm/viewvc.php/mpfr/trunk/src/gamma_inc.c?view=markup 121 - get: --> http://mpfr.loria.fr/gforge.html explains: 122 - svn checkout https://scm.gforge.inria.fr/anonscm/svn/mpfr/trunk mpfr 123 124 125* Accomplished 126** DONE 1) R: character -> mpfr "3.14159265358979323846264" -> mpfr 127 128** DONE 3) "Arith" and "Compare" methods currently ``lose dim + dimnames'' 129 for "mpfrArray" (& "mpfrMatrix") 130 131 The solution is a bit tedious because the Ops do recycling pretty 132 generously for vectors, 133 but pretty stringently when one of the operands is a matrix. 134 135 If the other part is a matrix their dim() must be identical, 136 if just a vector, its length must be a divisor of length(<matrix>) 137 138** DONE 10b) a factorialMPFR() which automatically uses full precision for 139 integer-valued argument, notably using MPFR's mpfr_fac_ui; see also end 140 of man/mpfr-class.Rd 141 142** DONE 13) all the NOT_YET in src/Ops.c are implemented, *apart* from trigamma() 143 --> TODO 16) 144 145** DONE 14) Want to *change* 'precBits' of existing MPFR numbers; 146 MPFR has mpfr_set_prec(X, PREC) but that sets the value to NaN. 147 Manual: "In case you want to keep the previous value stored in X, use 148 `mpfr_prec_round' instead." 149 150 --> fulfilled via roundMpfr(x, precBits) 151 152 153 154** DONE 15) beta(.,.) and lbeta(.,.) .. using my_mpfr_beta() in C. 155 Interestingly, the speedup is not dramatical 156 (50% for length 200; 300% for length 1) 157 158** DONE 4) format() got more (optional) arguments, along the format.default() 159 example. 160 Note that an option to "round() after decimal" should not be needed, 161 rather format(round(..), digits= ., drop0trailing=TRUE) does work. 162 163 164** DONE 12) crossprod(), tcrossprod() (and more?) methods for "mpfrMatrix". 165 166** DONE 10) chooseMpfr(a,n) is now implemented --- *NOT* based on gamma(), 167 but rather n. 168 169** DONE 11b) No longer --- problem were missing mpfr_clear() statements in src/utils.c : 170 format(<mpfr>) --> .mpfr2str() -> C mpfr2str() still suffers from a 171 memory bug, inspite of my efforts in src/convert.c 172 I think this is the MPFR library just allocating memory that's in use 173 by R, but it seems hard to prove and or fix that. 174 175 176** DONE 17a) as(1, "mpfr") & TRUE : no longer gives infinite recursion error 177 178** DONE 17a) Write a 'Rmpfr-package' help page that mentions that we have *many* 179 Math etc functions, including gamma, digamma, .... {which are not 180 quickly visible on the help pages now}. 181 182 183** DONE 20) integrateR( ... rel.tol, verbose= TRUE) : 184 the precision of the output should be increased a bit, 185 (still depending on rel.tol !) 186 187** DONE 21) What's the exponent range -- and possibly change it: R interface 188 in R/mpfr.R via .mpfr_erange() [and .mpfr_erange(.) <- v ] 189 to 190 - Function: mpfr_exp_t mpfr_get_emin (void) 191 - Function: mpfr_exp_t mpfr_get_emax (void) 192 193 - Function: int mpfr_set_emin (mpfr_exp_t exp) 194 - Function: int mpfr_set_emax (mpfr_exp_t exp) 195 196 - Function: mpfr_exp_t mpfr_get_emin_min (void) 197 - Function: mpfr_exp_t mpfr_get_emin_max (void) 198 - Function: mpfr_exp_t mpfr_get_emax_min (void) 199 - Function: mpfr_exp_t mpfr_get_emax_max (void) 200 201** DONE 22) apply(<mpfrArray>, .) --> S4 method 202 203** DONE 22x) *do* test rank(), sort(), order() : we claim they work in man/mpfr-class.Rd 204 quantile() should work from R 2.15.1 (June 2012) on ~~~~~~~~~~~~~~~~~ 205 206** DONE 23) quantile(<mpfr>) does not work ---- but will from R 2.15.1 with better quantile.default() 207 Reason: in stats::quantile.default(), 208 ifelse(h == 0, qs[i], (1 - h) * qs[i] + h * x[hi[i]]) 209 produces a list of 'mpfr1' instead a real 'mpfr' vector. 210 -> Fixed in ~/R/D/r-devel/R/src/library/stats/R/quantile.R 211 212** DONE 23x) sumBinomMpfr() accept f.x 213 214** DONE 28) determinant() and hence det() "fail" (with a reasonable error message). 215 Easy: go via asNumeric(.) with a warning() that can be suppressed 216 Standard R uses LAPACK's LU() decomposition. Would be cool to have 217 that in C (using MPFR). 218 Alternatively, do it for 2x2 and via recursion mxm for small m. 219 (but it is *really* inefficient: complexity \propto m! ) 220 221** DONE 30) pmin() and pmax(), for simple cases pmin(x,n) are now, 2013-02, quite "fast", 222 (see also --> ~/R/MM/Pkg-ex/Rmpfr/SLOW-pmin.R) 223** DONE 33) asNumeric(mpfr(10, 99)^500) gives 1.797693e+308 without warning. Should (warn or) give Inf 224** DONE 36) got mpfrImport(), reading mpfrXport() 225** DONE 38) Define a norm() method. For expm(*, "Higham08") need solve() ! 226** DONE 34) Can specify the rounding mode in mpfr() and roundMpfr(), since 2015-08-08 227** DONE 9) median(<mpfr>) now works, as does mean(<mpfr>, trim = *) and quantile() w/o warning. 228** DONE is.finite(.) etc: *Must* keep matrices -> src/utils.c 229** DONE 52) No longer use 'representation' but rather 'slots' in 'setClass(..)' 230** DONE 25) mpfr(<mpfr>, ...) should work: return argument unchanged *if* precision 231 and rounding mode are ok; otherwise use roundMpfr() *BUT* that needs an 232 additional rnd.mode argument -- as mpfr 233** DONE 27) Now have simple sapplyMpfr() 234** DONE 57) Fix pnorm() bugs reported by Jerry Lewis (>> ~/R/MM/Pkg-ex/Rmpfr/pnorm-* ) 235** DONE 58) print(<large Matrix>) somehow fails to use max.digits = 9999 [print.mpfr(*, max.digits = ...)] ??? 236*** Why did we introduce 'max.digits' *and* allow digits=NULL to use *all* digits before the decimal point? 237**** I now think max.digits is almost superfluous if digits=NULL behaved as promised: 238 We have in R code comments __digits = NULL : use as many digits "as needed"__ , 239 and in help man/formatMpfr.Rd 240 __The default, \code{NULL}, uses enough digits to represent the full precision, often one or two digits more than you would expect.__ 241** DONE 56) Add regr.tests for new dgamma(<mpfr>) !! 242** DONE 55a) .getPrec(), .mpfr_erange, .mpfr_erange_set, .mpfr_maxPrec, .mpfr_minPrec etc are documented now. 243 244