Path: csiph.com!eternal-september.org!feeder.eternal-september.org!nntp.eternal-september.org!.POSTED!not-for-mail From: peter Newsgroups: comp.lang.forth Subject: Re: Range Reduction Using Big Number Arithmetic Date: Fri, 11 Sep 2026 14:56:56 +0200 Organization: A noiseless patient Spider Lines: 165 Message-ID: <20260911145656.0000764d@tin.it> References: <115pm90$3cl9h$1@dont-email.me> <1168f21$bnk$1@dont-email.me> <116sa3b$2dpj8$4@dont-email.me> <117cuca$4f48$1@dont-email.me> <117rjv3$14ov3$1@dont-email.me> <20260909184929.0000509e@tin.it> <117vgfk$2eh6r$1@dont-email.me> MIME-Version: 1.0 Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Injection-Date: Fri, 11 Sep 2026 12:56:58 +0000 (UTC) Injection-Info: dont-email.me; logging-data="2931655"; mail-complaints-to="abuse@eternal-september.org"; posting-account="U2FsdGVkX1/zWsPsWadMahQoXTowbhB2CQb5ciLkFDk="; posting-host="8fd06948898d777b3b4c9a856d230eaf" Cancel-Lock: sha1:F+FhSk1U2bT9X6UJjEn65d4qn/s= sha256:aOjPSRstEF4XiFGdCO5vt8o8VpKVnhTok/DxiB5kQlk= sha1:uV+d9WLuFsfnLb07zqe9pPcZM2M= sha256:z0IgjDDE4t271+hvMR3yjQ1csE2BB9rvxJ/3LhGTz4o= X-Newsreader: Claws Mail 4.4.0 (GTK 3.24.51; x86_64-w64-mingw32) Xref: csiph.com comp.lang.forth:135660 On Thu, 10 Sep 2026 19:03:00 -0500 Krishna Myneni wrote: > On 9/9/26 11:49, peter wrote: > > On Wed, 9 Sep 2026 07:37:55 -0500 > ... > >=20 > > An alternative is to use the files from the Core-Math project > >=20 > > https://gitlab.inria.fr/core-math/core-math > >=20 > > 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. > > ... >=20 > The latest kForth-Win32 (v2.6.9) implements the glibc sin(x) and cos(x)=20 > functions. The words FSIN and FCOS call these functions which provide=20 > high accuracy results even at large double-precision angles. This makes=20 > consistent the accuracy of FSIN and FCOS across all kForth variants=20 > -32/64/Win32. >=20 > I have made a reference table of inputs and outputs for FCOS and FSIN=20 > over a range of angles. These may be useful for others who implement the= =20 > same algorithm for FSIN and FCOS, for double precision floating point.=20 > The reference table may be found at >=20 > https://ccreweb.org/software/glibc/sincos/kForth-Win32_sincos_ref.txt >=20 > I would be interested to know if the core-math sin(x) and cos(x)=20 > 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=20 worst cases file and transferred to Forth I have put them on dropbox at this link=20 https://www.dropbox.com/scl/fo/6j3md5qhll7sn4rtaf6mj/AFx3v5tB2Nwui9mb7Thw6o= A?rlkey=3Dsavdo7h9gawbrhz3ctn9j63h5&st=3D59deulw0&dl=3D0 (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=B4877 1=B4877 0 0 0 0 0 0 = 0 0 facos 1=B4569 1=B4569 0 0 0 0 0 0 = 0 0 fasinh 2=B4262 2=B4262 0 0 0 0 0 0 = 0 0 fsinh 1=B4655 1=B4655 0 0 0 0 0 0 = 0 0 fatanh 1=B4552 1=B4552 0 0 0 0 0 0 = 0 0 fatan 1=B4735 1=B4735 0 0 0 0 0 0 = 0 0 fcbrt 138 138 0 0 0 0 0 0 0 = 0 fcosh 2=B4026 2=B4026 0 0 0 0 0 0 = 0 0 fcos 1=B4576 1=B4576 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=B4668 1=B4668 0 0 0 0 0 0 = 0 0 fexpm1 7=B4578 7=B4578 0 0 0 0 0 0 = 0 0 f2** 1=B4145 1=B4145 0 0 0 0 0 0 = 0 0 fexp 2=B4268 2=B4268 0 0 0 0 0 0 = 0 0 flog 1=B4883 1=B4883 0 0 0 0 0 0 = 0 0 flnp1 7=B4550 7=B4550 0 0 0 0 0 0 = 0 0 flg2 929 929 0 0 0 0 0 0 0 = 0 fln 2=B4813 2=B4813 0 0 0 0 0 0 = 0 0 fsinh 2=B4215 2=B4215 0 0 0 0 0 0 = 0 0 fsin 1=B4611 1=B4611 0 0 0 0 0 0 = 0 0 ftan 1=B4706 1=B4706 0 0 0 0 0 0 = 0 0 ftanh 1=B4852 1=B4852 0 0 0 0 0 0 = 0 0 frsqr 2=B4202 2=B4202 0 0 0 0 0 0 = 0 0 frsqrt 2=B4360 2=B4360 0 0 0 0 0 0 = 0 0 Total 52=B4858 52=B4858 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 si= gn sin 10=B4613 10=B4613 0 0 0 0 0 0 0 = 0 cos 10=B4790 10=B4790 0 0 0 0 0 0 0 = 0 tan 5=B4715 5=B4715 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=B4618 5=B4618 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=B4000 10=B4000 0 0 0 0 0 0 0 = 0 exp 4=B4282 4=B4282 0 0 0 0 0 0 0 = 0 expm1 238 238 0 0 0 0 0 0 0 = 0 ln 1=B4127 1=B4127 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 >=20