search for: dbl_epsilon

...= 1) gave 1000000 intervals (on a machine), made me trace through the code to function 'R_pretty' in https://svn.r-project.org/R/trunk/src/appl/pretty.c . *lo is -1e300, *up is 1e300. cell = fmax2(fabs(*lo),fabs(*up)); 'cell' is 1e300. i_small = dx < cell * U * imax2(1,*ndiv) * DBL_EPSILON *3; When *ndiv is (int) 1e9, apparently cell * U * imax2(1,*ndiv) overflows to infinity and 'i_small' is 1 (true). It doesn't happen when *ndiv is (int) 1e6. Putting parentheses may avoid the floating point overflow. For example, i_small = dx < cell * (U * imax2(1,*ndiv) * DBL_EPSIL...

...de to function 'R_pretty' in https://svn.r-project.org/R/trunk/src/appl/pretty.c . ? ? > thank you. ? ? >> *lo is -1e300, *up is 1e300. ? ? >> cell = fmax2(fabs(*lo),fabs(*up)); ? ? >> 'cell' is 1e300. ? ? >> i_small = dx < cell * U * imax2(1,*ndiv) * DBL_EPSILON *3; ? ? >> When *ndiv is (int) 1e9, apparently cell * U * imax2(1,*ndiv) overflows to infinity and 'i_small' is 1 (true). It doesn't happen when *ndiv is (int) 1e6. [[elided Yahoo spam]] ? ? >> Putting parentheses may avoid the floating point overflow. For example, ? ? &gt...

...1.0D+16 [...] RTNSIG ENMTEN XLARGE [...] IEEE (IBM/XT, SUN, etc.) (D.P.) 1.0E-4 8.90D-308 1.0D+4 [...] Is there some reason for using different values in bessel_j.c than in bessel_l.c? Or is this just a bug? ;^) -Greg PS: Where is DBL_EPSILON defined, and what is its value? -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the su...

...e me trace through the code to function 'R_pretty' in https://svn.r-project.org/R/trunk/src/appl/pretty.c . thank you. > *lo is -1e300, *up is 1e300. > cell = fmax2(fabs(*lo),fabs(*up)); > 'cell' is 1e300. > i_small = dx < cell * U * imax2(1,*ndiv) * DBL_EPSILON *3; > When *ndiv is (int) 1e9, apparently cell * U * imax2(1,*ndiv) overflows to infinity and 'i_small' is 1 (true). It doesn't happen when *ndiv is (int) 1e6. well spotted! > Putting parentheses may avoid the floating point overflow. For example, > i_small = dx &...

...tiply w by the factor a b _ _ _ x (1-x) | (a+b) / ( a | (a) | (b) ) . */ y = a * log(x); t = b * log(xc); /* Resort to logarithms. */ y += t + lgamma(a + b) - lgamma(a) - lgamma(b); y += log(w / a); t = exp(y); done: if (flag == 1) { if (t <= DBL_EPSILON) t = 1 - DBL_EPSILON; else t = 1 - t; } return t; } /* Continued fraction expansion #1 * for incomplete beta integral */ static double incbcf(double a, double b, double x) { double xk, pk, pkm1, pkm2, qk, qkm1, qkm2; double k1, k2, k3, k4, k5, k6, k7, k8; double r, t, ans, thresh;...

...s the original NaN came from division of sinh and cosh that had reached machine infinity. static void z_tan(Rcomplex *r, Rcomplex *z) { double x2, y2, den; x2 = 2.0 * z->r; y2 = 2.0 * z->i; den = cos(x2) + cosh(y2); r->r = sin(x2)/den; /* any limit above about log(DBL_EPSILON) will do */ if (fabs(r->i) < 40.0) r->i = sinh(y2)/den; else r->i = 1.0; }

...> for (i = 0, j = 0; i < len; i++, j += size) { > ld1 = (long double) REAL(object)[i]; > format.c: > static const long double tbl[] = > ... but it doesn't. Perhaps the original code at issue: > arithmetic.c: > static LDOUBLE q_1_eps = 1 / LDBL_EPSILON; > only makes gcc unhappy because of the very large value trying to be stored > in the static long double, which would make it span the "folded double" on > that architecture. > ***** > It seems that the options are: > A) Patch the one place w...

...pecs=/usr/lib/rpm/redhat/redhat-annobin-cc1 -m64 -mcpu=power8 > -mtune=power8 -fasynchronous-unwind-tables -fstack-clash-protection -c > arithmetic.c -o arithmetic.o > arithmetic.c:180:26: error: initializer element is not constant > 180 | static LDOUBLE q_1_eps = (1 / LDBL_EPSILON); > | ^ > make[3]: *** [../../Makeconf:124: arithmetic.o] Error 1 > Took me a bit to figure out why, but this is happening because on > powerpc64, gcc is compiled with -mlong-double-128, and the long double > format used on PPC is IBM'...

Full_Name: Morten Welinder Version: 2 OS: Solaris/space/gcc2.95.2 Submission from: (NULL) (65.213.85.217) I changed src/nmath/standalone/test.c to read: --------------------------------------------------------------------------------- #define MATHLIB_STANDALONE 1 #include <Rmath.h> #include <stdio.h> int main() { double x; for (x = 99990; x <= 100009; x++) printf

...; i++, j += size) { >>> ld1 = (long double) REAL(object)[i]; > >>> format.c: >>> static const long double tbl[] = > >>> ... but it doesn't. Perhaps the original code at issue: > >>> arithmetic.c: >>> static LDOUBLE q_1_eps = 1 / LDBL_EPSILON; > >>> only makes gcc unhappy because of the very large value trying to be stored >>> in the static long double, which would make it span the "folded double" on >>> that architecture. > >>> ***** > >>> It seems that the options are:...

...d Nocedal, the 2nd one being one of the original L-BFGS-B(1997) paper, make two remarks, the 2nd one about the "machine epsilon" used, and I can assure you that R's optim() version never suffered from that; we've always been using a C translation of the fortran code, and then used DBL_EPSILON. R's (main) source file for that is in .../src/appl/lbfgsb.c, e.g., here https://svn.r-project.org/R/trunk/src/appl/lbfgsb.c OTOH, their remark 1 is very relevant and promising faster / more reliable convergence. I'd be "happy" if optim() could gain a new option, say, "L-BF...

.../rpm/redhat/redhat-hardened-cc1 -specs=/usr/lib/rpm/redhat/redhat-annobin-cc1 -m64 -mcpu=power8 -mtune=power8 -fasynchronous-unwind-tables -fstack-clash-protection -c arithmetic.c -o arithmetic.o arithmetic.c:180:26: error: initializer element is not constant 180 | static LDOUBLE q_1_eps = (1 / LDBL_EPSILON); | ^ make[3]: *** [../../Makeconf:124: arithmetic.o] Error 1 Took me a bit to figure out why, but this is happening because on powerpc64, gcc is compiled with -mlong-double-128, and the long double format used on PPC is IBM's 128bit long double which is two doub...

...1.0D+16 [...] RTNSIG ENMTEN XLARGE [...] IEEE (IBM/XT, SUN, etc.) (D.P.) 1.0E-4 8.90D-308 1.0D+4 [...] Is there some reason for using different values in bessel_j.c than in bessel_l.c? Or is this just a bug? ;^) -Greg PS: Where is DBL_EPSILON defined, and what is its value? -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the sub...

...om.c, qnbinom.c -- * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ if(!lower_tail || log_p) { p = R_DT_qIv(p); /* need check again (cancellation!): */ if (p == 0.) return 0; if (p == 1.) return ML_POSINF; } /* temporary hack --- FIXME --- */ if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF; /* y := approx.value (Cornish-Fisher expansion) : */ z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);

