From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 74431 invoked by alias); 27 Jun 2017 19:22:32 -0000 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org Received: (qmail 74400 invoked by uid 89); 27 Jun 2017 19:22:30 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-1.8 required=5.0 tests=AWL,BAYES_00,RP_MATCHES_RCVD,SPF_PASS autolearn=ham version=3.3.2 spammy=H*Ad:D*lanl.gov, gsl-discuss, gsldiscuss, twenty X-HELO: proofpoint8.lanl.gov Received: from proofpoint8.lanl.gov (HELO proofpoint8.lanl.gov) (204.121.3.47) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with ESMTP; Tue, 27 Jun 2017 19:22:29 +0000 Received: from pps.filterd (proofpoint8.lanl.gov [127.0.0.1]) by proofpoint8.lanl.gov (8.16.0.20/8.16.0.20) with SMTP id v5RJMNOk113850; Tue, 27 Jun 2017 13:22:27 -0600 Received: from mailrelay2.lanl.gov (mailrelay2.lanl.gov [128.165.4.103]) by proofpoint8.lanl.gov with ESMTP id 2bbcvm3ybu-1; Tue, 27 Jun 2017 13:22:27 -0600 Received: from localhost (localhost.localdomain [127.0.0.1]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 02FAEE8B92E; Tue, 27 Jun 2017 13:22:27 -0600 (MDT) X-NIE-2-Virus-Scanner: amavisd-new at mailrelay2.lanl.gov Received: from nostromo.lanl.gov (nostromo.lanl.gov [130.55.124.72]) by mailrelay2.lanl.gov (Postfix) with ESMTP id BAFA0E8B926; Tue, 27 Jun 2017 13:22:26 -0600 (MDT) To: gsl-discuss@sourceware.org, alken@colorado.edu From: Gerard Jungman Subject: sin_pi and cos_pi Message-ID: Date: Tue, 27 Jun 2017 19:22:00 -0000 User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.0 MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="------------BFD5AD6A2D28B183FEB96510" X-Proofpoint-Virus-Version: vendor=fsecure engine=2.50.10432:,, definitions=2017-06-27_11:,, signatures=0 X-Proofpoint-Spam-Details: rule=notspam policy=default score=0 spamscore=0 suspectscore=0 malwarescore=0 phishscore=0 adultscore=0 bulkscore=0 classifier=spam adjust=0 reason=mlx scancount=1 engine=8.0.1-1703280000 definitions=main-1706270308 X-SW-Source: 2017-q2/txt/msg00000.txt.bz2 This is a multi-part message in MIME format. --------------BFD5AD6A2D28B183FEB96510 Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit Content-length: 3392 > After Gerards suggestion I've implemented (see attached) the routines > sin_pi and cos_pi. Essentially, modf and fmod do the trick. I'm a bit > surprised at how simple the solution is, so please check if you see > any problems, but the numerical tests all pass. > > Best, > Konrad Since this is really a design and implementation discussion, and not a bug discussion, I wanted to move it to gsl-discuss. Also, I have no idea how to interact with the bug forum... Anyway, I do have some comments about the method. modf() is almost certainly the right choice for the basic argument reduction. But the combination of round() and fmod() in the "sign" step is probably not ideal. Unfortunately, there are something like a dozen different functions for rounding in the standard library, and some of them have surprising consequences. It's a bit of a minefield. round() is one of the nasty surprise functions in the set. The specification requires it to round away from zero, regardless of the current rounding direction. That sounds like a nice and easy-to-understand choice, until you learn that the way it does this is by _always_ resetting the rounding direction before pushing the instructions into the pipeline. Typically (and certainly on x86-ish cpus) this requires flushing the pipeline before and after the operation. Very expensive. I was looking at the C99 standard rounding functions a few months ago, and some tests showed that round() can easily run twenty times slower than other choices in the set. Some of the functions are specified to not raise an FE_ type hardware exception. Again, this sounds ok, until you learn that this is done by disabling floating-point exceptions before executing the instructions. This operation itself is expensive, combined with the fact that it also has to flush the pipeline. Basically, any of the functions which require changing the state of the floating-point unit will probably have terrible performance implications. Amongst all the choices, testing showed that rint() was relatively harmless. modf() is also ok, as far as I can tell. Attached is an implementation that I was working on earlier this year. It is probably far from perfect. But maybe it has better performance. I would be happy if anybody could prove me wrong. But you'll need to get out your measurement tools and read the docs for all those dopey rounding functions... In any case, any definitive statement about the correct choices here would be a real contribution. I am still unsure about most of it. Finally, if you look deep inside the glibc implementation, you will find an implementation of sin_pi. It is hidden in the gamma function implementation where it is used in the reflection formula. It is not exported, since it is not a standard function. These functions were written at some Sun subsidiary in the early 1990s. They are probably very good, in context. But they are also written in this kind of internal library implementation gibberish style, which is impossible to understand and impossible to export for use elsewhere. You can find one in the file e_lgamma_r.c, for example. There is more than one, because of the support for different platforms. Still, it might be worthwhile to evaluate this. With some work, it might be possible to lift it out of its strange little glibc ghetto and put it to use elsewhere. -- G. Jungman --------------BFD5AD6A2D28B183FEB96510 Content-Type: application/gzip; name="trig-pi.tar.gz" Content-Transfer-Encoding: base64 Content-Disposition: attachment; filename="trig-pi.tar.gz" Content-length: 2949 H4sIADaqUlkAA+1a/2/iOhLvz/wVc13dE/QoEChUr+29E9uvrFpaUXqrlU5C JnHAT8FhE6eF3tv//WbshC+Fdvu023b3nkeUJM7MeDzz8YztoiIx2B6L8sYL UgVpd7dOV2e3Xlm8ZrTh1JzqTr26s1urblScas1pbED9JY3KKIkViwA2fk/k YMTko3xfe/+Tkkrjr68RZ0HJ/e59UIAbjZ1H4u84Tq02i3+l3sD479R2MP6V 727JGvqLx7+8BW44nmL0hwr2IH9YgCpGB/pTOOURizz4YAaegy38QHcoYhhH 4SBiI8BbP+Ic4tBXdyzi+zANE3CZhIh7IkZM9RPFQShg0iuHEYxCT/hT0oNt ifR4BGrIQfFoFEPo64fT9g12LbHzAK6SfiBcOBculzEHhl1TSzzkHlpIekji hGy4Tm2AkxAVMyVCuQ9c4PsIbnkU4zPUsj5ShUUII1KSZ4osjyAck1wBzZ1C wNRctPTI8Oej9EBIrXsYjnFEQ1SJY7wTQQB9DknM/SQokgpkho+t7tnlTRea 7U/wsdnpNNvdT/vIrIYhvuW33KgSo3EgUDOOK2JSTdF80nBx3Dk8Q5Hm+9Z5 q/sJBwEnrW77+PoaTi470ISrZqfbOrw5b3bg6qZzdXl9XIJrTlZxkn/Cw74O EnrR44qJIM7G/QnjGqNxgQdDdssxvi4Xt2ga0/D5euxICQtCOdCjROa5H0vQ 8kGGqggx2ngwVGq8Vy7f3d2VBjIphdGgHBgdcfk3sqeceyekGyQeh01KW6Xh 5rzlwA9CpkrD33LY5sNZ89/HvetW+/DyOvfO476QHHpoZe/68qZzeLwgN2Jq SGLveIDWrmuXiN1cLoczponw8EIcH8cxcFdocMlk1EfA9MVgoGGtp4wZKH5Y okJUJVwWBFPElwmykIojdzqq1L7ux8v3rVPI500Xhbxzfg4HB3D0/rx3gWHv HbVOCwW0xLzP/U4m9vxEur1YyN5Y9CiPp9IwKeT+m4PMXLmP924oY5W13MM/ aVr6+UkRfpGF+XsdLIDP+D4PPuvHeZwYB5l1/wKJ6aIChRWFsRb4DH+HKnJt O6UKMtL3KquHrJVSHbaN/vtVjgQ5SHaBA1mEn0efS7QZ/QBA4wOEpEoiCe1m m7R8wT+KJPE6jVJlS8t7NALUN5fCaB6yMUrq2YGdRxhQhVMWYS5wMtCEuAfJ sUZsk6n4+A+8lrQwzQ19gZvYyMdD4etc4HGphJruZQwAaHB+LPD2voCDigcS RwOUfWPTnDee+OP+Dz0okirr65JDpih70btqoYC3v/qa4evp1nTNG3fmSexn i23lzcNTN+VfK+iqcr1B3zV971Txe3/R3zGYChILHJFbBBPEJfejzFbyLY4n e8jxeP1+jp81p8N9vuOTP+V4pd88w91zvzuOdvYuObu8o7+ruqXxpPfVQ+8/ nBtxOux0HBnzl7W5BIFpc8lqLqn+CLlkljQMpGdIftkU8uJIVmsz9xumjmU/ b89z9cunjJfI1dsrTl6XIWiUazLEbSi8B2sNkyJm2aGY2Y9ZZrz45I7/Mplj C0eO3GnmwGd38fnHWJUIHn/buuTNMlD0vSfGK+S3FBFLlTrKKvUMID9M+nsW OJ5aO71d2vyJ0fHIQm4GD525o5XM/XCLCyZkixm8aHIxpeCFHrcwmAsdmMd3 qd4lzK6sF5cg+6BWpJtkUzLe+mjL0jNo8fy3NHyZPp4+/63Uq435+X+t5mxU nNpOrWrPf1+D7PmvPf/96c5/EbTXASE2mAICjZxPhvPJOIwUQjIcgd6olGmj gg4XaLZCkJTg9MMedBNOoIbqLjjVvZ1f93Z24eKoa3CP2qmqIjR9+HDR7J71 Tm7ah71up3XaO+thlT2/OTo+mp3VPsGS04pIT6/njoMkpr8cnyCqJGweblL5 np8p45ouHI1ptszWVpNC5n0DO/wYGGVIcxlG7A4BzvWB8gTB4iLmxkLL0KFz EHHmTemQOsbFHckx8tQoCViR5qRRgxiTSst4PBbIa9QxXFQSMwWXRYNkRCpw YieuWpgRcBWF2IRCRgLFk0DFOiDp8bbuacgCfztt0HKZyhixmfg+zjMEAe73 Yp0JmOsmEXOnWk+AvPyBUNY9Ymeix7p8pF7MdrS4L6Kl7zJn6HkrjNsp4xJv m7VnDHi/rGgcxrhextmAm0AhBc3PCNfdA7bUqCe9drRxDUnyz4m4ZQE5VIX6 XSD6mG2nOvqTgpG5Gwp3mE27gZ5Wita/XjhiGEoeRaHxrTFQy6RGlp/7z4H9 B+ibrdUt+v4M+l4dJxSnF8LJ2oPfhzhZyFJa9QJu0NQBp7o789LNY8Dph1gL +izGOkELUMCxJjxOfc3RR9yjGiO9BxDTXp3DKq1rGbhK0DGRMOq0/Gx3HHHS gg5K3R2iMWn1EpiNb4WXYPEiRxDOsMJzFwFAVvTDWz4PvI4Pg3jEqMIntJzB 7lUite+lDKdMuulCICuNafBSsFMBibPtWsHUxWjJ8NRhm4lErQkV1TDCNdOm gddHKtdBEN6ZTnCzeks7+FAWSVBzEmoYHdzh29GYFkYCK3DmLhb1hUrRJP0k W+eUv/Wob3996fsC5TJgxYN5DZxVQHOlo42/PV5SybS3Xi//v1G2/6NM8QI/ /dH0lf1fpe4s/P6ravZ/lYrd/70G4Yw7wg3LFCtkhOnIVE5KFiLQqYKy8jwZ mhU29yh9PLHMXv6RxkGsPBHq31Cs/G4jl8OiDAS+vPkPwdKR20Qf6FXra0/3 H1nZTFaP791l7sX6NjFH+eMIzfDzm/TvgkbJaXCY32U3/5GbRUp8cRFcI5Um 68q+PeyyZMmSJUuWLFmyZMmSJUuWLFmyZMmSJUuWLFmyZMmSJUuWLFmyZMmS JUuvTP8DAK25BwBQAAA= --------------BFD5AD6A2D28B183FEB96510--