-
Notifications
You must be signed in to change notification settings - Fork 138
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Accuracy of single precision functions #212
Comments
While the max ULP is an important metric, two additional metrics may be more important: speed and distribution of ULP (if max ULP > 1). Consider % tlibm -DfEsP -x 1 -X 4 j0 Interval tested for j0f: [1,4] 1000000 calls, 0.090046 secs, 0.09005 usecs/call ulp <= 0.5: 49.102% 8237966 | 49.102% 8237966 0.5 < ulp < 0.6: 7.534% 1263940 | 56.636% 9501906 0.6 < ulp < 0.7: 6.538% 1096871 | 63.174% 10598777 0.7 < ulp < 0.8: 5.690% 954682 | 68.864% 11553459 0.8 < ulp < 0.9: 4.883% 819227 | 73.747% 12372686 0.9 < ulp < 1.0: 4.015% 673617 | 77.762% 13046303 1.0 < ulp < 1.5: 12.171% 2042032 | 89.933% 15088335 1.5 < ulp < 2.0: 5.488% 920737 | 95.422% 16009072 2.0 < ulp < 3.0: 3.391% 568873 | 98.812% 16577945 3.0 < ulp < : 1.188% 199271 | 100.000% 16777216 Max ulp: 900690.687500 at 2.40482545e+00 The reference function was % ./tlibm -DfEsP -x 1 -X 4 j0 Interval tested for j0f: [1,4] 1000000 calls, 0.045775 secs, 0.04578 usecs/call ulp <= 0.5: 57.112% 9581798 | 57.112% 9581798 0.5 < ulp < 0.6: 7.607% 1276254 | 64.719% 10858052 0.6 < ulp < 0.7: 5.920% 993139 | 70.639% 11851191 0.7 < ulp < 0.8: 4.656% 781177 | 75.295% 12632368 0.8 < ulp < 0.9: 3.957% 663922 | 79.252% 13296290 0.9 < ulp < 1.0: 3.321% 557242 | 82.573% 13853532 1.0 < ulp < 1.5: 10.080% 1691123 | 92.653% 15544655 1.5 < ulp < 2.0: 4.563% 765544 | 97.216% 16310199 2.0 < ulp < 3.0: 2.529% 424342 | 99.746% 16734541 3.0 < ulp < : 0.254% 42675 | 100.000% 16777216 Max ulp: 5.246863 at 3.98081970e+00 A much more accurate approximation that is twice as fast! Unfortunately, Bessel functions are not periodic so argument reduction cannot be used. Thus, one needs multiple polynomial approximations to cover a significant domain. J. Harrison [1] proposed an algorithm that uses multiple polynomial for |
Dear Steve, on my side I am more interested by the max ULP. I know other people are more interested by speed. I guess the Intel Math Library implements the algorithm from John Harrison. |
I have updated my results with OpenLibm 0.7.0. The results are almost the same, except for y1f whose largest error has improved. The url is still the same. Previous versions are available from https://members.loria.fr/PZimmermann/papers/ |
I have revised my note with double and quadruple precision, and bivariate functions like atan2, hypot, pow: https://homepages.loria.fr/PZimmermann/papers/#accuracy For double precision, the largest error I found for OpenLibm 0.7.0 is less than 4 ulps, except for the Bessel functions, and the lgamma and tgamma functions. If something has changed in newer versions of OpenLibm, please tell me and I will update my note. |
a new version of my note is now available, including the "long double" format: https://homepages.loria.fr/PZimmermann/papers/#accuracy I found a few issues with Openlibm, which I reported in separate new issues here. |
I wonder what the easiest way is to fix some of these. Maybe just pick the good implementations from other libm libraries and bring them in here? |
I am about to release of new version of the accuracy comparison. It will be with Openlibm 0.7.5, the last release. For single precision, the large error for powf has been fixed. It remains large errors for the Bessel functions (j0f, j1f, y0f, y1f) and lgammaf. All other functions have a maximal error less than 3.17 ulps. |
I said above that in OpenLibm 0.7.5, for single precision, the large error for powf has been fixed. But I have found a new pair of inputs giving a huge error: for x=0x1.ffffecp-1 and y=-0x1.000002p+27 OpenLibm 0.7.5 yields +Inf, whereas the correct result is 0x1.557a86p115 (rounding to nearest). |
This patch fixes the issue reported by Paul Zimmerman.
|
Summary: From Steve Kargl: Paul Zimmermann has identified a bug in Openlibm's powf(), which is identical to FreeBSD's libm. Both derived from fdlibm. JuliaMath/openlibm#212. Consider % cat h.c int main(void) { float x, y, z; x = 0x1.ffffecp-1F; y = -0x1.000002p+27F; z = 0x1.557a86p115F; printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z); return 0; } % cc -o h -fno-builtin h.c -lm && ./h 9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34 Reviewers: manu Subscribers: imp, andrew, emaste Differential Revision: https://reviews.freebsd.org/D31865
Fix in a090d3e |
Summary: From Steve Kargl: Paul Zimmermann has identified a bug in Openlibm's powf(), which is identical to FreeBSD's libm. Both derived from fdlibm. JuliaMath/openlibm#212. Consider % cat h.c int main(void) { float x, y, z; x = 0x1.ffffecp-1F; y = -0x1.000002p+27F; z = 0x1.557a86p115F; printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z); return 0; } % cc -o h -fno-builtin h.c -lm && ./h 9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34 Reviewers: manu Subscribers: imp, andrew, emaste Differential Revision: https://reviews.freebsd.org/D31865
within the CORE-MATH project, we now have a tool to measure the efficiency of math library functions. Here is a comparison between the CORE-MATH routines and OpenLibm: https://sympa.inria.fr/sympa/arc/core-math/2022-03/msg00011.html. The CORE-MATH rountines are under MIT license, and should be easy to integrate in OpenLibm. |
here are updated timings for the single precision functions, with revision d14da3d of CORE-MATH, on an Intel i5-4590:
First column is CORE-MATH's reciprocal throughput, 2nd is for GNU libc, and 3rd for OpenLibm. |
ae2d916 Correctly round double precision sqrt (#256) 81d5e16 Add fmod assembly version (#255) 465ca0a Update README.md 428e7af Support for riscv64 architecture (#254) ed7aea3 Bump version to 0.8 (#248) 69bb280 Another Windows ARM64 fix (#253) 3d4a902 Fixes for Windows ARM64 (#251) a9568fb [Windows] install import library to libdir (#249) f88e52a CI (Windows): set `msys2 {0}` as the default shell for all Windows steps (#247) b48a2f7 CI (Linux and macOS): Remove the `arch` variable, which currently has no effect (#246) 2a47fa5 CI: A variety of miscellaneous tweaks (#244) d0ef09a prefix symbols with _ for 32-bit x86 Windows (#242) 60dec83 msys2 ci (#243) 6ea5fa2 Merge pull request #240 from JuliaMath/vs/msys 437c139 Update ci.yml e993267 Update ci.yml 4a36c50 Update ci.yml 24cec17 Update ci.yml d26ed98 Update ci.yml 7b96025 Update ci.yml 7af65db Update ci.yml a2e053e Revert "Update ci.yml" 4a52bb0 Update ci.yml fb10fcf Update ci.yml abf5aaa Update ci.yml 98dcc48 Update ci.yml ff822f3 Update ci.yml ab8d1ad Update ci.yml 4d97e2d Update ci.yml 72caeab Update ci.yml 9dd3049 Create ci.yml 15119bc Merge pull request #239 from JuliaMath/revert-238-patch-1 4bca0f2 Revert "prefix symbols with _ for 32-bit x86 Windows" 3b9454f Merge pull request #238 from jeremyd2019/patch-1 6ae6318 Update src/cdefs-compat.h 71a8fd1 Merge pull request #233 from lephe/more-long-double-aliases 7a3ef59 prefix symbols with _ for 32-bit x86 Windows a871457 Merge pull request #230 from PetteriAimonen/master a090d3e Fix powf: JuliaMath/openlibm#212 (comment) 57dd0f9 add missing weak references for long double functions 327b1bd Replace remaining __strong_alias uses f052f42 Merge pull request #228 from JuliaMath/aa/hypotl 711654e Fix incorrect results in `hypotl` near underflow aeab19f Fix for #211 Co-authored by: @kargl 5449705 Merge pull request #227 from JuliaMath/vs/powf 98f8713 Fix #211 Patched by importing latest msun version 6a85b33 Merge pull request #225 from JuliaMath/vs/strict_assign 40dac9d Restore STRICT_ASSIGN on FreeBSD as suggested in #215 2d10c90 Merge pull request #218 from jcestibariz/fix-wasm32 5d70ac5 Merge pull request #221 from maleadt/tb/static_fenv 63aa875 Make fenv methods static on additional platforms. 9152b0d Fix compilation errors on wasm32 3cb8045 Merge pull request #217 from epsilon-0/master c856101 Merge pull request #219 from maleadt/tb/dont_export_fenv be31bff Revert "Export `fenv` functions on all platforms (#213)" eb21e8a don't alter toolchain vars if already provided b34f107 Fix Apple Silicon build (#214) 5a27b4c Export `fenv` functions on all platforms (#213) 878948d Update list of libm libraries 508603d Update index.html 0276147 Merge pull request #209 from embeddedartistry/master f2a8b36 Update download links to point to releases 1d6befd Merge pull request #208 from cndesantana/patch-1 4f559d4 Replace a few remaining __strong_reference uses (#210) 0418324 Refactor: OLM_DLLEXPORT definition now lives in a standalone header. 861b2ad Fix small typo 5b0e7e9 Update FUNDING.yml 382b8e9 Add musl-libc math library to resources. f6ad75a update openlibm website. 5efed30 Bump SONAME as discussed in #200 e6ac7d7 Update README with new OS and arch support f731481 Merge pull request #199 from llucinat/wasm32-weakref f952e16 Fix weak reference macro redefinition in wasm32 target 97de1a4 Merge pull request #198 from gufe44/netbsd-fix-openlibm_weak_reference c4dca1e Add files via upload d4077aa Suggestions ea065f9 Update src/cdefs-compat.h 2080b23 NetBSD fix 14bf902 Merge pull request #195 from ode33/patch-1 3bb2215 Update README.md 33c8313 Create FUNDING.yml 72f33a3 wasm32 support (#192) f24b1bf Fix compilation of gcc when using openlibm as system libm (#190) 0f22aeb Create CNAME c68e7d2 Delete CNAME b524581 Rename doc -> docs ebbba43 Move website to doc/ on master instead of gh-pages branch 4e3d709 update ULPs for s390 (#187) 65d7406 Merge pull request #185 from sharkcz/s390x 2a1e568 s390(x) port cca41bc Merge branch 'master' of github.com:JuliaMath/openlibm 74b54c7 Add MIPS ce69bf1 Update references to JuliaLang to point to JuliaMath (#182) a96f074 Merge pull request mirage#130 from ginggs/enable-optimization-again c782ca2 Merge pull request #177 from JuliaMath/aa/windows 52df60b Update appveyor.yml ce33de1 Add Windows testing with AppVeyor 4971b56 Update Make.inc ca996ac Merge pull request #180 from JuliaLang/ginggs-0.5.6 73b3d88 Merge pull request #181 from CDLuminate/mipsport 3aa5c3b Merge pull request #174 from iniserve/master 4b4b41c Merge pull request #178 from JuliaLang/aa/upstream a4b3fde travis: Add mips, mipsel, mips64el build. ad9673e Makefile: clean mips/*.o 4dcc76e Using cdefs-compat.h and stdint.h instead <sys/types.h> fenv-softfloat.h file added SOFTFLOAT code parts are not tested. 4f5112e Support for mips architectures a24a5eb Enable optimization again for *int.c and *intf.c a40570b Bump version to 0.5.6 787652b msun: signed overflow in atan2 8d91ecb Add TOOLPREFIX git-subtree-dir: openlibm git-subtree-split: ae2d91698508701c83cab83714d42a1146dccf85
Summary: From Steve Kargl: Paul Zimmermann has identified a bug in Openlibm's powf(), which is identical to FreeBSD's libm. Both derived from fdlibm. JuliaMath/openlibm#212. Consider % cat h.c int main(void) { float x, y, z; x = 0x1.ffffecp-1F; y = -0x1.000002p+27F; z = 0x1.557a86p115F; printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z); return 0; } % cc -o h -fno-builtin h.c -lm && ./h 9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34 Reviewers: manu Subscribers: imp, andrew, emaste Differential Revision: https://reviews.freebsd.org/D31865 (cherry picked from commit 292815e)
@zimmermann6 has a publication comparing the accuracy of single precision functions across several libms.
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
The tests are with openlibm 0.4.1, and it would be nice to be able to try the latest release. In general musl fares better.
The text was updated successfully, but these errors were encountered: