Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/sincosq_kernel.c @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents 561a7518be6b
children 1830386684a0
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
1 /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>. 1 /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>.
2 Copyright (C) 1999 Free Software Foundation, Inc. 2 Copyright (C) 1999-2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library. 3 This file is part of the GNU C Library.
4 Contributed by Jakub Jelinek <jj@ultra.linux.cz> 4 Contributed by Jakub Jelinek <jj@ultra.linux.cz>
5 5
6 The GNU C Library is free software; you can redistribute it and/or 6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public 7 modify it under the terms of the GNU Lesser General Public
108 if (tix < 0x3ffc3000) /* |x| < 0.1484375 */ 108 if (tix < 0x3ffc3000) /* |x| < 0.1484375 */
109 { 109 {
110 /* Argument is small enough to approximate it by a Chebyshev 110 /* Argument is small enough to approximate it by a Chebyshev
111 polynomial of degree 16(17). */ 111 polynomial of degree 16(17). */
112 if (tix < 0x3fc60000) /* |x| < 2^-57 */ 112 if (tix < 0x3fc60000) /* |x| < 2^-57 */
113 if (!((int)x)) /* generate inexact */ 113 {
114 { 114 math_check_force_underflow (x);
115 *sinx = x; 115 if (!((int)x)) /* generate inexact */
116 *cosx = ONE; 116 {
117 return; 117 *sinx = x;
118 } 118 *cosx = ONE;
119 return;
120 }
121 }
119 z = x * x; 122 z = x * x;
120 *sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+ 123 *sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
121 z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8))))))))); 124 z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
122 *cosx = ONE + (z*(COS1+z*(COS2+z*(COS3+z*(COS4+ 125 *cosx = ONE + (z*(COS1+z*(COS2+z*(COS3+z*(COS4+
123 z*(COS5+z*(COS6+z*(COS7+z*COS8)))))))); 126 z*(COS5+z*(COS6+z*(COS7+z*COS8))))))));
124 } 127 }
125 else 128 else
126 { 129 {
127 /* So that we don't have to use too large polynomial, we find 130 /* So that we don't have to use too large polynomial, we find
128 l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83 131 l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83
129 possible values for h. We look up cosl(h) and sinl(h) in 132 possible values for h. We look up cosq(h) and sinq(h) in
130 pre-computed tables, compute cosl(l) and sinl(l) using a 133 pre-computed tables, compute cosq(l) and sinq(l) using a
131 Chebyshev polynomial of degree 10(11) and compute 134 Chebyshev polynomial of degree 10(11) and compute
132 sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and 135 sinq(h+l) = sinq(h)cosq(l) + cosq(h)sinq(l) and
133 cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */ 136 cosq(h+l) = cosq(h)cosq(l) - sinq(h)sinq(l). */
134 index = 0x3ffe - (tix >> 16); 137 index = 0x3ffe - (tix >> 16);
135 hix = (tix + (0x200 << index)) & (0xfffffc00 << index); 138 hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
136 x = fabsq (x); 139 if (signbitq (x))
140 {
141 x = -x;
142 y = -y;
143 }
137 switch (index) 144 switch (index)
138 { 145 {
139 case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break; 146 case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
140 case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break; 147 case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
141 default: 148 default: