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

Groups > comp.lang.forth > #135346 > unrolled thread

Range Reduction Using Big Number Arithmetic

Started byKrishna Myneni <krishna.myneni@ccreweb.org>
First post2026-08-15 07:32 -0500
Last post2026-09-12 10:33 -0500
Articles 20 on this page of 27 — 4 participants

Back to article view | Back to comp.lang.forth

Contents

  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

Page 1 of 2  [1] 2  Next page →

#135346 — Range Reduction Using Big Number Arithmetic

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-15 07:32 -0500
SubjectRange Reduction Using Big Number Arithmetic
Message-ID<115pm90$3cl9h$1@dont-email.me>
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 ===

[toc] | [next] | [standalone]

#135347

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-15 07:48 -0500
Message-ID<115pn6e$3cout$1@dont-email.me>
In reply to#135346
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
> \      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 ===
> 

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

--
KM

[toc] | [prev] | [next] | [standalone]

#135379

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-20 22:01 -0500
Message-ID<1168f21$bnk$1@dont-email.me>
In reply to#135346
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 ===

[toc] | [prev] | [next] | [standalone]

#135392

Fromalbert@spenarnc.xs4all.nl
Date2026-08-22 19:31 +0200
Message-ID<nnd$75752e09$145a90a2@8eb1fe0acd09f684>
In reply to#135379
In article <1168f21$bnk$1@dont-email.me>,
Krishna Myneni  <krishna.myneni@ccreweb.org> wrote:
>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.
>>
>>
<SNIP>
>   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
>------------------------------------------------

Without range reduction on a 8087.
1E18 FSIN 30 FS.
-8.751442827729692093E-1

S[ ] OK 1E18

S[ ] OK 20 FS.
D.E0B6B3A763FFFF6000_E
S[ ] OK D.E0B6B3A763FFFF61_E
S[ ] OK 20 FS.
D.E0B6B3A763FFFF6000_E  \ So the number is exact.
S[ ] OK D.E0B6B3A763FFFF5_E FSIN DECIMAL 20 FS.
-8.432138792941086135E-1
S[ ] OK HEX D.E0B6B3A763FFFF7_E FSIN DECIMAL 20 FS.
-9.036572665558625520E-1

The difference is 0.3E-1 for the neighboring numbers.
Assuming that is approximately correct,
we arrive at
   -0.992 +- 0.015
assuming the original number was found by correctly rounding.

If I convert it to the 64 bit IEEE I get
D.E0B6B3A763FF800000_E
then the result would be much worse.
HEX  D.E0B6B3A763FF7_E FSIN DECIMAL 20 FS.
-8.061383098048763381E-1  OK
HEX  D.E0B6B3A763FF8_E FSIN DECIMAL 20 FS.
6.016457844727310581E-1  OK
HEX  D.E0B6B3A763FF9_E FSIN DECIMAL 20 FS.
-3.462052687889450736E-1  OK

That would amount to 0.6 +- 0.7 .

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

Groetjes Albert
-- 
The Chinese government is satisfied with its military superiority over USA.
The next 5 year plan has as primary goal to advance life expectancy
over 80 years, like Western Europe.

[toc] | [prev] | [next] | [standalone]

#135496

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-31 08:49 -0500
Message-ID<11740q0$vgtv$1@dont-email.me>
In reply to#135392
On 8/22/26 12:31, albert@spenarnc.xs4all.nl wrote:
> In article <1168f21$bnk$1@dont-email.me>,
> Krishna Myneni  <krishna.myneni@ccreweb.org> wrote:
>> 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.
>>>
>>>
> <SNIP>
>>    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
>> ------------------------------------------------
> 
> Without range reduction on a 8087.
> 1E18 FSIN 30 FS.
> -8.751442827729692093E-1
> 
> S[ ] OK 1E18
> 
> S[ ] OK 20 FS.
> D.E0B6B3A763FFFF6000_E
> S[ ] OK D.E0B6B3A763FFFF61_E
> S[ ] OK 20 FS.
> D.E0B6B3A763FFFF6000_E  \ So the number is exact.
> S[ ] OK D.E0B6B3A763FFFF5_E FSIN DECIMAL 20 FS.
> -8.432138792941086135E-1
> S[ ] OK HEX D.E0B6B3A763FFFF7_E FSIN DECIMAL 20 FS.
> -9.036572665558625520E-1
> 
> The difference is 0.3E-1 for the neighboring numbers.
> Assuming that is approximately correct,
> we arrive at
>     -0.992 +- 0.015
> assuming the original number was found by correctly rounding.
> 
> If I convert it to the 64 bit IEEE I get
> D.E0B6B3A763FF800000_E
> then the result would be much worse.
> HEX  D.E0B6B3A763FF7_E FSIN DECIMAL 20 FS.
> -8.061383098048763381E-1  OK
> HEX  D.E0B6B3A763FF8_E FSIN DECIMAL 20 FS.
> 6.016457844727310581E-1  OK
> HEX  D.E0B6B3A763FF9_E FSIN DECIMAL 20 FS.
> -3.462052687889450736E-1  OK
> 
> That would amount to 0.6 +- 0.7 .
> ...

The error due to finite precision of the argument is important. But, it 
is a separate issue from what we are discussing here. Here, the goal is 
to quantify the error due to not having range reduction, for an argument 
which can be *exactly represented* in double precision.

The argument, 1e18, is exactly representable.

20 set-precision
fvariable r  \ assumes fvariable creates a double-precision fp variable

1e18 r f!
r f@ fs.
1.0000000000000000000e+18  ok

55 set-precision
r f@ fs.
1.000000000000000000000000000000000000000000000000000000e+18  ok


To look at the interval between adjacent fp numbers using 64-bit Forth,

hex
r @ u.
43ABC16D674EC800  ok

\ next value above 1e18 in double precision fp
43ABC16D674EC801 r !  \ significand with incremented fraction
decimal
20 set-precision

r f@ fs.
1.0000000000000001280e+18  ok

r f@ 1e18 f- fs.
1.2800000000000000000e+02  ok

Thus, the next number above 1e18 in double precision fp is r = 1e18 + 128.

The FSIN of 1e18 and 1e18+128 will be different. Using FSIN which calls 
gcc sin(x) function,

r f@ fsin 1e18 fsin f- fs.
1.7663442831911686515e+00  ok

Using FSIN which uses x87 FSIN instruction gives
r f@ fsin-x87 1e18 fsin-x87 f- fs.
1.7670065804470365123e+00  ok

The above number, computed with the native FSIN-x87, is comparable but 
not the same as using the range-reduced evaluation of FSIN e.g. gcc 
sin(x). What is the error? First, obtain the reference value for 
sin(1e18) from a high-precision arithmetic source.

Reference value (from Wolfram alpha)
sin(1e18) =
-0.992969320740405076209553017263630270859984568678211649303457453...

1e18 fsin fs.  \ using gcc sin(x)
-9.9296932074040511473e-01  ok

1e18 fsin-x87 fs.  \ using x87 FSIN
-9.9281610405300346756e-01  ok