...nity. Why not use an exit condition based on how far away from 'a' you are? qbeta.c:174 (prevent excess precision) Might this not just get optimized out? I doubt modern compilers respect the volatile keyword. You should probably replace this with a safe comparison. if(fabs(tx-xinbta)<=DBL_EPSILON*fabs(xinbta))... ** This should be checked with someone familiar with floating point arithmetic; I'm not and may have got the correct expression wrong. The 'infinite loop' is due to the speed of the pbeta routine, rather than the initial approximation from the qbeta. There is another a...

...= 0; i < len; i++, j += size) { >> ld1 = (long double) REAL(object)[i]; >> format.c: >> static const long double tbl[] = >> ... but it doesn't. Perhaps the original code at issue: >> arithmetic.c: >> static LDOUBLE q_1_eps = 1 / LDBL_EPSILON; >> only makes gcc unhappy because of the very large value trying to be stored >> in the static long double, which would make it span the "folded double" on >> that architecture. >> ***** >> It seems that the options are: >>...

...pile times go down from 0.36s to 0.13s on x86 and 0.35s to 0.19s on x86-64 as expected. Here's the code: #include <float.h> #include <stdio.h> #include <stdlib.h> #include <math.h> #define real float #define epsilon FLT_EPSILON //#define real double //#define epsilon DBL_EPSILON real delta, INFINITY; typedef struct { real x, y, z; } Vec; Vec vec(real x, real y, real z) { Vec r; r.x=x; r.y=y; r.z=z; return r; } Vec add(const Vec a, const Vec b) { return vec(a.x+b.x, a.y+b.y, a.z+b.z); } Vec sub(const Vec a, const Vec b) { return vec(a.x-b.x, a.y-b.y, a.z-b.z); } Vec scale...

...ty. Why not use an exit condition based on how far away from 'a' you are? qbeta.c:174 (prevent excess precision) Might this not just get optimized out? I doubt modern compilers respect the volatile keyword. You should probably replace this with a safe comparison. if(fabs(tx-xinbta)<=DBL_EPSILON*fabs(xinbta))... ** This should be checked with someone familiar with floating point arithmetic; I'm not and may have got the correct expression wrong. The 'infinite loop' is due to the speed of the pbeta routine, rather than the initial approximation from the qbeta. There is anothe...

FYI. http://leonardo-m.livejournal.com/73732.html If anyone is motivated, please file bugs for the losing cases. Also, it might make sense to incorporate the tests into our nightly tester test suite. Thanks, Evan

Hi, Mark et al.: Thanks, Mark. Three comments: 1. Rvmmin was one of the methods I tried after Ravi directed me to optimx. It returned NAs for essentially everything. See my email of this subject stamped 4:43 PM Central time = 21:43 UTC. 2. It would be interesting to know if the current algorithm behind optim and optimx with