Re: Range Reduction Using Big Number Arithmetic

From	Krishna Myneni <krishna.myneni@ccreweb.org>
Newsgroups	comp.lang.forth
Subject	Re: Range Reduction Using Big Number Arithmetic
Date	2026-09-03 18:03 -0500
Organization	A noiseless patient Spider
Message-ID	<117cuca$4f48$1@dont-email.me> (permalink)
References	<115pm90$3cl9h$1@dont-email.me> <1168f21$bnk$1@dont-email.me> <116sa3b$2dpj8$4@dont-email.me>

Show all headers | View raw

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

Thread

Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-15 07:32 -0500
  Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-15 07:48 -0500
  Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-20 22:01 -0500
    Re: Range Reduction Using Big Number Arithmetic albert@spenarnc.xs4all.nl - 2026-08-22 19:31 +0200
      Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-31 08:49 -0500
    Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-28 10:39 -0500
      Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-08-28 10:42 -0500
      Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-03 18:03 -0500
        Re: Range Reduction Using Big Number Arithmetic marcel hendrix <mhx@iae.nl> - 2026-09-04 07:29 +0200
          Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-04 07:45 -0500
            Re: Range Reduction Using Big Number Arithmetic marcel hendrix <mhx@iae.nl> - 2026-09-04 15:29 +0200
              Re: Range Reduction Using Big Number Arithmetic albert@spenarnc.xs4all.nl - 2026-09-05 12:27 +0200
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-05 08:39 -0500
        Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-09 07:37 -0500
          Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-09 18:49 +0200
            Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-09 15:06 -0500
            Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-10 19:03 -0500
              Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 14:56 +0200
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 09:17 -0500
                Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 17:04 +0200
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 10:46 -0500
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 13:24 -0500
                Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-11 23:21 +0200
                Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-12 00:42 +0200
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 18:33 -0500
                Re: Range Reduction Using Big Number Arithmetic peter <peter.noreply@tin.it> - 2026-09-12 11:04 +0200
                Re: Range Reduction Using Big Number Arithmetic Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-12 10:33 -0500

csiph-web