From Newsgroup: comp.lang.forth

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 ===


--- Synchronet 3.22a-Linux NewsLink 1.2

From Newsgroup: comp.lang.forth

On 8/15/26 07:32, Krishna Myneni wrote:

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.

...

=== 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
\-a-a-a-a-a native x87 trig functions are accurate.
\
\ 2. Range reduction is valid over the range,
\-a-a-a |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-a-a DP_RR_CUTOFF
9.2e18 fconstant-a-a DP_RR_MAX_ANGLE

\ BIG number constants and variables
create-a 5^52-a-a-a BIG_DBL CELLS allot-a-a 5 52 big_s^n-a 5^52 big-move create 10^52-a-a-a BIG_DBL CELLS allot-a 10 52 big_s^n 10^52 big-move
create BIG_2PI-a BIG_DBL CELLS allot
big 62831853071795864769252867665590057683943387987502116
BIG_2PI big-move

: exponent ( F: r -- ) ( -- n ) fexponent 1023 - ;

variable m
: freduce-range ( F: r1 -- r2)
-a-a-a \ Store sign, perform error checks
-a-a-a fdup f0< >r-a ( R: bsign)
-a-a-a fabs
-a-a-a fdup DP_RR_CUTOFF f< IF
-a-a-a-a-a r> IF fnegate THEN EXIT-a \ no rr correction
-a-a-a THEN
-a-a-a fdup DP_RR_MAX_ANGLE f> IF-a \ exceed limit of algorithm
-a-a-a-a-a cr ." Angle out of range." cr
-a-a-a-a-a -46 throw
-a-a-a THEN
-a-a-a fdup PI 2e f* f/ ftrunc>s m ! \ multiple of 2pi

-a-a-a \ Perform range reduction
-a-a-a fdup ffraction drop >r big-here dup 5^52 big>here
-a-a-a r> big*s 10^52 big+ >r-a ( R: abig bsign) ( F: r)
-a-a-a exponent 1 swap lshift
-a-a-a r@ swap big*s
-a-a-a r> r> IF bignegate THEN
-a-a-a big-here dup BIG_2PI big>here m @ big*s
-a-a-a big-
-a-a-a dup big0< swap dup bigabs <big# big#s swap bigsign #big>
-a-a-a >float IF
-a-a-a-a-a-a 1.0e-52 f*
-a-a-a ELSE
-a-a-a-a-a ." Conversion to double-precision float error!"
-a-a-a-a-a 0e
-a-a-a THEN
;
=== end rr-big.4th ===

I forgot to mention that the code above assumes a 64-bit Forth system
with a separate fp stack.

--
KM
--- Synchronet 3.22a-Linux NewsLink 1.2