On Saturday, July 2, 2022 at 10:37:38 PM UTC+2, Krishna Myneni wrote:
The test code appears to show this version,
running under kForth, exceeds the accuracy reported in the paper for
double precision calculations. One should note that kForth's FSIN and
FCOS make use of the C math library routines for computing sin(x) and
cos(x), instead of the f.p.u. processor instruction FSINCOS -- the
latter has a well-known degradation of performance for large arguments,
x, while the C math library functions modulo the argument into a range
with full accuracy.
Thanks for publishing this. It was a breeze to port to 4tH and ran straight away. That having said, the accuracy I achieved was MUCH lower - and probably due to the Remez algorithm of my high level FSIN + FCOS functions.
But if these ranges are that large, why not use your own "fsin near pi" function?
: fsin_near_npi ( rdelta n -- r )
2 mod if fnegate then
fdup 7 fpow 5040 s>f f/
fover 5 fpow 120 s>f f/ f-
fover 3 fpow 6 s>f f/ f+
f-
;
The test code appears to show this version,
running under kForth, exceeds the accuracy reported in the paper for
double precision calculations. One should note that kForth's FSIN and
FCOS make use of the C math library routines for computing sin(x) and
cos(x), instead of the f.p.u. processor instruction FSINCOS -- the
latter has a well-known degradation of performance for large arguments,
x, while the C math library functions modulo the argument into a range
with full accuracy.
The range of the argument to the different j_l(x) and n_l(x) does notI dunno - YOU tell ME ;-)
exceed 300 in the test code, which is likely too small to see a
difference between using the native FSINCOS fpu instruction vs the C
library FSIN and FCOS functions. A comparison against values for large x should show the problem.
I'm not familiar with the algorithm above. Doesn't range reduction
require higher precision arithmetic?
\ Krishna Myneni, <krishn...@ccreweb.org> <<===
On Sunday, July 3, 2022 at 1:40:18 PM UTC+2, Krishna Myneni wrote:
The range of the argument to the different j_l(x) and n_l(x) does notI dunno - YOU tell ME ;-)
exceed 300 in the test code, which is likely too small to see a
difference between using the native FSINCOS fpu instruction vs the C
library FSIN and FCOS functions. A comparison against values for large x
should show the problem.
I'm not familiar with the algorithm above. Doesn't range reduction
require higher precision arithmetic?
HB
\ High accuracy calculation of double-precision sin(x),
\ using Taylor Series approximation when x is near
\ a multiple of pi.
\ Due to finite precision in the representation of x,
\ the calculation of sin(x) can suffer a loss of
\ significant digits when x is near pi.
\ The Intel and AMD FPUs also suffer a loss of accuracy in
\ the FSIN instruction, for an argument, x, which is
\ near a multiple of pi, n*pi, for n>1. This function uses
\ an 10th order Taylor Series approximation to compute sin(x)
\ for those cases:
\ sin(pi + delta) = -delta + delta^3/6 - delta^5/120
\ + delta^7/5040 + ...
\ Higher order terms are not used.
\ For |delta| <= 1e-5, the accuracy in double precision
\ will be 16 significant digits.
\ Krishna Myneni, <krishn...@ccreweb.org> <<===
\ Compute sin(n*pi + rdelta) to high accuracy, where
\ rdelta is an offset and n = 0,1,2, ...
[UNDEFINED] fsin_near_pi [IF]
[UNDEFINED] fpow [IF] include lib/fpow.4th [THEN]
: fsin_near_npi ( rdelta n -- r )
2 mod if fnegate then
fdup 7 fpow 5040 s>f f/
fover 5 fpow 120 s>f f/ f-
fover 3 fpow 6 s>f f/ f+
f-
;
[THEN]
The above can be considered a special case of range reduction only whenOh! I didn't know that! I thought it was +/-0.5 pi ;-)
the argument is very close to pi for sin(x) (delta <= 1e-5).
Range reduction of a general argument into the interval -pi,piI know all about that. And it gets worse when precision is pretty low.
requires higher precision to avoid degrading the precision of the
original argument.
On Sunday, July 3, 2022 at 4:37:20 PM UTC+2, Krishna Myneni wrote:
The above can be considered a special case of range reduction only whenOh! I didn't know that! I thought it was +/-0.5 pi ;-)
the argument is very close to pi for sin(x) (delta <= 1e-5).
Range reduction of a general argument into the interval -pi,piI know all about that. And it gets worse when precision is pretty low.
requires higher precision to avoid degrading the precision of the
original argument.
So, SIN(1000 PI) is always a problem (when given as an FP number).
But many thanks for your clarifications! I learned something!
The above results may be compared with high precision reference values
for sin(x), e.g. using maxima, Wolfram Alpha, etc.
It occurred to me that we should be able to implement proper range
reduction in Forth using Julian Noble's double-double arithmetic words,
for cases where it is important, without resorting to C math library >functions. I have this issue in kForth-Win32 -- the Digital Mars C++ >compiler's math library does not do proper range reduction for the trig >functions.
--
Krishna
In article <t9uli5$3ces2$1@dont-email.me>,
Krishna Myneni <krishna.myneni@ccreweb.org> wrote:
<SNIP>
The above results may be compared with high precision reference values
for sin(x), e.g. using maxima, Wolfram Alpha, etc.
It occurred to me that we should be able to implement proper range
reduction in Forth using Julian Noble's double-double arithmetic words,
for cases where it is important, without resorting to C math library
functions. I have this issue in kForth-Win32 -- the Digital Mars C++
compiler's math library does not do proper range reduction for the trig
functions.
I looked at range reduction working on tforth's floating point
packet.
For a sufficient large x the uncertainly in x is +/-pi , and the
sin and cosine may be all over the place. in the range [-1,+1].
How can range reduction help, unless by accident it is known that
x is mathematically exact (which I guess would be rare for physics calculations) ?
On 7/3/22 10:40, Hans Bezemer wrote:
On Sunday, July 3, 2022 at 4:37:20 PM UTC+2, Krishna Myneni wrote:
The above can be considered a special case of range reduction only whenOh! I didn't know that! I thought it was +/-0.5 pi ;-)
the argument is very close to pi for sin(x) (delta <= 1e-5).
Range reduction of a general argument into the interval -pi,piI know all about that. And it gets worse when precision is pretty low.
requires higher precision to avoid degrading the precision of the
original argument.
So, SIN(1000 PI) is always a problem (when given as an FP number).
But many thanks for your clarifications! I learned something!
Some examples of reduced accuracy without range reduction may be
illustrative for others who are not aware of the issue we are discussing.
In kForth FSINCOS standard floating point word uses the native fpu instruction for computing the sin(x) and cos(x) functions -- this
provides a fast way of computing these functions when it is known in
advance that the range of x will be limited. In contrast the standard
FSIN and FCOS words call GCC's math library functions, which use proper
range reduction before calling the native fpu instructions. These words
are useful for obtaining the best accuracy over an unknown range of arguments.
Below is a comparison of using the library version of FSIN vs the fpu instruction for FSIN in kForth:
17 set-precision
ok
\ x = 10,000.
1.0e4 fsin fs.
-3.0561438888825215e-01 ok \ with range reduction
1.0e4 fsincos fdrop fs.
-3.0561438888825215e-01 ok \ native fpu instruction
\ x = 100,000.
1.0e5 fsin fs.
3.5748797972016508e-02 ok \ with range reduction
1.0e5 fsincos fdrop fs.
3.5748797972016383e-02 ok \ native fpu instruction
\ x = 1,000,000.
1.0e6 fsin fs.
-3.4999350217129294e-01 ok \ with range reduction
1.0e6 fsincos fdrop fs.
-3.4999350217129177e-01 ok \ native fpu instruction
The above results may be compared with high precision reference values
for sin(x), e.g. using maxima, Wolfram Alpha, etc.
Forth code for calculation over a large set of spherical Bessel and
spherical Neumann functions, with standard normalization, is available
in the kForth repos. -- see, for example,
https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/sph_bes_neu.4th
While the FSL provides a spherical Bessel function calculator for the
first ten functions, j_0(x) -- j_9(x), the present code extends the calculation range to beyond j_100(x). The spherical Neumann functions,
n_l(x) (also often denoted by y_l(x)) are also computed over the same
range of l.
Test code is included, with reference values obtained from Wolfram
Alpha. The original source for this code is a paper in Computers in
Physics vol. 2, pp. 62--72 (1988), linked below, in which they provide
an accuracy analysis. The test code appears to show this version,
running under kForth, exceeds the accuracy reported in the paper for
double precision calculations. One should note that kForth's FSIN and
FCOS make use of the C math library routines for computing sin(x) and
cos(x), instead of the f.p.u. processor instruction FSINCOS -- the
latter has a well-known degradation of performance for large arguments,
x, while the C math library functions modulo the argument into a range
with full accuracy.
The original 1988 paper is available at
https://aip.scitation.org/doi/pdf/10.1063/1.168296 >
On Monday, July 4, 2022 at 2:13:59 PM UTC+2, Krishna Myneni wrote:...
Below is a comparison of using the library version of FSIN vs the fpu
instruction for FSIN in kForth:
17 set-precision
ok
\ x = 10,000.
1.0e4 fsin fs.
-3.0561438888825215e-01 ok \ with range reduction
1.0e4 fsincos fdrop fs.
-3.0561438888825215e-01 ok \ native fpu instruction
\ x = 100,000.
1.0e5 fsin fs.
3.5748797972016508e-02 ok \ with range reduction
1.0e5 fsincos fdrop fs.
3.5748797972016383e-02 ok \ native fpu instruction
\ x = 1,000,000.
1.0e6 fsin fs.
-3.4999350217129294e-01 ok \ with range reduction
1.0e6 fsincos fdrop fs.
-3.4999350217129177e-01 ok \ native fpu instruction
The above results may be compared with high precision reference values
for sin(x), e.g. using maxima, Wolfram Alpha, etc.
: TEST ( -- )
CR 1e4 fdup +e. fsin -0.3056143888882521413609100352325069742318500438618062391101551450e f- space +e.
CR 1e5 fdup +e. fsin 0.0357487979720165093164705006958088290090456925781088968546167365e f- space +e.
CR 1e6 fdup +e. fsin -0.3499935021712929521176524867807714690614066053287162738570590540e f- space +e. ;
FORTH> TEST
1.0000000000000000000e+0004 -1.2251484549086200104e-0017
1.0000000000000000000e+0005 -1.2865752842435018710e-0016
1.0000000000000000000e+0006 1.2060122865642508572e-0015 ok
For up to 1e5 we are still ok with double precision. Above that, accuracy decreases linearly.
It might be more important that the approximation is continuous and smooth than that it is highly accurate.
Krishna Myneni schrieb am Samstag, 2. Juli 2022 um 22:37:38 UTC+2:
Forth code for calculation over a large set of spherical Bessel and
spherical Neumann functions, with standard normalization, is available
in the kForth repos. -- see, for example,
https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/sph_bes_neu.4th
While the FSL provides a spherical Bessel function calculator for the
first ten functions, j_0(x) -- j_9(x), the present code extends the
calculation range to beyond j_100(x). The spherical Neumann functions,
n_l(x) (also often denoted by y_l(x)) are also computed over the same
range of l.
Test code is included, with reference values obtained from Wolfram
Alpha. The original source for this code is a paper in Computers in
Physics vol. 2, pp. 62--72 (1988), linked below, in which they provide
an accuracy analysis. The test code appears to show this version,
running under kForth, exceeds the accuracy reported in the paper for
double precision calculations. One should note that kForth's FSIN and
FCOS make use of the C math library routines for computing sin(x) and
cos(x), instead of the f.p.u. processor instruction FSINCOS -- the
latter has a well-known degradation of performance for large arguments,
x, while the C math library functions modulo the argument into a range
with full accuracy.
The original 1988 paper is available at
https://aip.scitation.org/doi/pdf/10.1063/1.168296 >
Just being curious: I stumbled over spherical Bessel functions only once
when modelling heat propagation for thick-walled steam generator
components. FWIW there were always physically measured initial
conditions along with their naturally limited low measuring accuracies,
like pressures and temperatures at component walls. Measurement
error would always outweigh theoretical calculation precisions.
So what practical applications do really require those discussed enormous calculation accuracies for Bessel functions??
On 7/4/22 12:36, Marcel Hendrix wrote:
On Monday, July 4, 2022 at 2:13:59 PM UTC+2, Krishna Myneni wrote:...
Below is a comparison of using the library version of FSIN vs the fpu
instruction for FSIN in kForth:
17 set-precision
ok
\ x = 10,000.
1.0e4 fsin fs.
-3.0561438888825215e-01 ok \ with range reduction
1.0e4 fsincos fdrop fs.
-3.0561438888825215e-01 ok \ native fpu instruction
\ x = 100,000.
1.0e5 fsin fs.
3.5748797972016508e-02 ok \ with range reduction
1.0e5 fsincos fdrop fs.
3.5748797972016383e-02 ok \ native fpu instruction
\ x = 1,000,000.
1.0e6 fsin fs.
-3.4999350217129294e-01 ok \ with range reduction
1.0e6 fsincos fdrop fs.
-3.4999350217129177e-01 ok \ native fpu instruction
The above results may be compared with high precision reference values
for sin(x), e.g. using maxima, Wolfram Alpha, etc.
: TEST ( -- )
CR 1e4 fdup +e. fsin
-0.3056143888882521413609100352325069742318500438618062391101551450e
f- space +e.
CR 1e5 fdup +e. fsin
0.0357487979720165093164705006958088290090456925781088968546167365e f-
space +e.
CR 1e6 fdup +e. fsin
-0.3499935021712929521176524867807714690614066053287162738570590540e
f- space +e. ;
FORTH> TEST
1.0000000000000000000e+0004 -1.2251484549086200104e-0017
1.0000000000000000000e+0005 -1.2865752842435018710e-0016
1.0000000000000000000e+0006 1.2060122865642508572e-0015 ok
For up to 1e5 we are still ok with double precision. Above that, accuracy
decreases linearly.
See the following link. The graph showing "abs(error in ULPs)" vs
"Biased Exponent of Extended Precision input Argument" indicates that
the log-log plot of error vs magnitude is linear, not the error vs
magnitude.
...
| Sysop: | Keyop |
|---|---|
| Location: | Huddersfield, West Yorkshire, UK |
| Users: | 300 |
| Nodes: | 16 (3 / 13) |
| Uptime: | 50:21:56 |
| Calls: | 6,712 |
| Calls today: | 5 |
| Files: | 12,243 |
| Messages: | 5,354,922 |
| Posted today: | 1 |