Maths constants
in [Fortran]
|
From: Gordon Sande on 22 Jan 2010 18:09 On 2010-01-22 18:15:43 -0400, "Dr Ivan D. Reid" <Ivan.Reid(a)brunel.ac.uk> said: > On Fri, 22 Jan 2010 21:04:01 GMT, Gordon Sande <g.sande(a)worldnet.att.net> > wrote in <2010012217040116807-gsande(a)worldnetattnet>: > >>> I've done similar w/ my old HP-97 (the desktop version w/ printer, card >>> reader, etc.) It all still functions except the rubber rollers for the >>> paper feed are hard and slick so the paper doesn't feed well any longer. > >> You can buy a small bottle of numble-numble-alcohol (isoproponol?) which will >> remove the glaze and soften the rubber rollers. > > Close. Isopropanol, or iso-propyl alcohol. Used as swabs before > injections, etc, etc. I thought you guys called it rubbing alcohol, but my > Chambers suggests rubbing alcohol is methylated spirits (denatured ethanol); > maybe there's a transAtlantic schism there. Methanol used to be called rubbing alcohol in my experience. So there is clearly a problem with common names. Sort of like benzine, benzene, etc! Never took enough organic chemistry to know the difference between the various alcohols and other things that now seem to called surfactants. In any case rubber rollers can be nicely and easily restored with the suitable potion whatever is its correct name.
From: Richard Maine on 22 Jan 2010 18:35 Paul Hirose <jvcmz89uwf(a)earINVALIDthlink.net> wrote: > The only source code for acos() at hand here is P.J. Plauger, "The > Standard C Library", 1992. His code treats -1 as a special case, and > produces an accurate value for pi. For arguments other than special > cases, he uses polynomials from one of the standard reference books Seems to me that the necessity for special-casing illustrates the point. I presume this means that if someone implements it without such special-casing, the value is likely to be inaccurate. No, I haven't run into concrete examples of this one failing. Of course, that's partly because I don't use it, so I wouldn't have run into such failures in any case. However, I've also long learned not to regard "it hasn't failed for me before" as a good reason to use an inherently fragile method. If it requires special-casing to get an accurate answer, I'd call that support of the notion that it is inherently fragile. But the theory pretty much told me that story anyway. Yes, I regard such theoretical basis as important. That's the kind of thing that helps give you confidence that "it hasn't failed for me before" doesn't translate into "the first time it will bite me is tomorrow, which will be a really inconvenient time according to Murphy." The odds of it being a big issue for this particular usage? Probably low, I'd guess. But the odds of a general attitude of ignoring issues of numerical sensitivity will get someone in trouble? Much higher. I have run into such things a lot in helping people understand why their programs are giving garbage results. Well, trying to help anyway - there have been cases where people just could not see, whether because my explanation was inadequate or for whatever other reason. There was the fellow who was sure that the problems he had would be solved by going to quad precision (and wanted a grant to do that conversion). I could see at a glance that the problem was ill posed although it looked like there might be some interesting ideas in his work if he reformulated to avoid that. I couldn't get him to see that. He didn't get the grant. Although I don't use the trig functions methods for pi myself, I thought that the 4*atan(1) method mentioned elsethread was more the "usual" one to use. That is at a numerically "nice" point in the atan curve. -- Richard Maine | Good judgment comes from experience; email: last name at domain . net | experience comes from bad judgment. domain: summertriangle | -- Mark Twain
From: William Clodius on 22 Jan 2010 19:28 Richard Maine <nospam(a)see.signature> wrote: > <snip> > Although I don't use the trig functions methods for pi myself, I thought > that the 4*atan(1) method mentioned elsethread was more the "usual" one > to use. That is at a numerically "nice" point in the atan curve. However for ATAN(10 the implementor of the library for some implementations of floating point (IEEE single?) has to decide whether to obey the constraint atan(x) .le. pi/2, if the closest approximation to pi/2 is biger than pi/2. -- Bill Clodius los the lost and net the pet to email
From: Ron Shepard on 22 Jan 2010 19:54 In article <0_KdnUzqPYbKgcfWnZ2dnUVZ_rmdnZ2d(a)earthlink.com>, "Paul Hirose" <jvcmz89uwf(a)earINVALIDthlink.net> wrote: > "Richard Maine" <nospam(a)see.signature> wrote in message > news:1jckfrv.w6yp2rh7f8q7N%nospam(a)see.signature... > > The arc cosine of -1 is just about the worst case I can think of off > > the > > top of my head in terms of numerical robustness. Hint: look at the > > arc > > cosine of numbers very close to -1. If one is going to use trig > > expressions for things like that, at least use them at decently > > robust > > points. That's independent of the language and the hardware - just > > plain > > numerics. > > I hear you. That's why I formerly used 2*acos(0d0) to obtain pi. [...] This expression has the same problem as acos(-1.0) and 2*asin(1.0) that I discussed previously. The exact derivative of COS(theta) for theta=PI/2 is zero, so it is a second-order function in that neighborhood. Any floating point value of theta in the range PI/2-sqrt(eps) to PI/2+sqrt(eps) will result in COS(theta)=0.0, so the ACOS function might return any one of them. That is a LOT of numbers, right? Of course, the library writer might make a special case of a 0.0 argument and return the right result, but that is not guaranteed by the language standard. Something else that has not been mentioned in this context is that the results of the algorithms used to compute these function values might depend on the rounding mode for floating point arithmetic. In this case, changing the rounding mode setting might result in different values of the intrinsic functions, which would mean that the returned values for PI, or other "constants", might change depending on when they are computed during the programs execution or on other details of the computational environment. This is yet another reason to specify these values with literal constants. $.02 -Ron Shepard
From: steve on 22 Jan 2010 20:19
On Jan 22, 3:35 pm, nos...(a)see.signature (Richard Maine) wrote: > Paul Hirose <jvcmz89...(a)earINVALIDthlink.net> wrote: > > The only source code for acos() at hand here is P.J. Plauger, "The > > Standard C Library", 1992. His code treats -1 as a special case, and > > produces an accurate value for pi. For arguments other than special > > cases, he uses polynomials from one of the standard reference books > > Seems to me that the necessity for special-casing illustrates the point. > I presume this means that if someone implements it without such > special-casing, the value is likely to be inaccurate. > No. It means that it is faster to check for x = -1 than to use an approximation. http://svn.freebsd.org/viewvc/base/head/lib/msun/src/e_acosf.c?view=markup I was hoping to stay out of this history driven thread. 30 or 40 years one may have had issues with acos(-1.0), but troutmask:sgk[210] make acos cc -o acos -O2 -g -pipe -fno-builtin -march=native -static -I/usr/ local/include -I../mp acos.c -L/usr/local/lib -L../mp -lsgk -lmpfr - lgmp -lm_p sleep 1 touch acos.c troutmask:sgk[211] ./acos Max ULP: 0.513475 Check for line wrapping if you want to run the following. It checks the max err in ULPs for acos(-1.0) and its closest 10000 neighbors. #include <stdio.h> #include <float.h> #include <math.h> #include <limits.h> #include <gmp.h> #include <mpfr.h> #include "sgk.h" #if 1 #define DIGITS FLT_MANT_DIG #define ACOS acosf #define EPS FLT_EPSILON #define NEXT nextafterf #else #define DIGITS DBL_MANT_DIG #define ACOS acos #define EPS DBL_EPSILON #define NEXT nextafter #endif int main(int argc, char *argv[]) { #if 1 float x, y, z; float u, um = 0.; #else double x, y, z; double u, um = 0.; #endif mpfr_t xacc, xapp; mpfr_init2(xacc, 128); mpfr_init2(xapp, DIGITS); x = -1.; y = -1. + 10000 * EPS; do { z = ACOS(x); mpfr_set_d(xapp, (double)z, RD); mpfr_set_d(xacc, x, RD); mpfr_acos(xacc, xacc, RD); u = mp_ulp(xapp, xacc, DIGITS); if (u > um) um = u; x = NEXT(x, 1.); } while (x < y); fprintf(stderr, "Max ULP: %f\n", um); mpfr_clears(xacc, xapp, NULL); return 0; } -- steve |