We see that gcc sin(x) provides 16 significant digits, while the native 
FSIN instruction provides four significant digits.

This is also why having a high quality decimal string to binary floating 
point conversion is important! If the conversion is off by +/-1 ulp of 
the binary representation, then large arguments, which are exactly 
represented, to a properly range-reduced forth FSIN will have huge loss 
of accuracy. The loss in accuracy for FSIN is easy to observe even with 
exact arguments |x| ~ 10^6.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135460

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-28 10:39 -0500
Message-ID<116sa3b$2dpj8$4@dont-email.me>
In reply to#135379
On 8/20/26 22:01, Krishna Myneni wrote:
> 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.
>>
...
One can examine the glibc code for double precision sin(x) at the link 
below.

https://github.com/bminor/glibc/blob/6059938728a98270b9706488887f43baa0471eba/sysdeps/ieee754/dbl-64/s_sin.c

=== clip from above file ===
...
/*******************************************************************/
/* An ultimate sin routine. Given an IEEE double machine number x  */
/* it computes the rounded value of sin(x).			   */
/*******************************************************************/
#ifndef IN_SINCOS
double
SECTION
__sin (double x)
{
...
=== end of clip ===

There are several ranges of the input parameter x, each of which are 
handled differently:

I. 2^-26 < |x| < 0.855469;  no range reduction needed

II. 0.855469  < |x| < 2.426265; ?

III. 2.426265 < |x| < 105414350; uses function reduce_sincos()

IV.  105414350 < |x| < 2^1024; uses function _branred()

V. 105414350 < |x| < 2^1024; set error condition


The only ranges where it looks like it is performing range reduction on 
the argument are III and IV, above which the double-precision number is 
out of range. In II, the argument is being shifted by a constant, 
probably pi/4, and its cosine is being evaluated and the sign is adjusted.

It will be interesting to see the details of branred().

--
KM

[toc] | [prev] | [next] | [standalone]

#135461

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-08-28 10:42 -0500
Message-ID<116sa8m$2dpj8$5@dont-email.me>
In reply to#135460
On 8/28/26 10:39, Krishna Myneni wrote:
> On 8/20/26 22:01, Krishna Myneni wrote:
>> 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.
>>>
> ...
> One can examine the glibc code for double precision sin(x) at the link 
> below.
> 
> https://github.com/bminor/glibc/ 
> blob/6059938728a98270b9706488887f43baa0471eba/sysdeps/ieee754/dbl-64/ 
> s_sin.c
> 
> === clip from above file ===
> ...
> /*******************************************************************/
> /* An ultimate sin routine. Given an IEEE double machine number x  */
> /* it computes the rounded value of sin(x).               */
> /*******************************************************************/
> #ifndef IN_SINCOS
> double
> SECTION
> __sin (double x)
> {
> ...
> === end of clip ===
> 
> There are several ranges of the input parameter x, each of which are 
> handled differently:
> 
...
> V. 105414350 < |x| < 2^1024; set error condition
> 
> 

Copy and paste error:

V. |x| > 2^1024;  set error condition

[toc] | [prev] | [next] | [standalone]

#135545

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-03 18:03 -0500
Message-ID<117cuca$4f48$1@dont-email.me>
In reply to#135460
On 8/28/26 10:39, Krishna Myneni wrote:
> On 8/20/26 22:01, Krishna Myneni wrote:
>> 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.
>>>
> ...
> One can examine the glibc code for double precision sin(x) at the link 
> below.
> 
> https://github.com/bminor/glibc/ 
> blob/6059938728a98270b9706488887f43baa0471eba/sysdeps/ieee754/dbl-64/ 
> s_sin.c
> 
> === clip from above file ===
> ...
> /*******************************************************************/
> /* An ultimate sin routine. Given an IEEE double machine number x  */
> /* it computes the rounded value of sin(x).               */
> /*******************************************************************/
> #ifndef IN_SINCOS
> double
> SECTION
> __sin (double x)
> {
> ...
> === end of clip ===
> 
> There are several ranges of the input parameter x, each of which are 
> handled differently:
> 
> I. 2^-26 < |x| < 0.855469;  no range reduction needed
> 
> II. 0.855469  < |x| < 2.426265; ?
> 
> III. 2.426265 < |x| < 105414350; uses function reduce_sincos()
> 
> IV.  105414350 < |x| < 2^1024; uses function _branred()
> 
> V. 105414350 < |x| < 2^1024; set error condition
> 
> 
> The only ranges where it looks like it is performing range reduction on 
> the argument are III and IV, above which the double-precision number is 
> out of range. In II, the argument is being shifted by a constant, 
> probably pi/4, and its cosine is being evaluated and the sign is adjusted.
> 
> It will be interesting to see the details of branred().
> 

I've studied the branred.c code from glibc (version 2.4.2.9000, from 
Github repo bminor/glibc), and I've made a stand-alone set of files 
which compile to provide accurate, range-reduced sin(x) and cos(x) 
functions -- yes, it works at x=1e100 rad!

The interesting part is that branred() uses double-double floating point 
arithmetic (which provides a 106-bit significand, over an exponent range 
the same as ordinary double-precision). Consequently, the glibc code for 
sin(x) and cos(x) functions set the fpu precision and rounding mode to 
double-precision and round to nearest, in order for double-double 
precision arithmetic to work.

kForth's startup code always initializes the fpu to double-precision, 
round to nearest mode, but apparently this is not the default for glibc.

Example: the initialization code in kForth-32 vm32-common.s is,

# set kForth's default fpu settings
L_initfpu:
         LDSP
         fnstcw NDPcw           # save the NDP control word
         movl NDPcw, %ecx
         andb $240, %ch         # mask the high byte
         orb  $2,  %ch          # set double precision, round near
         mov  %ecx, (%ebx)
         fldcw (%ebx)
         ret

Note the '$' prefix for numbers in the GNU assembler indicates base 10, 
decimal not hex! NDPcw is the local 16-bit storage for saving the x87 
control word, in case it needs to be restored.

I am porting the range reduction code in branred.c to Forth.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135549

Frommarcel hendrix <mhx@iae.nl>
Date2026-09-04 07:29 +0200
Message-ID<117dl09$u1dl$1@dont-email.me>
In reply to#135545
On 9/4/2026 1:03 AM, Krishna Myneni wrote:
[..]
> I am porting the range reduction code in branred.c to Forth.
Please note that setting the FPU rounding mode may affect calls to/by 
foreign code (e.g. callbacks).

-marcel

[toc] | [prev] | [next] | [standalone]

#135553

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-04 07:45 -0500
Message-ID<117eei7$kpil$1@dont-email.me>
In reply to#135549
On 9/4/26 00:29, marcel hendrix wrote:
> On 9/4/2026 1:03 AM, Krishna Myneni wrote:
> [..]
>> I am porting the range reduction code in branred.c to Forth.
> Please note that setting the FPU rounding mode may affect calls to/by 
> foreign code (e.g. callbacks).
> 
> -marcel

That would be a problem for signal handlers and interrupt service 
routines, but only if they deal with floating point arithmetic. That's 
probably pretty rare.

In general, if double-precision and round-to-nearest are not your 
default mode, words which change the fpu control word should save the 
prior control word and restore it before they exit.

JVN's double-double arithmetic Forth code assumes the fpu control 
settings and does not set or restore the fpu control word.

It is also important that the Forth interpreter's conversion of decimal 
input to double precision floating point must be accurate to within 1 
ulp for double-double arithmetic to work. These issues we have been 
talking about, decimal string to floating point binary and vice versa 
and range reduction, are connected to the accuracy of floating point 
calculations.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135554

Frommarcel hendrix <mhx@iae.nl>
Date2026-09-04 15:29 +0200
Message-ID<117eh3s$27m81$1@dont-email.me>
In reply to#135553
On 9/4/2026 2:45 PM, Krishna Myneni wrote:
> On 9/4/26 00:29, marcel hendrix wrote:
>> On 9/4/2026 1:03 AM, Krishna Myneni wrote:
>> [..]
>>> I am porting the range reduction code in branred.c to Forth.
>> Please note that setting the FPU rounding mode may affect calls to/by 
>> foreign code (e.g. callbacks).
>>
>> -marcel
> 
> That would be a problem for signal handlers and interrupt service 
> routines, but only if they deal with floating point arithmetic. That's 
> probably pretty rare.
I have ported all NRC functions (iForth's gaussj vsn 1.79 library). Many 
of those functions (e.g. the root and ODE solvers) require function 
pointers. In iForth use that means C calls Forth vv. through callbacks, 
i.e., the code for the C callback is written in iForth. It is not 
uncommon that a callback is nested, e.g., when a Forth word X uses a few 
NRC calls and then X's address is passed as an argument to one of the 
NRC ODE solvers which is written in C.

The problem with DD+E.R I showed here recently seems to be related to 
the words 53-bits! and 80-bits! setting the mantissa size but assuming 
the rounding mode is even. For some (bitrot) reason I changed the global 
rounding mode at some point without realizing (or knowing) that it would 
affect the double-double package. ISTR that setting the FPU rounding 
mode in assembler was quite a hassle and decided that a shortcut would 
be a good idea.

-marcel

[toc] | [prev] | [next] | [standalone]

#135565

Fromalbert@spenarnc.xs4all.nl
Date2026-09-05 12:27 +0200
Message-ID<nnd$208d7f49$2131b7cd@4402cfd8cda4c5d6>
In reply to#135554
In article <117eh3s$27m81$1@dont-email.me>, marcel hendrix  <mhx@iae.nl> wrote:
>On 9/4/2026 2:45 PM, Krishna Myneni wrote:
>> On 9/4/26 00:29, marcel hendrix wrote:
>>> On 9/4/2026 1:03 AM, Krishna Myneni wrote:
>>> [..]
>>>> I am porting the range reduction code in branred.c to Forth.
>>> Please note that setting the FPU rounding mode may affect calls to/by
>>> foreign code (e.g. callbacks).
>>>
>>> -marcel
>>
>> That would be a problem for signal handlers and interrupt service
>> routines, but only if they deal with floating point arithmetic. That's
>> probably pretty rare.
>I have ported all NRC functions (iForth's gaussj vsn 1.79 library). Many
>of those functions (e.g. the root and ODE solvers) require function
>pointers. In iForth use that means C calls Forth vv. through callbacks,
>i.e., the code for the C callback is written in iForth. It is not
>uncommon that a callback is nested, e.g., when a Forth word X uses a few
>NRC calls and then X's address is passed as an argument to one of the
>NRC ODE solvers which is written in C.
>
>The problem with DD+E.R I showed here recently seems to be related to
>the words 53-bits! and 80-bits! setting the mantissa size but assuming
>the rounding mode is even. For some (bitrot) reason I changed the global
>rounding mode at some point without realizing (or knowing) that it would
>affect the double-double package. ISTR that setting the FPU rounding
>mode in assembler was quite a hassle and decided that a shortcut would
>be a good idea.
>
>-marcel

This set me thinking ...

: FROUND   nearest-mode rounding: (FROUND) ;
: FLOOR   down-mode rounding: (FROUND) ;

Apparently I set truncate-mode  as default rounding.

\ Set  rounding-mode  for the duration of a definition.
: rounding:   set-rounding-mode  CO  truncate-mode set-rounding-mode ;
CODE (FROUND)   FRNDINT, NEXT, END-CODE

Why did I do that? Is nearest-mode not more appropriate?

default-mode 0 10 LSHIFT + CONSTANT nearest-mode
default-mode 1 10 LSHIFT + CONSTANT down-mode
default-mode 2 10 LSHIFT + CONSTANT up-mode
default-mode 3 10 LSHIFT + CONSTANT truncate-mode

Groetjes Albert
-- 
The Chinese government is satisfied with its military superiority over USA.
The next 5 year plan has as primary goal to advance life expectancy
over 80 years, like Western Europe.

[toc] | [prev] | [next] | [standalone]

#135568

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-05 08:39 -0500
Message-ID<117h62o$1hpi9$1@dont-email.me>
In reply to#135565
On 9/5/26 05:27, albert@spenarnc.xs4all.nl wrote:
> In article <117eh3s$27m81$1@dont-email.me>, marcel hendrix  <mhx@iae.nl> wrote:
>> On 9/4/2026 2:45 PM, Krishna Myneni wrote:
>>> On 9/4/26 00:29, marcel hendrix wrote:
>>>> On 9/4/2026 1:03 AM, Krishna Myneni wrote:
>>>> [..]
>>>>> I am porting the range reduction code in branred.c to Forth.
>>>> Please note that setting the FPU rounding mode may affect calls to/by
>>>> foreign code (e.g. callbacks).
>>>>
>>>> -marcel
>>>
>>> That would be a problem for signal handlers and interrupt service
>>> routines, but only if they deal with floating point arithmetic. That's
>>> probably pretty rare.
>> I have ported all NRC functions (iForth's gaussj vsn 1.79 library). Many
>> of those functions (e.g. the root and ODE solvers) require function
>> pointers. In iForth use that means C calls Forth vv. through callbacks,
>> i.e., the code for the C callback is written in iForth. It is not
>> uncommon that a callback is nested, e.g., when a Forth word X uses a few
>> NRC calls and then X's address is passed as an argument to one of the
>> NRC ODE solvers which is written in C.
>>
>> The problem with DD+E.R I showed here recently seems to be related to
>> the words 53-bits! and 80-bits! setting the mantissa size but assuming
>> the rounding mode is even. For some (bitrot) reason I changed the global
>> rounding mode at some point without realizing (or knowing) that it would
>> affect the double-double package. ISTR that setting the FPU rounding
>> mode in assembler was quite a hassle and decided that a shortcut would
>> be a good idea.
>>
>> -marcel
> 
> This set me thinking ...
> 
> : FROUND   nearest-mode rounding: (FROUND) ;
> : FLOOR   down-mode rounding: (FROUND) ;
> 
> Apparently I set truncate-mode  as default rounding.
> 
> \ Set  rounding-mode  for the duration of a definition.
> : rounding:   set-rounding-mode  CO  truncate-mode set-rounding-mode ;
> CODE (FROUND)   FRNDINT, NEXT, END-CODE
> 
> Why did I do that? Is nearest-mode not more appropriate?
> 
> default-mode 0 10 LSHIFT + CONSTANT nearest-mode
> default-mode 1 10 LSHIFT + CONSTANT down-mode
> default-mode 2 10 LSHIFT + CONSTANT up-mode
> default-mode 3 10 LSHIFT + CONSTANT truncate-mode
> 

Personally, I think rounding to nearest mode is a better default. In any 
case, for changing the rounding mode during a definition, it would be 
better if you save the current rounding mode before changing it and then 
restore it at the end of the definition instead of going back to your 
default setting. For example, if FROUND is called from a word which uses 
another rounding mode other than your default, you want to restore that 
mode when returning from FROUND.

It is more expensive to do so since you have to fetch the current fpu 
control word, but doing so would avoid a nasty surprise later.

--
KM

[toc] | [prev] | [next] | [standalone]

#135629

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-09 07:37 -0500
Message-ID<117rjv3$14ov3$1@dont-email.me>
In reply to#135545
On 9/3/26 6:03 PM, Krishna Myneni wrote:
> On 8/28/26 10:39, Krishna Myneni wrote:
...
>>>> 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.
>>>>
>> ...
>> One can examine the glibc code for double precision sin(x) at the link 
>> below.
>>
>> https://github.com/bminor/glibc/ 
>> blob/6059938728a98270b9706488887f43baa0471eba/sysdeps/ieee754/dbl-64/ 
>> s_sin.c
>>
>> === clip from above file ===
>> ...
>> /*******************************************************************/
>> /* An ultimate sin routine. Given an IEEE double machine number x  */
>> /* it computes the rounded value of sin(x).               */
>> /*******************************************************************/
>> #ifndef IN_SINCOS
>> double
>> SECTION
>> __sin (double x)
>> {
>> ...
>> === end of clip ===
>>
>> There are several ranges of the input parameter x, each of which are 
>> handled differently:
>>
>> I. 2^-26 < |x| < 0.855469;  no range reduction needed
>>
>> II. 0.855469  < |x| < 2.426265; ?
>>
>> III. 2.426265 < |x| < 105414350; uses function reduce_sincos()
>>
>> IV.  105414350 < |x| < 2^1024; uses function _branred()
>>
>> V. 105414350 < |x| < 2^1024; set error condition
>>
>>
...
> I've studied the branred.c code from glibc (version 2.4.2.9000, from 
> Github repo bminor/glibc), and I've made a stand-alone set of files 
> which compile to provide accurate, range-reduced sin(x) and cos(x) 
> functions -- yes, it works at x=1e100 rad!
> 
> The interesting part is that branred() uses double-double floating point 
> arithmetic (which provides a 106-bit significand, over an exponent range 
> the same as ordinary double-precision). Consequently, the glibc code for 
> sin(x) and cos(x) functions set the fpu precision and rounding mode to 
> double-precision and round to nearest, in order for double-double 
> precision arithmetic to work.
> ...

I have extracted the necessary files from glibc and modified them to 
make standalone object files for high accuracy sin(x) and cos(x) 
computation with full-range double-precision args and tested them under 
gcc linux x86 and Digital Mars C++ compilers.

kForth-Win32 will soon get an upgrade for its FSIN and FCOS words.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135631

Frompeter <peter.noreply@tin.it>
Date2026-09-09 18:49 +0200
Message-ID<20260909184929.0000509e@tin.it>
In reply to#135629
On Wed, 9 Sep 2026 07:37:55 -0500
Krishna Myneni <krishna.myneni@ccreweb.org> wrote:

> On 9/3/26 6:03 PM, Krishna Myneni wrote:
> > On 8/28/26 10:39, Krishna Myneni wrote:
> ...
> >>>> 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.
> >>>>
> >> ...
> >> One can examine the glibc code for double precision sin(x) at the link 
> >> below.
> >>
> >> https://github.com/bminor/glibc/ 
> >> blob/6059938728a98270b9706488887f43baa0471eba/sysdeps/ieee754/dbl-64/ 
> >> s_sin.c
> >>
> >> === clip from above file ===
> >> ...
> >> /*******************************************************************/
> >> /* An ultimate sin routine. Given an IEEE double machine number x  */
> >> /* it computes the rounded value of sin(x).               */
> >> /*******************************************************************/
> >> #ifndef IN_SINCOS
> >> double
> >> SECTION
> >> __sin (double x)
> >> {
> >> ...
> >> === end of clip ===
> >>
> >> There are several ranges of the input parameter x, each of which are 
> >> handled differently:
> >>
> >> I. 2^-26 < |x| < 0.855469;  no range reduction needed
> >>
> >> II. 0.855469  < |x| < 2.426265; ?
> >>
> >> III. 2.426265 < |x| < 105414350; uses function reduce_sincos()
> >>
> >> IV.  105414350 < |x| < 2^1024; uses function _branred()
> >>
> >> V. 105414350 < |x| < 2^1024; set error condition
> >>
> >>
> ...
> > I've studied the branred.c code from glibc (version 2.4.2.9000, from 
> > Github repo bminor/glibc), and I've made a stand-alone set of files 
> > which compile to provide accurate, range-reduced sin(x) and cos(x) 
> > functions -- yes, it works at x=1e100 rad!
> > 
> > The interesting part is that branred() uses double-double floating point 
> > arithmetic (which provides a 106-bit significand, over an exponent range 
> > the same as ordinary double-precision). Consequently, the glibc code for 
> > sin(x) and cos(x) functions set the fpu precision and rounding mode to 
> > double-precision and round to nearest, in order for double-double 
> > precision arithmetic to work.
> > ...
> 
> I have extracted the necessary files from glibc and modified them to 
> make standalone object files for high accuracy sin(x) and cos(x) 
> computation with full-range double-precision args and tested them under 
> gcc linux x86 and Digital Mars C++ compilers.
> 
> kForth-Win32 will soon get an upgrade for its FSIN and FCOS words.

An alternative is to use the files from the Core-Math project

https://gitlab.inria.fr/core-math/core-math

They are claimed to be correctly rounded, my tests also show this.
They are bloated and will compile to a large object file.
I have seen that also glibc has started to use selected files from them.

BR
Peter


> Krishna
>

[toc] | [prev] | [next] | [standalone]

#135633

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-09 15:06 -0500
Message-ID<117se80$1eerr$1@dont-email.me>
In reply to#135631
On 9/9/26 11:49 AM, peter wrote:
> On Wed, 9 Sep 2026 07:37:55 -0500
...
>> I have extracted the necessary files from glibc and modified them to
>> make standalone object files for high accuracy sin(x) and cos(x)
>> computation with full-range double-precision args and tested them under
>> gcc linux x86 and Digital Mars C++ compilers.
>>
>> kForth-Win32 will soon get an upgrade for its FSIN and FCOS words.
> 
> An alternative is to use the files from the Core-Math project
> 
> https://gitlab.inria.fr/core-math/core-math
> 
> They are claimed to be correctly rounded, my tests also show this.
> They are bloated and will compile to a large object file.
> I have seen that also glibc has started to use selected files from them.
> 
Below is a link to my modified files. I've included only the bare 
minimum to make it work under gcc and sc.exe (Digital Mars C/C++). The 
src/ directory contains the subdirs gcc/ and sc/ as well as the original 
files from glibc in orig/ .

The readme.txt files in gcc/ and sc/ provide build instructions, under 
those compilers, for a test program called test_sincos .

In the base folder, there is a comparison of sin(x) and cos(x) outputs 
over a range of angles for gcc and sc.

The object files are fairly small. Under sc,

$ ls -l *.obj
... 2533 Sep  9 08:49 branred.obj
... 3728 Sep  9 08:49 sincostab.obj
... 2315 Sep  9 08:49 s_sin.obj
...  546 Sep  9 08:54 test_sincos.obj

The standalone test program is linked with the C libraries,

... 64028 Sep  9 08:54 test_sincos.exe

--
KM

https://ccreweb.org/software/glibc/sincos/

[toc] | [prev] | [next] | [standalone]

#135655

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-10 19:03 -0500
Message-ID<117vgfk$2eh6r$1@dont-email.me>
In reply to#135631
On 9/9/26 11:49, peter wrote:
> On Wed, 9 Sep 2026 07:37:55 -0500
...
> 
> An alternative is to use the files from the Core-Math project
> 
> https://gitlab.inria.fr/core-math/core-math
> 
> They are claimed to be correctly rounded, my tests also show this.
> They are bloated and will compile to a large object file.
> I have seen that also glibc has started to use selected files from them.
> ...

The latest kForth-Win32 (v2.6.9) implements the glibc sin(x) and cos(x) 
functions. The words FSIN and FCOS call these functions which provide 
high accuracy results even at large double-precision angles. This makes 
consistent the accuracy of FSIN and FCOS across all kForth variants 
-32/64/Win32.

I have made a reference table of inputs and outputs for FCOS and FSIN 
over a range of angles. These may be useful for others who implement the 
same algorithm for FSIN and FCOS, for double precision floating point. 
The reference table may be found at

https://ccreweb.org/software/glibc/sincos/kForth-Win32_sincos_ref.txt

I would be interested to know if the core-math sin(x) and cos(x) 
functions give the same results for double precision.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135660

Frompeter <peter.noreply@tin.it>
Date2026-09-11 14:56 +0200
Message-ID<20260911145656.0000764d@tin.it>
In reply to#135655
On Thu, 10 Sep 2026 19:03:00 -0500
Krishna Myneni <krishna.myneni@ccreweb.org> wrote:

> On 9/9/26 11:49, peter wrote:
> > On Wed, 9 Sep 2026 07:37:55 -0500
> ...
> > 
> > An alternative is to use the files from the Core-Math project
> > 
> > https://gitlab.inria.fr/core-math/core-math
> > 
> > They are claimed to be correctly rounded, my tests also show this.
> > They are bloated and will compile to a large object file.
> > I have seen that also glibc has started to use selected files from them.
> > ...
> 
> The latest kForth-Win32 (v2.6.9) implements the glibc sin(x) and cos(x) 
> functions. The words FSIN and FCOS call these functions which provide 
> high accuracy results even at large double-precision angles. This makes 
> consistent the accuracy of FSIN and FCOS across all kForth variants 
> -32/64/Win32.
> 
> I have made a reference table of inputs and outputs for FCOS and FSIN 
> over a range of angles. These may be useful for others who implement the 
> same algorithm for FSIN and FCOS, for double precision floating point. 
> The reference table may be found at
> 
> https://ccreweb.org/software/glibc/sincos/kForth-Win32_sincos_ref.txt
> 
> I would be interested to know if the core-math sin(x) and cos(x) 
> functions give the same results for double precision.

I have collected close to 100000 tests for different functions.
They come from 2 sources

- the crmlib project (correctly rounded libm). This looks now to be dead.

- This site https://www.vinc17.net/research/testlibm/ I have taken the 
  worst cases file and transferred to Forth

I have put them on dropbox at this link 
https://www.dropbox.com/scl/fo/6j3md5qhll7sn4rtaf6mj/AFx3v5tB2Nwui9mb7Thw6oA?rlkey=savdo7h9gawbrhz3ctn9j63h5&st=59deulw0&dl=0

(I hope this long url works now when i copied it)

In the zip file there are 3 forth files ulptest.4 ulptest2.4 ulptest2h.4
The first one has the worst cases. The other 2 use the .dat files
and just show them in different way. These are based on crlibm data
All rounding should be set to nearest. They will work also on x87 80 bit
functions but you should set the fpu to operate at 53 bits. Expect most
trig functions to fail under x87 as they test large arguments.
All arguments and results are stored as 64 bit hex values. This avoids all
conversion errors but make it difficult to see what you actually test!

Here is the output from my system using core-math

include ulptest.4
Defined : f1/10** fnegate f10** ;
                 ulps error
Function    Total      0      1      2      3      4      5   6-10    >10   sign
facosh      1´877  1´877      0      0      0      0      0      0      0      0
facos       1´569  1´569      0      0      0      0      0      0      0      0
fasinh      2´262  2´262      0      0      0      0      0      0      0      0
fsinh       1´655  1´655      0      0      0      0      0      0      0      0
fatanh      1´552  1´552      0      0      0      0      0      0      0      0
fatan       1´735  1´735      0      0      0      0      0      0      0      0
fcbrt         138    138      0      0      0      0      0      0      0      0
fcosh       2´026  2´026      0      0      0      0      0      0      0      0
fcos        1´576  1´576      0      0      0      0      0      0      0      0
fcube         150    150      0      0      0      0      0      0      0      0
f1/10**       538    538      0      0      0      0      0      0      0      0
f10**       1´668  1´668      0      0      0      0      0      0      0      0
fexpm1      7´578  7´578      0      0      0      0      0      0      0      0
f2**        1´145  1´145      0      0      0      0      0      0      0      0
fexp        2´268  2´268      0      0      0      0      0      0      0      0
flog        1´883  1´883      0      0      0      0      0      0      0      0
flnp1       7´550  7´550      0      0      0      0      0      0      0      0
flg2          929    929      0      0      0      0      0      0      0      0
fln         2´813  2´813      0      0      0      0      0      0      0      0
fsinh       2´215  2´215      0      0      0      0      0      0      0      0
fsin        1´611  1´611      0      0      0      0      0      0      0      0
ftan        1´706  1´706      0      0      0      0      0      0      0      0
ftanh       1´852  1´852      0      0      0      0      0      0      0      0
frsqr       2´202  2´202      0      0      0      0      0      0      0      0
frsqrt      2´360  2´360      0      0      0      0      0      0      0      0
Total      52´858 52´858      0      0      0      0      0      0      0      0
ok

include ulptest2.4

                 ulps error
Function Total      0      1      2      3      4      5   6-10    >10   sign
sin     10´613 10´613      0      0      0      0      0      0      0      0
cos     10´790 10´790      0      0      0      0      0      0      0      0
tan      5´715  5´715      0      0      0      0      0      0      0      0
asin       726    726      0      0      0      0      0      0      0      0
acos       110    110      0      0      0      0      0      0      0      0
atan     5´618  5´618      0      0      0      0      0      0      0      0
sinh       416    416      0      0      0      0      0      0      0      0
cosh       549    549      0      0      0      0      0      0      0      0
pow     10´000 10´000      0      0      0      0      0      0      0      0
exp      4´282  4´282      0      0      0      0      0      0      0      0
expm1      238    238      0      0      0      0      0      0      0      0
ln       1´127  1´127      0      0      0      0      0      0      0      0
log10       52     52      0      0      0      0      0      0      0      0
logp1      227    227      0      0      0      0      0      0      0      0
To see records with a specific ulp difference run :
n ulpshow ulpxxx.dat
where n is the difference xxx the function sin asin etc
ok

If I switch  fdlibm instead about 50% of the test show 1 ulp error.
glibc libm is a bit better then that.

BR
Peter

> Krishna
>

[toc] | [prev] | [next] | [standalone]

#135661

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2026-09-11 09:17 -0500
Message-ID<11812i3$2v7oj$1@dont-email.me>
In reply to#135660
On 9/11/26 07:56, peter wrote:
> On Thu, 10 Sep 2026 19:03:00 -0500
> Krishna Myneni <krishna.myneni@ccreweb.org> wrote:
> 
>> On 9/9/26 11:49, peter wrote:
>>> On Wed, 9 Sep 2026 07:37:55 -0500
>> ...
>>>
>>> An alternative is to use the files from the Core-Math project
>>>
>>> https://gitlab.inria.fr/core-math/core-math
>>>
>>> They are claimed to be correctly rounded, my tests also show this.
>>> They are bloated and will compile to a large object file.
>>> I have seen that also glibc has started to use selected files from them.
>>> ...
>>
>> The latest kForth-Win32 (v2.6.9) implements the glibc sin(x) and cos(x)
>> functions. The words FSIN and FCOS call these functions which provide
>> high accuracy results even at large double-precision angles. This makes
>> consistent the accuracy of FSIN and FCOS across all kForth variants
>> -32/64/Win32.
>>
>> I have made a reference table of inputs and outputs for FCOS and FSIN
>> over a range of angles. These may be useful for others who implement the
>> same algorithm for FSIN and FCOS, for double precision floating point.
>> The reference table may be found at
>>
>> https://ccreweb.org/software/glibc/sincos/kForth-Win32_sincos_ref.txt
>>
>> I would be interested to know if the core-math sin(x) and cos(x)
>> functions give the same results for double precision.
> 
> I have collected close to 100000 tests for different functions.
> They come from 2 sources
> 
> - the crmlib project (correctly rounded libm). This looks now to be dead.
> 
> - This site https://www.vinc17.net/research/testlibm/ I have taken the
>    worst cases file and transferred to Forth
> 
> I have put them on dropbox at this link
> https://www.dropbox.com/scl/fo/6j3md5qhll7sn4rtaf6mj/AFx3v5tB2Nwui9mb7Thw6oA?rlkey=savdo7h9gawbrhz3ctn9j63h5&st=59deulw0&dl=0
> 
> (I hope this long url works now when i copied it)
> 
> In the zip file there are 3 forth files ulptest.4 ulptest2.4 ulptest2h.4
> The first one has the worst cases. The other 2 use the .dat files
> and just show them in different way. These are based on crlibm data
> All rounding should be set to nearest. They will work also on x87 80 bit
> functions but you should set the fpu to operate at 53 bits. Expect most
> trig functions to fail under x87 as they test large arguments.
> All arguments and results are stored as 64 bit hex values. This avoids all
> conversion errors but make it difficult to see what you actually test!
> 
> Here is the output from my system using core-math
> 
> include ulptest.4
> Defined : f1/10** fnegate f10** ;
>                   ulps error
> Function    Total      0      1      2      3      4      5   6-10    >10   sign
> facosh      1´877  1´877      0      0      0      0      0      0      0      0
> facos       1´569  1´569      0      0      0      0      0      0      0      0
> fasinh      2´262  2´262      0      0      0      0      0      0      0      0
> fsinh       1´655  1´655      0      0      0      0      0      0      0      0
> fatanh      1´552  1´552      0      0      0      0      0      0      0      0
> fatan       1´735  1´735      0      0      0      0      0      0      0      0
> fcbrt         138    138      0      0      0      0      0      0      0      0
> fcosh       2´026  2´026      0      0      0      0      0      0      0      0
> fcos        1´576  1´576      0      0      0      0      0      0      0      0
> fcube         150    150      0      0      0      0      0      0      0      0
> f1/10**       538    538      0      0      0      0      0      0      0      0
> f10**       1´668  1´668      0      0      0      0      0      0      0      0
> fexpm1      7´578  7´578      0      0      0      0      0      0      0      0
> f2**        1´145  1´145      0      0      0      0      0      0      0      0
> fexp        2´268  2´268      0      0      0      0      0      0      0      0
> flog        1´883  1´883      0      0      0      0      0      0      0      0
> flnp1       7´550  7´550      0      0      0      0      0      0      0      0
> flg2          929    929      0      0      0      0      0      0      0      0
> fln         2´813  2´813      0      0      0      0      0      0      0      0
> fsinh       2´215  2´215      0      0      0      0      0      0      0      0
> fsin        1´611  1´611      0      0      0      0      0      0      0      0
> ftan        1´706  1´706      0      0      0      0      0      0      0      0
> ftanh       1´852  1´852      0      0      0      0      0      0      0      0
> frsqr       2´202  2´202      0      0      0      0      0      0      0      0
> frsqrt      2´360  2´360      0      0      0      0      0      0      0      0
> Total      52´858 52´858      0      0      0      0      0      0      0      0
> ok
> 
> include ulptest2.4
> 
>                   ulps error
> Function Total      0      1      2      3      4      5   6-10    >10   sign
> sin     10´613 10´613      0      0      0      0      0      0      0      0
> cos     10´790 10´790      0      0      0      0      0      0      0      0
> tan      5´715  5´715      0      0      0      0      0      0      0      0
> asin       726    726      0      0      0      0      0      0      0      0
> acos       110    110      0      0      0      0      0      0      0      0
> atan     5´618  5´618      0      0      0      0      0      0      0      0
> sinh       416    416      0      0      0      0      0      0      0      0
> cosh       549    549      0      0      0      0      0      0      0      0
> pow     10´000 10´000      0      0      0      0      0      0      0      0
> exp      4´282  4´282      0      0      0      0      0      0      0      0
> expm1      238    238      0      0      0      0      0      0      0      0
> ln       1´127  1´127      0      0      0      0      0      0      0      0
> log10       52     52      0      0      0      0      0      0      0      0
> logp1      227    227      0      0      0      0      0      0      0      0
> To see records with a specific ulp difference run :
> n ulpshow ulpxxx.dat
> where n is the difference xxx the function sin asin etc
> ok
> 
> If I switch  fdlibm instead about 50% of the test show 1 ulp error.
> glibc libm is a bit better then that.
> 


Thanks, Peter!

I successfully ran ulptest.4th under kForth-64 by adding definitions of 
ON and OFF. The output is shown below, but I don't know how to interpret 
the different columns:

=== begin ===
$ kforth64-0.8.0
kForth-64 v 0.8.0	 (Build: 2026-05-03)
Copyright (c) 1998--2026 Krishna Myneni
Contributions by: dpw gd mu bk abs tn cmb bg dnw imss al
Provided under the GNU Affero General Public License, v3.0 or later


Ready!
include ulptest

/home/krishna/kforth/ans-words.4th

Defined : on true swap ! ;
Defined : off false swap ! ;
Defined : h. base @ >r hex . r> base ! ;
Defined : fcbrt 1e 3e f/    f** ;
Defined : fcube    3e       f** ;
Defined : f10**   falog         ;
Defined : f2**     2e fswap f** ;
Defined : f1/10** fnegate f10** ;
Defined : frsqrt  -0.5e     f** ;
Defined : frsqr   fdup f* 1e fswap f/ ;
Defined : flg2  flog 2e flog f/ ;

                  ulps error
Function    Total      0      1      2      3      4      5   6-10 
 >10   sign
facosh       1877    963    911      3      0      0      0      0 
0      0
facos        1569   1115    454      0      0      0      0      0 
0      0
fasinh       2262   1129   1133      0      0      0      0      0 
0      0
fsinh        1655    802    853      0      0      0      0      0 
0      0
fatanh       1552    758    794      0      0      0      0      0 
0      0
fatan        1735    963    772      0      0      0      0      0 
0      0
fcbrt         138     60     78      0      0      0      0      0 
0      0
fcosh        2026   1081    945      0      0      0      0      0 
0      0
fcos         1576    884    692      0      0      0      0      0 
0      0
fcube         150     79     71      0      0      0      0      0 
0      0
f1/10**       538    269    269      0      0      0      0      0 
0      0
f10**        1668    829    828     11      0      0      0      0 
0      0
fexpm1       7578   4220   3358      0      0      0      0      0 
0      0
f2**         1145    560    585      0      0      0      0      0 
0      0
fexp         2268   1142   1126      0      0      0      0      0 
0      0
flog         1883    991    891      1      0      0      0      0 
0      0
flnp1        7550   3912   3638      0      0      0      0      0 
0      0
flg2          929    432    479     18      0      0      0      0 
0      0
fln          2813   2035    778      0      0      0      0      0 
0      0
fsinh        2215   1087   1125      3      0      0      0      0 
0      0
fsin         1611    857    754      0      0      0      0      0 
0      0
ftan         1706    919    787      0      0      0      0      0 
0      0
ftanh        1852    929    913     10      0      0      0      0 
0      0
frsqr        2202   1959    243      0      0      0      0      0 
0      0
frsqrt       2360   2323     37      0      0      0      0      0 
0      0
Total       52858  30298  22514     46      0      0      0      0 
0      0
  ok
=== end ===

Btw, I had to convert file formats to unix before anything will load -- 
this is a deficiency in my Forth system. I would like it to be 
insensitive to text file format.

--
Krishna

[toc] | [prev] | [next] | [standalone]

#135663

Frompeter <peter.noreply@tin.it>
Date2026-09-11 17:04 +0200
Message-ID<20260911170435.000065b2@tin.it>
In reply to#135661
On Fri, 11 Sep 2026 09:17:37 -0500
Krishna Myneni <krishna.myneni@ccreweb.org> wrote:

> On 9/11/26 07:56, peter wrote:
> > On Thu, 10 Sep 2026 19:03:00 -0500
> > Krishna Myneni <krishna.myneni@ccreweb.org> wrote:
> > 
> >> On 9/9/26 11:49, peter wrote:
> >>> On Wed, 9 Sep 2026 07:37:55 -0500
> >> ...
> >>>
> >>> An alternative is to use the files from the Core-Math project
> >>>
> >>> https://gitlab.inria.fr/core-math/core-math
> >>>
> >>> They are claimed to be correctly rounded, my tests also show this.
> >>> They are bloated and will compile to a large object file.
> >>> I have seen that also glibc has started to use selected files from them.
> >>> ...
> >>
> >> The latest kForth-Win32 (v2.6.9) implements the glibc sin(x) and cos(x)
> >> functions. The words FSIN and FCOS call these functions which provide
> >> high accuracy results even at large double-precision angles. This makes
> >> consistent the accuracy of FSIN and FCOS across all kForth variants
> >> -32/64/Win32.
> >>
> >> I have made a reference table of inputs and outputs for FCOS and FSIN
> >> over a range of angles. These may be useful for others who implement the
> >> same algorithm for FSIN and FCOS, for double precision floating point.
> >> The reference table may be found at
> >>
> >> https://ccreweb.org/software/glibc/sincos/kForth-Win32_sincos_ref.txt
> >>
> >> I would be interested to know if the core-math sin(x) and cos(x)
> >> functions give the same results for double precision.
> > 
> > I have collected close to 100000 tests for different functions.
> > They come from 2 sources
> > 
> > - the crmlib project (correctly rounded libm). This looks now to be dead.
> > 
> > - This site https://www.vinc17.net/research/testlibm/ I have taken the
> >    worst cases file and transferred to Forth
> > 
> > I have put them on dropbox at this link
> > https://www.dropbox.com/scl/fo/6j3md5qhll7sn4rtaf6mj/AFx3v5tB2Nwui9mb7Thw6oA?rlkey=savdo7h9gawbrhz3ctn9j63h5&st=59deulw0&dl=0
> > 
> > (I hope this long url works now when i copied it)
> > 
> > In the zip file there are 3 forth files ulptest.4 ulptest2.4 ulptest2h.4
> > The first one has the worst cases. The other 2 use the .dat files
> > and just show them in different way. These are based on crlibm data
> > All rounding should be set to nearest. They will work also on x87 80 bit
> > functions but you should set the fpu to operate at 53 bits. Expect most
> > trig functions to fail under x87 as they test large arguments.
> > All arguments and results are stored as 64 bit hex values. This avoids all
> > conversion errors but make it difficult to see what you actually test!
> > 
> > Here is the output from my system using core-math
> > 
> > include ulptest.4
> > Defined : f1/10** fnegate f10** ;
> >                   ulps error
> > Function    Total      0      1      2      3      4      5   6-10    >10   sign
> > facosh      1´877  1´877      0      0      0      0      0      0      0      0
> > facos       1´569  1´569      0      0      0      0      0      0      0      0
> > fasinh      2´262  2´262      0      0      0      0      0      0      0      0
> > fsinh       1´655  1´655      0      0      0      0      0      0      0      0
> > fatanh      1´552  1´552      0      0      0      0      0      0      0      0
> > fatan       1´735  1´735      0      0      0      0      0      0      0      0
> > fcbrt         138    138      0      0      0      0      0      0      0      0
> > fcosh       2´026  2´026      0      0      0      0      0      0      0      0
> > fcos        1´576  1´576      0      0      0      0      0      0      0      0
> > fcube         150    150      0      0      0      0      0      0      0      0
> > f1/10**       538    538      0      0      0      0      0      0      0      0
> > f10**       1´668  1´668      0      0      0      0      0      0      0      0
> > fexpm1      7´578  7´578      0      0      0      0      0      0      0      0
> > f2**        1´145  1´145      0      0      0      0      0      0      0      0
> > fexp        2´268  2´268      0      0      0      0      0      0      0      0
> > flog        1´883  1´883      0      0      0      0      0      0      0      0
> > flnp1       7´550  7´550      0      0      0      0      0      0      0      0
> > flg2          929    929      0      0      0      0      0      0      0      0
> > fln         2´813  2´813      0      0      0      0      0      0      0      0
> > fsinh       2´215  2´215      0      0      0      0      0      0      0      0
> > fsin        1´611  1´611      0      0      0      0      0      0      0      0
> > ftan        1´706  1´706      0      0      0      0      0      0      0      0
> > ftanh       1´852  1´852      0      0      0      0      0      0      0      0
> > frsqr       2´202  2´202      0      0      0      0      0      0      0      0
> > frsqrt      2´360  2´360      0      0      0      0      0      0      0      0
> > Total      52´858 52´858      0      0      0      0      0      0      0      0
> > ok
> > 
> > include ulptest2.4
> > 
> >                   ulps error
> > Function Total      0      1      2      3      4      5   6-10    >10   sign
> > sin     10´613 10´613      0      0      0      0      0      0      0      0
> > cos     10´790 10´790      0      0      0      0      0      0      0      0
> > tan      5´715  5´715      0      0      0      0      0      0      0      0
> > asin       726    726      0      0      0      0      0      0      0      0
> > acos       110    110      0      0      0      0      0      0      0      0
> > atan     5´618  5´618      0      0      0      0      0      0      0      0
> > sinh       416    416      0      0      0      0      0      0      0      0
> > cosh       549    549      0      0      0      0      0      0      0      0
> > pow     10´000 10´000      0      0      0      0      0      0      0      0
> > exp      4´282  4´282      0      0      0      0      0      0      0      0
> > expm1      238    238      0      0      0      0      0      0      0      0
> > ln       1´127  1´127      0      0      0      0      0      0      0      0
> > log10       52     52      0      0      0      0      0      0      0      0
> > logp1      227    227      0      0      0      0      0      0      0      0
> > To see records with a specific ulp difference run :
> > n ulpshow ulpxxx.dat
> > where n is the difference xxx the function sin asin etc
> > ok
> > 
> > If I switch  fdlibm instead about 50% of the test show 1 ulp error.
> > glibc libm is a bit better then that.
> > 
> 
> 
> Thanks, Peter!
> 
> I successfully ran ulptest.4th under kForth-64 by adding definitions of 
> ON and OFF. The output is shown below, but I don't know how to interpret 
> the different columns:
> 
> === begin ===
> $ kforth64-0.8.0
> kForth-64 v 0.8.0	 (Build: 2026-05-03)
> Copyright (c) 1998--2026 Krishna Myneni
> Contributions by: dpw gd mu bk abs tn cmb bg dnw imss al
> Provided under the GNU Affero General Public License, v3.0 or later
> 
> 
> Ready!
> include ulptest
> 
> /home/krishna/kforth/ans-words.4th
> 
> Defined : on true swap ! ;
> Defined : off false swap ! ;
> Defined : h. base @ >r hex . r> base ! ;
> Defined : fcbrt 1e 3e f/    f** ;
> Defined : fcube    3e       f** ;
> Defined : f10**   falog         ;
> Defined : f2**     2e fswap f** ;
> Defined : f1/10** fnegate f10** ;
> Defined : frsqrt  -0.5e     f** ;
> Defined : frsqr   fdup f* 1e fswap f/ ;
> Defined : flg2  flog 2e flog f/ ;
> 
>                   ulps error
> Function    Total      0      1      2      3      4      5   6-10 
>  >10   sign
> facosh       1877    963    911      3      0      0      0      0 
> 0      0
> facos        1569   1115    454      0      0      0      0      0 
> 0      0
> fasinh       2262   1129   1133      0      0      0      0      0 
> 0      0
> fsinh        1655    802    853      0      0      0      0      0 
> 0      0
> fatanh       1552    758    794      0      0      0      0      0 
> 0      0
> fatan        1735    963    772      0      0      0      0      0 
> 0      0
> fcbrt         138     60     78      0      0      0      0      0 
> 0      0
> fcosh        2026   1081    945      0      0      0      0      0 
> 0      0
> fcos         1576    884    692      0      0      0      0      0 
> 0      0
> fcube         150     79     71      0      0      0      0      0 
> 0      0
> f1/10**       538    269    269      0      0      0      0      0 
> 0      0
> f10**        1668    829    828     11      0      0      0      0 
> 0      0
> fexpm1       7578   4220   3358      0      0      0      0      0 
> 0      0
> f2**         1145    560    585      0      0      0      0      0 
> 0      0
> fexp         2268   1142   1126      0      0      0      0      0 
> 0      0
> flog         1883    991    891      1      0      0      0      0 
> 0      0
> flnp1        7550   3912   3638      0      0      0      0      0 
> 0      0
> flg2          929    432    479     18      0      0      0      0 
> 0      0
> fln          2813   2035    778      0      0      0      0      0 
> 0      0
> fsinh        2215   1087   1125      3      0      0      0      0 
> 0      0
> fsin         1611    857    754      0      0      0      0      0 
> 0      0
> ftan         1706    919    787      0      0      0      0      0 
> 0      0
> ftanh        1852    929    913     10      0      0      0      0 
> 0      0
> frsqr        2202   1959    243      0      0      0      0      0 
> 0      0
> frsqrt       2360   2323     37      0      0      0      0      0 
> 0      0
> Total       52858  30298  22514     46      0      0      0      0 
> 0      0
>   ok
> === end ===
> 
> Btw, I had to convert file formats to unix before anything will load -- 
> this is a deficiency in my Forth system. I would like it to be 
> insensitive to text file format.
> 
> --
> Krishna
> 

If you look at the total, 52858  30298  22514     46 

you ran 52858 tests
30298 was correct
22514 had an error in the last bit, 1 ulp error
46 had 2 ulps error.
This is a good result and I could guess that you link to libm in glibc
ulps is calculated as abs(correct-calculated) doing the calculation
on the 64 bit integer representation of the double float.
I basically do fsin tmp df! tmp @ correct - 

If you run the other file ulptest2.4 you will have more trig functions
tested. There are 2 files as they come from different sources

BR
Peter

[toc] | [prev] | [next] | [standalone]

Page 1 of 2  [1] 2  Next page →

Back to top | Article view | comp.lang.forth

csiph-web