Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/erfq.c @ 145:1830386684a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc-9.2.0
author anatofuz
date Thu, 13 Feb 2020 11:34:05 +0900
parents 04ced10e8804
children
comparison
equal deleted inserted replaced
131:84e7813d76e9 145:1830386684a0
25 but WITHOUT ANY WARRANTY; without even the implied warranty of 25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details. 27 Lesser General Public License for more details.
28 28
29 You should have received a copy of the GNU Lesser General Public 29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, write to the Free Software 30 License along with this library; if not, see
31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ 31 <http://www.gnu.org/licenses/>. */
32 32
33 /* __float128 erfq(__float128 x) 33 /* double erf(double x)
34 * __float128 erfcq(__float128 x) 34 * double erfc(double x)
35 * x 35 * x
36 * 2 |\ 36 * 2 |\
37 * erf(x) = --------- | exp(-t*t)dt 37 * erf(x) = --------- | exp(-t*t)dt
38 * sqrt(pi) \| 38 * sqrt(pi) \|
39 * 0 39 * 0
94 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, 94 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
95 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, 95 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
96 * erfc/erf(NaN) is NaN 96 * erfc/erf(NaN) is NaN
97 */ 97 */
98 98
99 #include <errno.h>
100 #include "quadmath-imp.h" 99 #include "quadmath-imp.h"
101
102
103
104 __float128 erfcq (__float128);
105
106 100
107 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ 101 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
108 102
109 static __float128 103 static __float128
110 neval (__float128 x, const __float128 *p, int n) 104 neval (__float128 x, const __float128 *p, int n)
141 135
142 136
143 137
144 static const __float128 138 static const __float128
145 tiny = 1e-4931Q, 139 tiny = 1e-4931Q,
146 one = 1.0Q, 140 one = 1,
147 two = 2.0Q, 141 two = 2,
148 /* 2/sqrt(pi) - 1 */ 142 /* 2/sqrt(pi) - 1 */
149 efx = 1.2837916709551257389615890312154517168810E-1Q; 143 efx = 1.2837916709551257389615890312154517168810E-1Q;
150 144
151 145
152 /* erf(x) = x + x R(x^2) 146 /* erf(x) = x + x R(x^2)
808 802
809 803
810 __float128 804 __float128
811 erfcq (__float128 x) 805 erfcq (__float128 x)
812 { 806 {
813 __float128 y = 0.0Q, z, p, r; 807 __float128 y, z, p, r;
814 int32_t i, ix, sign; 808 int32_t i, ix, sign;
815 ieee854_float128 u; 809 ieee854_float128 u;
816 810
817 u.value = x; 811 u.value = x;
818 sign = u.words32.w0; 812 sign = u.words32.w0;
866 z = x - 0.875Q; 860 z = x - 0.875Q;
867 y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18); 861 y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
868 y += C18a; 862 y += C18a;
869 break; 863 break;
870 case 8: 864 case 8:
871 z = x - 1.0Q; 865 z = x - 1;
872 y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); 866 y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
873 y += C19a; 867 y += C19a;
874 break; 868 break;
875 default: /* i == 9. */ 869 default: /* i == 9. */
876 z = x - 1.125Q; 870 z = x - 1.125Q;
877 y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); 871 y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
878 y += C20a; 872 y += C20a;
879 break; 873 break;
880 } 874 }
881 if (sign & 0x80000000) 875 if (sign & 0x80000000)
882 y = 2.0Q - y; 876 y = 2 - y;
883 return y; 877 return y;
884 } 878 }
885 /* 1.25 < |x| < 107 */ 879 /* 1.25 < |x| < 107 */
886 if (ix < 0x4005ac00) 880 if (ix < 0x4005ac00)
887 { 881 {
922 } 916 }
923 u.value = x; 917 u.value = x;
924 u.words32.w3 = 0; 918 u.words32.w3 = 0;
925 u.words32.w2 &= 0xfe000000; 919 u.words32.w2 &= 0xfe000000;
926 z = u.value; 920 z = u.value;
927 r = expq (-z * z - 0.5625) * expq ((z - x) * (z + x) + p); 921 r = expq (-z * z - 0.5625) *
922 expq ((z - x) * (z + x) + p);
928 if ((sign & 0x80000000) == 0) 923 if ((sign & 0x80000000) == 0)
929 { 924 {
930 __float128 ret = r / x; 925 __float128 ret = r / x;
931 if (ret == 0) 926 if (ret == 0)
932 errno = ERANGE; 927 errno = ERANGE;