Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/j1q.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
4 * 4 *
5 * 5 *
6 * 6 *
7 * SYNOPSIS: 7 * SYNOPSIS:
8 * 8 *
9 * __float128 x, y, j1q(); 9 * long double x, y, j1l();
10 * 10 *
11 * y = j1q( x ); 11 * y = j1l( x );
12 * 12 *
13 * 13 *
14 * 14 *
15 * DESCRIPTION: 15 * DESCRIPTION:
16 * 16 *
50 * 50 *
51 * 51 *
52 * 52 *
53 * SYNOPSIS: 53 * SYNOPSIS:
54 * 54 *
55 * __float128, y, y1q(); 55 * double x, y, y1l();
56 * 56 *
57 * y = y1q( x ); 57 * y = y1l( x );
58 * 58 *
59 * 59 *
60 * 60 *
61 * DESCRIPTION: 61 * DESCRIPTION:
62 * 62 *
90 but WITHOUT ANY WARRANTY; without even the implied warranty of 90 but WITHOUT ANY WARRANTY; without even the implied warranty of
91 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 91 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
92 Lesser General Public License for more details. 92 Lesser General Public License for more details.
93 93
94 You should have received a copy of the GNU Lesser General Public 94 You should have received a copy of the GNU Lesser General Public
95 License along with this library; if not, write to the Free Software 95 License along with this library; if not, see
96 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ 96 <http://www.gnu.org/licenses/>. */
97 97
98 #include <errno.h>
99 #include "quadmath-imp.h" 98 #include "quadmath-imp.h"
100 99
101 /* 1 / sqrt(pi) */ 100 /* 1 / sqrt(pi) */
102 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; 101 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
103 /* 2 / pi */ 102 /* 2 / pi */
104 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; 103 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
105 static const __float128 zero = 0.0Q; 104 static const __float128 zero = 0;
106 105
107 /* J1(x) = .5x + x x^2 R(x^2) 106 /* J1(x) = .5x + x x^2 R(x^2)
108 Peak relative error 1.9e-35 107 Peak relative error 1.9e-35
109 0 <= x <= 2 */ 108 0 <= x <= 2 */
110 #define NJ0_2N 6 109 #define NJ0_2N 6
124 5.934143516050192600795972192791775226920E13Q, 123 5.934143516050192600795972192791775226920E13Q,
125 2.168000911950620999091479265214368352883E11Q, 124 2.168000911950620999091479265214368352883E11Q,
126 5.673775894803172808323058205986256928794E8Q, 125 5.673775894803172808323058205986256928794E8Q,
127 1.080329960080981204840966206372671147224E6Q, 126 1.080329960080981204840966206372671147224E6Q,
128 1.411951256636576283942477881535283304912E3Q, 127 1.411951256636576283942477881535283304912E3Q,
129 /* 1.000000000000000000000000000000000000000E0Q */ 128 /* 1.000000000000000000000000000000000000000E0L */
130 }; 129 };
131 130
132 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), 131 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
133 0 <= 1/x <= .0625 132 0 <= 1/x <= .0625
134 Peak relative error 3.6e-36 */ 133 Peak relative error 3.6e-36 */
688 if (! finiteq (x)) 687 if (! finiteq (x))
689 { 688 {
690 if (x != x) 689 if (x != x)
691 return x + x; 690 return x + x;
692 else 691 else
693 return 0.0Q; 692 return 0;
694 } 693 }
695 if (x == 0.0Q) 694 if (x == 0)
696 return x; 695 return x;
697 xx = fabsq (x); 696 xx = fabsq (x);
698 if (xx <= 0x1p-58Q) 697 if (xx <= 0x1p-58Q)
699 { 698 {
700 __float128 ret = x * 0.5Q; 699 __float128 ret = x * 0.5Q;
701 math_check_force_underflow (ret); 700 math_check_force_underflow (ret);
702 if (ret == 0) 701 if (ret == 0)
703 errno = ERANGE; 702 errno = ERANGE;
704 return ret; 703 return ret;
705 } 704 }
706 if (xx <= 2.0Q) 705 if (xx <= 2)
707 { 706 {
708 /* 0 <= x <= 2 */ 707 /* 0 <= x <= 2 */
709 z = xx * xx; 708 z = xx * xx;
710 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); 709 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
711 p += 0.5Q * xx; 710 p += 0.5Q * xx;
721 = -1/sqrt(2) * (sin(x) + cos(x)) 720 = -1/sqrt(2) * (sin(x) + cos(x))
722 cf. Fdlibm. */ 721 cf. Fdlibm. */
723 sincosq (xx, &s, &c); 722 sincosq (xx, &s, &c);
724 ss = -s - c; 723 ss = -s - c;
725 cc = s - c; 724 cc = s - c;
726 if (xx <= FLT128_MAX / 2.0Q) 725 if (xx <= FLT128_MAX / 2)
727 { 726 {
728 z = cosq (xx + xx); 727 z = cosq (xx + xx);
729 if ((s * c) > 0) 728 if ((s * c) > 0)
730 cc = z / ss; 729 cc = z / ss;
731 else 730 else
738 if (x < 0) 737 if (x < 0)
739 z = -z; 738 z = -z;
740 return z; 739 return z;
741 } 740 }
742 741
743 xinv = 1.0Q / xx; 742 xinv = 1 / xx;
744 z = xinv * xinv; 743 z = xinv * xinv;
745 if (xinv <= 0.25) 744 if (xinv <= 0.25)
746 { 745 {
747 if (xinv <= 0.125) 746 if (xinv <= 0.125)
748 { 747 {
796 { 795 {
797 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); 796 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
798 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); 797 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
799 } 798 }
800 } 799 }
801 p = 1.0Q + z * p; 800 p = 1 + z * p;
802 q = z * q; 801 q = z * q;
803 q = q * xinv + 0.375Q * xinv; 802 q = q * xinv + 0.375Q * xinv;
804 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); 803 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
805 if (x < 0) 804 if (x < 0)
806 z = -z; 805 z = -z;
807 return z; 806 return z;
808 } 807 }
809 808
810 809
810
811 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) 811 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
812 Peak relative error 6.2e-38 812 Peak relative error 6.2e-38
813 0 <= x <= 2 */ 813 0 <= x <= 2 */
814 #define NY0_2N 7 814 #define NY0_2N 7
815 static __float128 Y0_2N[NY0_2N + 1] = { 815 static const __float128 Y0_2N[NY0_2N + 1] = {
816 -6.804415404830253804408698161694720833249E19Q, 816 -6.804415404830253804408698161694720833249E19Q,
817 1.805450517967019908027153056150465849237E19Q, 817 1.805450517967019908027153056150465849237E19Q,
818 -8.065747497063694098810419456383006737312E17Q, 818 -8.065747497063694098810419456383006737312E17Q,
819 1.401336667383028259295830955439028236299E16Q, 819 1.401336667383028259295830955439028236299E16Q,
820 -1.171654432898137585000399489686629680230E14Q, 820 -1.171654432898137585000399489686629680230E14Q,
821 5.061267920943853732895341125243428129150E11Q, 821 5.061267920943853732895341125243428129150E11Q,
822 -1.096677850566094204586208610960870217970E9Q, 822 -1.096677850566094204586208610960870217970E9Q,
823 9.541172044989995856117187515882879304461E5Q, 823 9.541172044989995856117187515882879304461E5Q,
824 }; 824 };
825 #define NY0_2D 7 825 #define NY0_2D 7
826 static __float128 Y0_2D[NY0_2D + 1] = { 826 static const __float128 Y0_2D[NY0_2D + 1] = {
827 3.470629591820267059538637461549677594549E20Q, 827 3.470629591820267059538637461549677594549E20Q,
828 4.120796439009916326855848107545425217219E18Q, 828 4.120796439009916326855848107545425217219E18Q,
829 2.477653371652018249749350657387030814542E16Q, 829 2.477653371652018249749350657387030814542E16Q,
830 9.954678543353888958177169349272167762797E13Q, 830 9.954678543353888958177169349272167762797E13Q,
831 2.957927997613630118216218290262851197754E11Q, 831 2.957927997613630118216218290262851197754E11Q,
843 { 843 {
844 __float128 xx, xinv, z, p, q, c, s, cc, ss; 844 __float128 xx, xinv, z, p, q, c, s, cc, ss;
845 845
846 if (! finiteq (x)) 846 if (! finiteq (x))
847 return 1 / (x + x * x); 847 return 1 / (x + x * x);
848 if (x <= 0.0Q) 848 if (x <= 0)
849 { 849 {
850 if (x < 0.0Q) 850 if (x < 0)
851 return (zero / (zero * x)); 851 return (zero / (zero * x));
852 return -1 / zero; /* -inf and divide by zero exception. */ 852 return -1 / zero; /* -inf and divide by zero exception. */
853 } 853 }
854 xx = fabsq (x); 854 xx = fabsq (x);
855 if (xx <= 0x1p-114) 855 if (xx <= 0x1p-114)
857 z = -TWOOPI / x; 857 z = -TWOOPI / x;
858 if (isinfq (z)) 858 if (isinfq (z))
859 errno = ERANGE; 859 errno = ERANGE;
860 return z; 860 return z;
861 } 861 }
862 if (xx <= 2.0Q) 862 if (xx <= 2)
863 { 863 {
864 /* 0 <= x <= 2 */ 864 /* 0 <= x <= 2 */
865 /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */ 865 SET_RESTORE_ROUNDF128 (FE_TONEAREST);
866 z = xx * xx; 866 z = xx * xx;
867 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); 867 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
868 p = -TWOOPI / xx + p; 868 p = -TWOOPI / xx + p;
869 p = TWOOPI * logq (x) * j1q (x) + p; 869 p = TWOOPI * logq (x) * j1q (x) + p;
870 return p; 870 return p;
877 = -1/sqrt(2) * (sin(x) + cos(x)) 877 = -1/sqrt(2) * (sin(x) + cos(x))
878 cf. Fdlibm. */ 878 cf. Fdlibm. */
879 sincosq (xx, &s, &c); 879 sincosq (xx, &s, &c);
880 ss = -s - c; 880 ss = -s - c;
881 cc = s - c; 881 cc = s - c;
882 if (xx <= FLT128_MAX / 2.0Q) 882 if (xx <= FLT128_MAX / 2)
883 { 883 {
884 z = cosq (xx + xx); 884 z = cosq (xx + xx);
885 if ((s * c) > 0) 885 if ((s * c) > 0)
886 cc = z / ss; 886 cc = z / ss;
887 else 887 else
889 } 889 }
890 890
891 if (xx > 0x1p256Q) 891 if (xx > 0x1p256Q)
892 return ONEOSQPI * ss / sqrtq (xx); 892 return ONEOSQPI * ss / sqrtq (xx);
893 893
894 xinv = 1.0Q / xx; 894 xinv = 1 / xx;
895 z = xinv * xinv; 895 z = xinv * xinv;
896 if (xinv <= 0.25) 896 if (xinv <= 0.25)
897 { 897 {
898 if (xinv <= 0.125) 898 if (xinv <= 0.125)
899 { 899 {
947 { 947 {
948 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); 948 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
949 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); 949 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
950 } 950 }
951 } 951 }
952 p = 1.0Q + z * p; 952 p = 1 + z * p;
953 q = z * q; 953 q = z * q;
954 q = q * xinv + 0.375Q * xinv; 954 q = q * xinv + 0.375Q * xinv;
955 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx); 955 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
956 return z; 956 return z;
957 } 957 }