Maths constants
in [Fortran]
|
From: Richard Maine on 19 Jan 2010 15:31 Paul Hirose <jvcmz89uwf(a)earINVALIDthlink.net> wrote: > Nowadays even a budget Windows machine does floating point in hardware > per the IEEE standard, so I believe arc cosine of -1 is a safe way to > compute pi, for example. However, I'd be interested to hear of any > cases where this would have caused trouble. Ouch! Ouch! Ouch! Hey, Fortran is often used for numeric stuff and us folk are supposed to have at least a surface acquaintance with issues of numeric accuracy. 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 use literals. -- Richard Maine | Good judgment comes from experience; email: last name at domain . net | experience comes from bad judgment. domain: summertriangle | -- Mark Twain
From: glen herrmannsfeldt on 19 Jan 2010 17:44 Paul Hirose <jvcmz89uwf(a)earinvalidthlink.net> wrote: (snip) > I calculate constants like those at run time. That way, someone > maintaining the code doesn't have to wonder where I got the values, > whether they contain a typo, etc. > Nowadays even a budget Windows machine does floating point in hardware > per the IEEE standard, so I believe arc cosine of -1 is a safe way to > compute pi, for example. However, I'd be interested to hear of any > cases where this would have caused trouble. x87 has the FLDPI instruction to load the internal value of PI in extended real (64 bit significand) into a register. It would be nice to have a way to do that from Fortran. Software that I know of computes: * ARCSIN-ARCCOSINE FUNCTION (SHORT) * 1. IF X BETWEEN 0 AND 1/2, COMPUTE ARCSIN BY POLYNOMIAL * 2. IF X BETWEEN 1/2 AND 1, * ARSIN(X) = PI/2-2*ARSIN(SQRT((1-X)/2)) * 3. IF X NEGATIVE, ARSIN(X) = -ARSIN(/X/) * 4. ARCOS(X) = PI/2-ARSIN(X) So ACOS(-1.) = PI/2-ASIN(-1.) = PI/2+ASIN(1.) = PI/2+PI/2-2*ASIN(0.) I would expect the polynomial used to give 0. for ASIN(0.), but I don't know that it guaranteed. I find no official documentation on the algorithms used by the intel hardware. It could be similar argument reduction and polynomials in microcode, or maybe CORDIC internally in hardware or microcode. > Sometimes the representation of a type is not specified by the > language. It's up to the implementer. With these languages, computing > math constants at runtime automatically upgrades their accuracy if, > say, "double precision" someday changes to a 128 bit representation. > Simply recompile your old code under the new architecture. Generating the appropriately rounded value at run time often requires a little more precision than the final result. That isn't easy to do in a system independent way. -- glen
From: robin on 19 Jan 2010 19:15 "James Van Buskirk" <not_valid(a)comcast.net> wrote in message news:hj2294$kmj$1(a)news.eternal-september.org... | Actually, f03 makes initialization expressions much more of a free- | fire zone. In f95 there was a separate class of constant expressions | which lacked self-consistency and required the vendor to provide 4 | effective compilers (one for ordinary expressions, one for | specification expressions, one for initialization expressions, and | one for constant expressions) Four expression analyzers, not four compilers, three being subsets of the most gereral one. Shouldn't be difficult to implement. | but the first corregendum to f95 | eliminated consant expressions by editing out their effects | throughout the standard. Later, an f95 constant expression was made | to be an f03 initialization expression, and transcendental functions | and transformational functions became fair game at this point. | There were some rough edges involved in these changes and in fact | as a result you can write code that should behave differently under | the f95 and f03 standards, but the differences are esoteric enough | that I don't think any compiler actually catches them when the f95 | standard is enforced. | | gfortran is the most agile cross compiler out there and it does use | a gpl software package to compute initialization expressions. For | it at least the initialization expression results are often more | accurate than the ordinary expression results: | | C:\gfortran\clf\accurate_init>type accurate_init.f90 | module funcs | implicit none | integer, parameter :: dp = selected_real_kind(15,300) | contains | subroutine sub(x) | real(dp) x | real(dp) y | | y = cos(cmplx(atan(real(1,kind(x))),x,kind(x))) | write(*,*) y | end subroutine sub | end module funcs | | program accurate_init | use funcs | implicit none | real(dp), parameter :: x = & | cos(cmplx(atan(real(1,kind(x))),log(huge(x))+1,kind(x))) | real(dp) y | | write(*,*) x | y = log(huge(x))+1 | call sub(y) | end program accurate_init | | C:\gfortran\clf\accurate_init>gfortran accurate_init.f90 -oaccurate_init | | C:\gfortran\clf\accurate_init>accurate_init | 1.72768693203654562E+308 | +Infinity | | This is one where it's tricky to avoid overflow and the initialization | expression package is powerful enough (and maybe slow enough) to get | the right result. Maybe, but it won't compile in F95.
From: Tim Prince on 19 Jan 2010 19:48 glen herrmannsfeldt wrote: > > x87 has the FLDPI instruction to load the internal value of PI > in extended real (64 bit significand) into a register. It would > be nice to have a way to do that from Fortran. > > Software that I know of computes: > > * ARCSIN-ARCCOSINE FUNCTION (SHORT) > * 1. IF X BETWEEN 0 AND 1/2, COMPUTE ARCSIN BY POLYNOMIAL > * 2. IF X BETWEEN 1/2 AND 1, > * ARSIN(X) = PI/2-2*ARSIN(SQRT((1-X)/2)) > * 3. IF X NEGATIVE, ARSIN(X) = -ARSIN(/X/) > * 4. ARCOS(X) = PI/2-ARSIN(X) > > So ACOS(-1.) = PI/2-ASIN(-1.) = PI/2+ASIN(1.) = PI/2+PI/2-2*ASIN(0.) > > I would expect the polynomial used to give 0. for ASIN(0.), but > I don't know that it guaranteed. > > I find no official documentation on the algorithms used by the > intel hardware. It could be similar argument reduction and polynomials > in microcode, or maybe CORDIC internally in hardware or microcode. Among the problems here is that x87 hasn't been a usual mode of operation for Fortran since SSE2 was introduced 8 years ago. If you use an inverse trig function to get Pi, you're at the mercy of the writer of the math library. I studied math back in the days when calculus was considered important to a much larger number of people than numerical analysis. From an error sensitivity point of view, acos(-1.) does seem an egregiously bad choice. Yet, you may luck into a situation such as the above where it gives you the library writer's version of Pi, while, in a substandard imlementation, atan(1.) might give you a polynomial evaluation with roundoff error. > >> Sometimes the representation of a type is not specified by the >> language. It's up to the implementer. With these languages, computing >> math constants at runtime automatically upgrades their accuracy if, >> say, "double precision" someday changes to a 128 bit representation. >> Simply recompile your old code under the new architecture. > > Generating the appropriately rounded value at run time often requires > a little more precision than the final result. That isn't easy to > do in a system independent way. > Other than by putting sufficient digits in your constant (and specifying the data type). IEEE754 laid down the law about requiring accurate use of the appropriate number of digits when converting a constant.
From: Aris on 20 Jan 2010 04:49
nospam(a)see.signature (Richard Maine) wrote: > Paul Hirose <jvcmz89uwf(a)earINVALIDthlink.net> wrote: > >> Nowadays even a budget Windows machine does floating point in hardware >> per the IEEE standard, so I believe arc cosine of -1 is a safe way to >> compute pi, for example. However, I'd be interested to hear of any >> cases where this would have caused trouble. > > Ouch! Ouch! Ouch! > > Hey, Fortran is often used for numeric stuff and us folk are supposed to > have at least a surface acquaintance with issues of numeric accuracy. > 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. Right at the singularity ;-) > 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 use literals. If I would have to implement the arccos, I would return the literal for that input. |