Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #135346
| From | Krishna Myneni <krishna.myneni@ccreweb.org> |
|---|---|
| Newsgroups | comp.lang.forth |
| Subject | Range Reduction Using Big Number Arithmetic |
| Date | 2026-08-15 07:32 -0500 |
| Organization | A noiseless patient Spider |
| Message-ID | <115pm90$3cl9h$1@dont-email.me> (permalink) |
A simple Forth implementation for reducing large angles (|x| > 500000
rad) for accurate evaluation by the x87 native fpu trig instructions
FSIN FCOS FSINCOS FTAN is given below. The Forth code relies on the
Forth Scientific Library (FSL #47) big number arithmetic module.
The reduction of a large angle by subtracting an integer multiple of 2pi,
|x'| = |x| - m*2pi
to bring the angle x' into the range minus pi to plus pi, where the x87
fpu instructions have highest accuracy in double precision floating
point evaluation of the trig functions. Peforming this reduction with
double-precision floating point arithmetic will actually degrade the
accuracy of evaluating the trig functions! The reduction requires
greater arithmetic precision.
The implementation with big number arithmetic is
FREDUCE-RANGE ( F: x -- x' )
With the present implementation, the argument x is limited to a range,
|x| < 9.2e18 ( radians )
Test results of the effect of range reduction implementation on the
accuracy of the native FSIN instruction are provided. Shown is the rapid
loss of accuracy for |x| > 10^6 rad when providing x as the argument to
the native FSIN fpu instruction (only 3 significant digits at x = 1e18
rad). In contrast, the reduced-range angle, x', computed by
FREDUCE-RANGE maintains 14 to 15 significant digits in the result of
FSIN over the range of x stated above. A vertical separator "|" shows
the decimal place where the native FSIN result deviates from the
range-reduced result. Note how quickly the separator moves to the left
as the angle argument is increased by factors of 10.
The test results have three sections. In the first section, the angle
defined by the constant DP_RR_CUTOFF is set to a low value (100 rad) in
order to determine a cutoff limit for using the native FSIN fpu
instruction without range-reduction. Below this cutoff, performing range
reduction does not help and may provide a slightly less accurate answer
for the present algorithm. The second section shows results for angle
arguments above the value of DP_RR_CUTOFF actually used in the code. The
third section shows what happens when the max angle limit, given by the
fconstant DP_RR_MAX_ANGLE, is exceeded.
If you are using a C-based Forth, the libc function may or may not
provide range-reduction. Comparison with the GNU library function
sin( x ),
which does perform range-reduction is shown in the test results.
--
Krishna Myneni
=== begin test results ===
5 August 2026
\
\ 53 digits of 2*pi
\ 6 283 185 307 179 586 476 925 286 766 559 005 768 394 338 798 750 211 6
\
Test performance of big number implementation of FREDUCE-RANGE
x = large angle in radians
x' is the reduced angle returned by "<x> FREDUCE-RANGE"
FSIN = native x87 instruction
kForth-64/32 specific mmands:
: S1 ( F: x -- sin[x]) fsincos fdrop ; \ FSIN( x )
: S2 ( F: x -- sin[x]) freduce-range S1 ; \ reduced x FSIN
: S3 ( F: x -- sin[x]) fsin ; \ use gcc sin( x )
----------------------------------------------
DP_RR_CUTOFF = 1e02
x = 1e03
8.26879540532002|52158e-01 FSIN( x )
8.26879540532002|63260e-01 FSIN( x' )
8.26879540532002|52158e-01 gcc sin( x )
8.26879540532002|56026e-01 Wolfram alpha
x = 1e04
-3.05614388888252|15310e-01 FSIN( x )
-3.05614388888252|31963e-01 FSIN( x' )
-3.05614388888252|15310e-01 gcc sin( x )
-3.05614388888252|14136e-01 Wolfram alpha
x = 5e04
-9.99840189089789|55498e-01 FSIN( x )
-9.99840189089789|55498e-01 FSIN( x' )
-9.99840189089789|55498e-01 gcc sin( x )
-9.99840189089789|60183e-01 Wolfram alpha
x = 6e04
9.57466750100169|57241e-01 FSIN( x )
9.57466750100169|68344e-01 FSIN( x' )
9.57466750100169|68344e-01 gcc sin( x )
9.57466750100169|63469e-01 Wolfram alpha
x = 1e05
3.5748797972016|382873e-02 FSIN( x )
3.5748797972016|674307e-02 FSIN( x' )
3.5748797972016|507773e-02 gcc sin( x )
3.5748797972016|509316e-02 Wolfram alpha
x = 2e05
-7.14518952125|19642153e-02 FSIN( x )
-7.14518952125|20225020e-02 FSIN( x' )
-7.14518952125|19905831e-02 gcc sin( x )
-7.14518952125|19901159e-02 Wolfram alpha
-------------------------------------------------
x = 5e05
1.77831201518258|26394e-01 FSIN( x )
1.77831201518258|79129e-01 FSIN( x' )
1.77831201518258|90232e-01 gcc sin( x )
1.77831201518258|90009e-01 Wolfram Alpha
x = 1e06
-3.4999350217129|177043e-01 FSIN( x )
-3.4999350217129|271412e-01 FSIN( x' )
-3.4999350217129|293616e-01 gcc sin( x )
-3.4999350217129|295212e-01 Wolfram alpha
x = 1e07
4.205477931907|7079114e-01 FSIN( x )
4.205477931907|8250399e-01 FSIN( x' )
4.205477931907|8250399e-01 gcc sin( x )
4.205477931907|8249130e-01 Wolfram alpha
x = 1e08
9.31639027109|67927412e-01 FSIN( x )
9.31639027109|72590349e-01 FSIN( x' )
9.31639027109|72601451e-01 gcc sin( x )
9.31639027109|72600803e-01 Wolfram alpha
x = 1e09
5.4584344944|977825076e-01 FSIN( x )
5.4584344944|869955807e-01 FSIN( x' )
5.4584344944|869955807e-01 gcc sin( x )
5.4584344944|869956424e-01 Wolfram alpha
x = 1e10
-4.875060250|7626999982e-01 FSIN( x )
-4.875060250|8751067488e-01 FSIN( x' )
-4.875060250|8751067488e-01 gcc sin( x )
-4.875060250|8751069153e-01 Wolfram alpha
x = 1e11
9.28693660|54433544200e-01 FSIN( x )
9.28693660|49659196616e-01 FSIN( x' )
9.28693660|49659196616e-01 gcc sin( x )
9.28693660|49659195286e-01 Wolfram alpha
x = 1e12
-6.1123870|135796987135e-01 FSIN( x )
-6.1123870|237688904261e-01 FSIN( x' )
-6.1123870|237688948670e-01 gcc sin( x )
-6.1123870|237688949819e-01 Wolfram alpha
x = 1e13
-2.888852|8249228406786e-01 FSIN( x )
-2.888852|9481752533989e-01 FSIN( x' )
-2.888852|9481752511785e-01 gcc sin( x )
-2.888852|9481752512226e-01 Wolfram alpha
x = 1e14
-2.09408|43338349970915e-01 FSIN( x )
-2.09408|30749645225617e-01 FSIN( x' )
-2.09408|30749645231168e-01 gcc sin( x )
-2.09408|30749645230269e-01 Wolfram alpha
x = 1e15
8.58272|13247637335058e-01 FSIN( x )
8.58272|79317023585925e-01 FSIN( x' )
8.58272|79317023585925e-01 gcc sin( x )
8.58272|79317023583552e-01 Wolfram alpha
x = 1e16
7.796|7994516106697844e-01 FSIN( x )
7.796|8800660697878957e-01 FSIN( x' )
7.796|8800660697878957e-01 gcc sin( x )
7.796|8800660697875024e-01 Wolfram alpha
x = 1e17
-4.64|64410893576429951e-01 FSIN( x )
-4.64|53010483537254816e-01 FSIN( x' )
-4.64|53010483537271469e-01 gcc sin( x )
-4.64|53010483537269615e-01 Wolfram alpha
x = 1e18
-9.92|81610405300346756e-01 FSIN( x )
-9.92|96932074040522576e-01 FSIN( x' )
-9.92|96932074040511473e-01 gcc sin( x )
-9.92|96932074040507621e-01 Wolfram alpha
------------------------------------------------
DP_RR_MAX_ANGLE = 1e19
x = 1e19
-nan FSIN( x )
VM ERROR(-270): Division overflow FSIN( x' )
-9.2706316604865035558e-01 gcc sin( x )
-9.2706316604865038523e-01 Wolfram alpha
------------------------------------------------
=== end test results ===
=== begin rr-big.4th ===
\ rr-big.4th
\
\ Perform range reduction of large angle x (radians) for
\ accurate trigonometric functions with x87 using the FSL
\ big number arithmetic module.
\
\ K. Myneni, 2026-08-15
\
\ Notes:
\ 1. Angles |x| < DP_RR_CUTOFF are not reduced. The
\ native x87 trig functions are accurate.
\
\ 2. Range reduction is valid over the range,
\ |x| < DP_RR_MAX_ANGLE
\
include ans-words
include modules
include ieee-754
include fsl/big
include fsl/extras/big-extras
8 constant BIG_DBL
5.0e05 fconstant DP_RR_CUTOFF
9.2e18 fconstant DP_RR_MAX_ANGLE
\ BIG number constants and variables
create 5^52 BIG_DBL CELLS allot 5 52 big_s^n 5^52 big-move
create 10^52 BIG_DBL CELLS allot 10 52 big_s^n 10^52 big-move
create BIG_2PI BIG_DBL CELLS allot
big 62831853071795864769252867665590057683943387987502116
BIG_2PI big-move
: exponent ( F: r -- ) ( -- n ) fexponent 1023 - ;
variable m
: freduce-range ( F: r1 -- r2)
\ Store sign, perform error checks
fdup f0< >r ( R: bsign)
fabs
fdup DP_RR_CUTOFF f< IF
r> IF fnegate THEN EXIT \ no rr correction
THEN
fdup DP_RR_MAX_ANGLE f> IF \ exceed limit of algorithm
cr ." Angle out of range." cr
-46 throw
THEN
fdup PI 2e f* f/ ftrunc>s m ! \ multiple of 2pi
\ Perform range reduction
fdup ffraction drop >r big-here dup 5^52 big>here
r> big*s 10^52 big+ >r ( R: abig bsign) ( F: r)
exponent 1 swap lshift
r@ swap big*s
r> r> IF bignegate THEN
big-here dup BIG_2PI big>here m @ big*s
big-
dup big0< swap dup bigabs <big# big#s swap bigsign #big>
>float IF
1.0e-52 f*
ELSE
." Conversion to double-precision float error!"
0e
THEN
;
=== end rr-big.4th ===
Back to comp.lang.forth | Previous | Next — Next in thread | Find similar | Unroll thread
Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-15 07:32 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-15 07:48 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-20 22:01 -0500
Re: Range Reduction Using Big Number Arithmetic albert@spenarnc.xs4all.nl - 2026-08-22 19:31 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-31 08:49 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-28 10:39 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-28 10:42 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-03 18:03 -0500
Re: Range Reduction Using Big Number Arithmetic marcel hendrix <mhx@iae.nl> - 2026-09-04 07:29 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-04 07:45 -0500
Re: Range Reduction Using Big Number Arithmetic marcel hendrix <mhx@iae.nl> - 2026-09-04 15:29 +0200
Re: Range Reduction Using Big Number Arithmetic albert@spenarnc.xs4all.nl - 2026-09-05 12:27 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-05 08:39 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-09 07:37 -0500
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-09 18:49 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-09 15:06 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-10 19:03 -0500
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 14:56 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 09:17 -0500
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 17:04 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 10:46 -0500
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 13:24 -0500
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 23:21 +0200
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-12 00:42 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 18:33 -0500
Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-12 11:04 +0200
Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-12 10:33 -0500
csiph-web