Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #135675
| From | peter <peter.noreply@tin.it> |
|---|---|
| Newsgroups | comp.lang.forth |
| Subject | Re: Range Reduction Using Big Number Arithmetic |
| Date | 2026-09-12 00:42 +0200 |
| Organization | A noiseless patient Spider |
| Message-ID | <20260912004252.00007f04@tin.it> (permalink) |
| References | (5 earlier) <20260909184929.0000509e@tin.it> <117vgfk$2eh6r$1@dont-email.me> <20260911145656.0000764d@tin.it> <1181h08$34q9g$1@dont-email.me> <20260911232156.00004f2b@tin.it> |
On Fri, 11 Sep 2026 23:21:56 +0200
peter <peter.noreply@tin.it> wrote:
> On Fri, 11 Sep 2026 13:24:06 -0500
> Krishna Myneni <krishna.myneni@ccreweb.org> wrote:
>
> > On 9/11/26 7:56 AM, 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.
> > >
> >
> > In ulptest.4th:
> >
> > The 32-bit T2 test does not work the same way as the 64-bit T2 test.
> > Specifically, the 64-bit T2 test checks for same sign or nans while the
> > 32-bit T2 tests only checks for same sign.
> >
> > I think there is also somewhere in the 32-bit section where the order of
> > the bits is not correct, because I'm getting a lot of >10 ulp errors.
> > When I turn verbose on, I see that for some tests the byte order of the
> > 32 msbits and the 32 lsbits of the result are reversed. It does not
> > occur for every test index number.
>
> I will have to look into this. I have not used the 32 bit version much
> the last 10 years! so probably I have forgotten to keep 64 and 32 bit
> aligned in implementation. my 32 bit lxf/ntf uses the 80 bits float.
> They will have a larger number of problems and I have not notised
>
> Peter
>
> >
> > Example:
> > \ Single line is broken into multiple lines for readability
> > 52245 tanh $3F062B4331584D27. T1 ftanh $3F062B43311F8E20. T2
> > \ 0x1.62b4331584d27p-15 0x1.62b43311f8e2p-15
> > \ Got $311F8E203F062B43
> >
> >
> >
> > Under kForth-32 for Linux, I obtain the following:
> >
> > === begin output ===
> > $ kforth32-2.8.0
> > kForth-32 v 2.8.0 (Build: 2026-05-02)
> > Copyright (c) 1998--2026 Krishna Myneni
> > Contributions by: dpw gd mu bk abs tn cmb bg dnw
> > 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 947 0 0 0 0 0 0 930
> > 0
> > facos 1569 784 0 0 0 0 0 0 785
> > 0
> > fasinh 2262 1130 0 0 0 0 0 0 132
> > 0
> > fsinh 1655 832 0 0 0 0 0 0 82 0
> > fatanh 1552 748 0 0 0 0 0 0 803
> > 1
> > fatan 1735 869 0 0 0 0 0 0 865
> > 1
> > fcbrt 138 58 0 0 0 0 0 0 80
> > 0
> > fcosh 2026 947 0 0 0 0 0 0 079
> > 0
> > fcos 1576 889 0 0 0 0 0 0 687
> > 0
> > fcube 150 73 0 0 0 0 0 0 77
> > 0
> > f1/10** 538 86 0 0 0 0 0 0 452
> > 0
> > f10** 1668 600 0 0 0 0 0 0 068
> > 0
> > fexpm1 7578 3770 0 0 0 0 0 0 807
> > 1
> > f2** 1145 776 0 0 0 0 0 0 369
> > 0
> > fexp 2268 716 0 0 0 0 0 0 552
> > 0
> > flog 1883 1014 0 0 0 0 0 0 869
> > 0
> > flnp1 7550 3612 0 0 0 0 0 0 936
> > 2
> > flg2 929 461 0 0 0 0 0 0 468
> > 0
> > fln 2813 1598 0 0 0 0 0 0 215
> > 0
> > fsinh 2215 840 0 0 0 0 0 0 375
> > 0
> > fsin 1611 859 0 0 0 0 0 0 752
> > 0
> > ftan 1706 893 0 0 0 0 0 0 813
> > 0
> > ftanh 1852 901 0 0 0 0 0 0 951
> > 0
> > frsqr 2202 1959 0 0 0 0 0 0 243
> > 0
> > frsqrt 2360 2262 0 0 0 0 0 0 98
> > 0
> > Total 52858 27624 0 0 0 0 0 0 5229
> > 5
> > ok
> > === end output ===
> >
> > --
> > KM
I have looked at the file and can not find an obvious error in swaping
high and low parts. But I can see also in your example that they are swaped
when printed in the .error function.
I also tried it on the 32 bit ntf and get
ulps error
Function Total 0 1 2 3 4 5 6-10 >10 sign
facosh 1877 401 431 23 19 15 8 28 952 0
facos 1569 255 313 34 14 6 7 25 915 0
fasinh 2262 725 766 65 36 27 17 57 569 0
fsinh 1655 848 807 0 0 0 0 0 0 0
fatanh 1552 778 774 0 0 0 0 0 0 0
fatan 1735 869 866 0 0 0 0 0 0 0
fcbrt 138 70 68 0 0 0 0 0 0 0
fcosh 2026 1018 1008 0 0 0 0 0 0 0
fcos 1576 806 770 0 0 0 0 0 0 0
fcube 150 76 74 0 0 0 0 0 0 0
f1/10** 538 266 272 0 0 0 0 0 0 0
f10** 1668 811 857 0 0 0 0 0 0 0
fexpm1 7578 3779 3799 0 0 0 0 0 0 0
f2** 1145 592 553 0 0 0 0 0 0 0
fexp 2268 1150 1118 0 0 0 0 0 0 0
flog 1883 997 886 0 0 0 0 0 0 0
flnp1 7550 3623 3927 0 0 0 0 0 0 0
flg2 929 526 403 0 0 0 0 0 0 0
fln 2813 1598 1215 0 0 0 0 0 0 0
fsinh 2215 1120 1095 0 0 0 0 0 0 0
fsin 1611 790 821 0 0 0 0 0 0 0
ftan 1706 816 890 0 0 0 0 0 0 0
ftanh 1852 916 936 0 0 0 0 0 0 0
frsqr 2202 2169 33 0 0 0 0 0 0 0
frsqrt 2360 2318 42 0 0 0 0 0 0 0
Total 52858 27317 22724 122 69 48 32 110 2436 0
ok
I had to redefine f2** and flg2 that are present in ntf but with completly
different meaning. the large number of problems in facosh, facos and fasinh
are not unexpected. They have very simplistic implementations.
Are you using glibc also on the 32 bit Kforth. Is it using the FPU or
other algorithms?
I got this for the same line
52545 tanh $3F28EE5F1F66AC0A. T1
ftanh $3F28EE5F1A5B5193. T2 \ 0x1.8ee5f1f66ac0ap-13
0x1.8ee5f1a5b5193p-13 Got $3F28EE5F1A5B5193
Wich is strange as they look identical and should not have been signaled
I see now that .error tests the ulp variable but it is never set anywhere
I will have to fix that
BR
Peter
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll 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