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


Groups > comp.lang.forth > #135663

Re: Range Reduction Using Big Number Arithmetic

From peter <peter.noreply@tin.it>
Newsgroups comp.lang.forth
Subject Re: Range Reduction Using Big Number Arithmetic
Date 2026-09-11 17:04 +0200
Organization A noiseless patient Spider
Message-ID <20260911170435.000065b2@tin.it> (permalink)
References (4 earlier) <117rjv3$14ov3$1@dont-email.me> <20260909184929.0000509e@tin.it> <117vgfk$2eh6r$1@dont-email.me> <20260911145656.0000764d@tin.it> <11812i3$2v7oj$1@dont-email.me>

Show all headers | View raw


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

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


Thread

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

csiph-web