From: Markus M. <mar...@gm...> - 2022-07-26 14:45:30
|
Hello, We noticed in downstream Octave that the results of complex inverse trigonometric functions are off (sometimes by pi/2) for large input on MinGW. See also this downstream report: https://savannah.gnu.org/bugs/index.php?49091 Afaict, this is caused by the following: The implementation of cacos uses cacosh, and casin uses casinh. "cacosh" is implemented using the identity: cacosh(z) = log(z + sqrt(z*z - 1)) "casinh" is implemented using the identity: casinh(z) = log(z + sqrt(z*z + 1)) IIUC, "z*z" might overflow or become numerically unstable for large floating point "z" in both cases. The proposed patches use approximations for large z that avoid squaring "z" for large values of "z". For large "z", "z + sqrt(z*z + 1)" or "z + sqrt(z*z - 1)" is approximately "2*z". Additionally, it uses symmetries to get the correct results for input in any quadrant. I compared the results for some large z (z=±1e150; z=±1e150i; z=±1e150±1e150i; for cacos, cacosh, casin, and casinh) to the results of Wolfram Alfa. E.g.: https://www.wolframalpha.com/input?i=acos%281e150%29 With the proposed approximations, they coincide to 16 significant decimal digits in the real and imaginary parts for double precision floating point. Without the approximations (the current base line), the deviations are about pi/2 for some (e.g., the real part of cacos(1e150+0i) or cacos(-1e150+0i)). This is the first time I'm writing to this list. So, please let me know if this is not the correct place or format to propose patches. I'm looking forward to any feedback and hope that something similar could be accepted. Markus |
From: Martin S. <ma...@ma...> - 2022-08-02 13:09:44
|
On Tue, 26 Jul 2022, Markus Mützel wrote: > Hello, > > We noticed in downstream Octave that the results of complex inverse trigonometric functions are off (sometimes by pi/2) for large input on MinGW. > See also this downstream report: > https://savannah.gnu.org/bugs/index.php?49091 > > Afaict, this is caused by the following: > The implementation of cacos uses cacosh, and casin uses casinh. > "cacosh" is implemented using the identity: cacosh(z) = log(z + sqrt(z*z - 1)) > "casinh" is implemented using the identity: casinh(z) = log(z + sqrt(z*z + 1)) > > IIUC, "z*z" might overflow or become numerically unstable for large floating point "z" in both cases. > > The proposed patches use approximations for large z that avoid squaring "z" for large values of "z". For large "z", "z + sqrt(z*z + 1)" or "z + sqrt(z*z - 1)" is approximately "2*z". Additionally, it uses symmetries to get the correct results for input in any quadrant. > > I compared the results for some large z (z=±1e150; z=±1e150i; z=±1e150±1e150i; for cacos, cacosh, casin, and casinh) to the results of Wolfram Alfa. > E.g.: https://www.wolframalpha.com/input?i=acos%281e150%29 > With the proposed approximations, they coincide to 16 significant decimal digits in the real and imaginary parts for double precision floating point. Without the approximations (the current base line), the deviations are about pi/2 for some (e.g., the real part of cacos(1e150+0i) or cacos(-1e150+0i)). > > This is the first time I'm writing to this list. So, please let me know if this is not the correct place or format to propose patches. This is the right place for such patches. The patches look good to me (although it took some time to work through the math - it's been a while since I used my complex math); I didn't try to check how the symmetries that are used work out here, but it seems reasonable and I presume you tested those aspects. Minor nitpick; in cacosh in the comment, "half plain" should be "half plane". I can push the patch later if there's no other disagreement about it. // Martin |
From: Markus M. <mar...@gm...> - 2022-08-02 15:13:56
Attachments:
0002-cacosh-Use-approximation-for-large-input.patch
|
Am 02. August 2022 um 15:09 Uhr schrieb "Martin Storsjö": > On Tue, 26 Jul 2022, Markus Mützel wrote: > > > Hello, > > > > We noticed in downstream Octave that the results of complex inverse trigonometric functions are off (sometimes by pi/2) for large input on MinGW. > > See also this downstream report: > > https://savannah.gnu.org/bugs/index.php?49091 > > > > Afaict, this is caused by the following: > > The implementation of cacos uses cacosh, and casin uses casinh. > > "cacosh" is implemented using the identity: cacosh(z) = log(z + sqrt(z*z - 1)) > > "casinh" is implemented using the identity: casinh(z) = log(z + sqrt(z*z + 1)) > > > > IIUC, "z*z" might overflow or become numerically unstable for large floating point "z" in both cases. > > > > The proposed patches use approximations for large z that avoid squaring "z" for large values of "z". For large "z", "z + sqrt(z*z + 1)" or "z + sqrt(z*z - 1)" is approximately "2*z". Additionally, it uses symmetries to get the correct results for input in any quadrant. > > > > I compared the results for some large z (z=±1e150; z=±1e150i; z=±1e150±1e150i; for cacos, cacosh, casin, and casinh) to the results of Wolfram Alfa. > > E.g.: https://www.wolframalpha.com/input?i=acos%281e150%29 > > With the proposed approximations, they coincide to 16 significant decimal digits in the real and imaginary parts for double precision floating point. Without the approximations (the current base line), the deviations are about pi/2 for some (e.g., the real part of cacos(1e150+0i) or cacos(-1e150+0i)). > > > > This is the first time I'm writing to this list. So, please let me know if this is not the correct place or format to propose patches. > > This is the right place for such patches. > > The patches look good to me (although it took some time to work through > the math - it's been a while since I used my complex math); I didn't try > to check how the symmetries that are used work out here, but it seems > reasonable and I presume you tested those aspects. Thank you for reviewing. The symmetries should be correct in theory. And I also tested them with a couple of different inputs. > Minor nitpick; in cacosh in the comment, "half plain" should be "half > plane". I'm attaching an updated patch with that spelling error corrected. > I can push the patch later if there's no other disagreement about it. Looking forward to that. Markus |
From: Martin S. <ma...@ma...> - 2022-08-02 22:01:45
|
On Tue, 2 Aug 2022, Markus Mützel wrote: > Am 02. August 2022 um 15:09 Uhr schrieb "Martin Storsjö": >> On Tue, 26 Jul 2022, Markus Mützel wrote: >> >>> Hello, >>> >>> We noticed in downstream Octave that the results of complex inverse trigonometric functions are off (sometimes by pi/2) for large input on MinGW. >>> See also this downstream report: >>> https://savannah.gnu.org/bugs/index.php?49091 >>> >>> Afaict, this is caused by the following: >>> The implementation of cacos uses cacosh, and casin uses casinh. >>> "cacosh" is implemented using the identity: cacosh(z) = log(z + sqrt(z*z - 1)) >>> "casinh" is implemented using the identity: casinh(z) = log(z + sqrt(z*z + 1)) >>> >>> IIUC, "z*z" might overflow or become numerically unstable for large floating point "z" in both cases. >>> >>> The proposed patches use approximations for large z that avoid squaring "z" for large values of "z". For large "z", "z + sqrt(z*z + 1)" or "z + sqrt(z*z - 1)" is approximately "2*z". Additionally, it uses symmetries to get the correct results for input in any quadrant. >>> >>> I compared the results for some large z (z=±1e150; z=±1e150i; z=±1e150±1e150i; for cacos, cacosh, casin, and casinh) to the results of Wolfram Alfa. >>> E.g.: https://www.wolframalpha.com/input?i=acos%281e150%29 >>> With the proposed approximations, they coincide to 16 significant decimal digits in the real and imaginary parts for double precision floating point. Without the approximations (the current base line), the deviations are about pi/2 for some (e.g., the real part of cacos(1e150+0i) or cacos(-1e150+0i)). >>> >>> This is the first time I'm writing to this list. So, please let me know if this is not the correct place or format to propose patches. >> >> This is the right place for such patches. >> >> The patches look good to me (although it took some time to work through >> the math - it's been a while since I used my complex math); I didn't try >> to check how the symmetries that are used work out here, but it seems >> reasonable and I presume you tested those aspects. > > Thank you for reviewing. > The symmetries should be correct in theory. And I also tested them with a couple of different inputs. > >> Minor nitpick; in cacosh in the comment, "half plain" should be "half >> plane". > > I'm attaching an updated patch with that spelling error corrected. > >> I can push the patch later if there's no other disagreement about it. > > Looking forward to that. Pushed now, thanks for the patch! // Martin |