Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #135379

Re: Range Reduction Using Big Number Arithmetic

From Krishna Myneni <krishna.myneni@ccreweb.org>
Newsgroups comp.lang.forth
Subject Re: Range Reduction Using Big Number Arithmetic
Date 2026-08-20 22:01 -0500
Organization A noiseless patient Spider
Message-ID <1168f21$bnk$1@dont-email.me> (permalink)
References <115pm90$3cl9h$1@dont-email.me>

Show all headers | View raw

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.
...
> 
> 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.
> 
...
> 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 )
> 
> ...
> 
>    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

I also tested the cosine function with FREDUCED-RANGE:

FCOS( x )
FCOS( x' )
gcc cos( x )

over the same set of large angle inputs. The gcc sin( x ) and cos( x ) 
functions are doing what they are supposed to do, which is to provide 
the nearest rounded number to the ideal accurate result which can be 
represented in the IEEE-754 binary 64, or double-precision format.

So why do the gcc sin( x ) and cos( x ) values shown in the tests (see 
below) differ from the numbers shown for Wolfram alpha sin( x ) and cos( 
x )? The Wolfram numbers are computed at much higher precision and the 
outputs have 60+ decimal digits. In the test table, they are rounded to 
provide 20 significant decimal digits for comparison with the above 
tests. Double precision only provides a 53 bit significand, which limits 
the number of significant decimal digits to 15 or 16, if the computation 
was accurate to within 1 ulp. The comparison between Wolfram alpha 
numbers and the gcc trig function numbers are good to 15 or 16 digits.

Our FREDUCE-RANGE implementation using big number arithmetic to reduce 
the large angle, x, to a smaller equivalent angle, x', for sin() and 
cos() computation using the x87 fpu functions, do not agree as well as 
gcc trig functions agree with the Wolfram alpha numbers.

FREDUCE-RANGE does provide a reduced angle, x', which gives far more 
accurate results than using the large angle, x, with the native fpu 
instructions, many orders of magnitude better. Sometimes the results of 
using FREDUCE-RANGE with native FSIN and FCOS are the same as those 
given by gcc; however, improvements in the implementation should be able 
to produce the same results as gcc, which are as accurate as possible 
within double-precision format.

--
Krishna


=== Updated test results ===
15 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

------------------------------------------------


15 August 2024

Tests for FCOS

   x = 5e05
   -9.84061006120338|31077e-01  FCOS( x )
   -9.84061006120338|31077e-01  FCOS( x' )
   -9.84061006120338|19974e-01  gcc cos( x )
   -9.84061006120338|24936e-01  Wolfram alpha

   x = 1e06
    9.3675212753314|518466e-01  FCOS( x )
    9.3675212753314|485159e-01  FCOS( x' )
    9.3675212753314|474057e-01  gcc cos( x )
    9.3675212753314|478694e-01  Wolfram alpha

   x = 1e07
   -9.072703861817|4493667e-01  FCOS( x )
   -9.072703861817|3960760e-01  FCOS( x' )
   -9.072703861817|3960760e-01  gcc cos( x )
   -9.072703861817|3956116e-01  Wolfram alpha

   x = 1e08
   -3.63385089355|81048664e-01  FCOS( x )
   -3.63385089355|69074909e-01  FCOS( x' )
   -3.63385089355|69052705e-01  gcc cos( x )
   -3.63385089355|69055387e-01  Wolfram alpha

   x = 1e09
    8.37887181363|19961793e-01  FCOS( x )
    8.37887181363|90227809e-01  FCOS( x' )
    8.37887181363|90238911e-01  gcc cos( x )
    8.37887181363|90233439e-01  Wolfram alpha

   x = 1e10
    8.731196226|8313225710e-01  FCOS( x )
    8.731196226|7685605532e-01  FCOS( x' )
    8.731196226|7685605532e-01  gcc cos( x )
    8.731196226|7685600118e-01  Wolfram alpha

   x = 1e11
    3.70847792|04514978224e-01  FCOS( x )
    3.70847792|16471114305e-01  FCOS( x' )
    3.70847792|16471114305e-01  gcc cos( x )
    3.70847792|16471115901e-01  Wolfram alpha

   x = 1e12
    7.9144630|263980797480e-01  FCOS( x )
    7.9144630|185289066571e-01  FCOS( x' )
    7.9144630|185289022162e-01  gcc cos( x )
    7.9144630|185289027005e-01  Wolfram alpha

   x = 1e13
    9.573637|2061999153171e-01  FCOS( x )
    9.573637|1690083987840e-01  FCOS( x' )
    9.573637|1690083998942e-01  gcc cos( x )
    9.573637|1690083993528e-01  Wolfram alpha

   x = 1e14
   -9.778282|6100899145061e-01  FCOS( x )
   -9.778282|8796853249465e-01  FCOS( x' )
   -9.778282|8796853249465e-01  gcc cos( x )
   -9.778282|8796853247834e-01  Wolfram alpha

   x = 1e15
   -5.1319|484273953741571e-01  FCOS( x )
   -5.1319|373778697019439e-01  FCOS( x' )
   -5.1319|373778697030541e-01  gcc cos( x )
   -5.1319|373778697025223e-01  Wolfram alpha

   x = 1e16
   -6.261|7823589904142434e-01  FCOS( x )
   -6.261|6819813308610687e-01  FCOS( x' )
   -6.261|6819813308621789e-01  gcc cos( x )
   -6.261|6819813308617176e-01  Wolfram alpha

   x = 1e17
   -8.85|49751667144138700e-01  FCOS( x )
   -8.85|55732829763078584e-01  FCOS( x' )
   -8.85|55732829763067482e-01  gcc cos( x )
   -8.85|55732829763068505e-01  Wolfram alpha

   x = 1e18
    1.1|965025504785124777e-01  FCOS( x )
    1.1|837199021870972726e-01  FCOS( x' )
    1.1|837199021871072646e-01  gcc cos( x )
    1.1|837199021871073261e-01  Wolfram alpha

-----------------------------------------------
=== end of updated test results ===

Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

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