1 /* -------------------------------------------------------------- */
2 /* (C)Copyright 2001,2008, */
3 /* International Business Machines Corporation, */
4 /* Sony Computer Entertainment, Incorporated, */
5 /* Toshiba Corporation, */
6 /* */
7 /* All Rights Reserved. */
8 /* */
9 /* Redistribution and use in source and binary forms, with or */
10 /* without modification, are permitted provided that the */
11 /* following conditions are met: */
12 /* */
13 /* - Redistributions of source code must retain the above copyright*/
14 /* notice, this list of conditions and the following disclaimer. */
15 /* */
16 /* - Redistributions in binary form must reproduce the above */
17 /* copyright notice, this list of conditions and the following */
18 /* disclaimer in the documentation and/or other materials */
19 /* provided with the distribution. */
20 /* */
21 /* - Neither the name of IBM Corporation nor the names of its */
22 /* contributors may be used to endorse or promote products */
23 /* derived from this software without specific prior written */
24 /* permission. */
25 /* */
26 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND */
27 /* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, */
28 /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
29 /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
30 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR */
31 /* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, */
32 /* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT */
33 /* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
34 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) */
35 /* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN */
36 /* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR */
37 /* OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, */
38 /* EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
39 /* -------------------------------------------------------------- */
40 /* PROLOG END TAG zYx */
41 #ifdef __SPU__
42
43 #ifndef _ACOSF4_H_
44 #define _ACOSF4_H_ 1
45
46 #include <spu_intrinsics.h>
47
48 #include "divf4.h"
49 #include "sqrtf4.h"
50
51 /*
52 * FUNCTION
53 * vector float _acosf4(vector float x)
54 *
55 * DESCRIPTION
56 * The _acosf4 function computes the arc cosine for a vector of values x;
57 * that is the values whose cosine is x. Results are undefined if x is
58 * outside the range [-1, 1].
59 *
60 * RETURNS
61 * The _acosf4 function returns the arc cosine in radians and the value is
62 * mathematically defined to be in the range [0, pi].
63 *
64 */
_acosf4(vector float x)65 static __inline vector float _acosf4(vector float x)
66 {
67 vec_float4 zero = spu_splats(0.0f);
68 vec_float4 half = spu_splats(0.5f);
69 vec_float4 one = spu_splats(1.0f);
70 vec_float4 two = spu_splats(2.0f);
71 vec_float4 pi = spu_splats(3.1415925026e+00f);
72 vec_float4 pio2_hi = spu_splats(1.5707962513e+00f);
73 vec_float4 pio2_lo = spu_splats(7.5497894159e-08f);
74
75 vec_float4 snan = (vec_float4)spu_splats((unsigned int)0x7FC00000);
76 vec_uint4 denorm_threshold = spu_splats((unsigned int)0x23000000);
77 vec_uint4 sign_mask = spu_splats((unsigned int)0x80000000);
78
79
80 vec_float4 p0 = (vec_float4)spu_splats((unsigned int)0x3E2AAAAB);
81 vec_float4 p1 = (vec_float4)spu_splats((unsigned int)0xBEA6B090);
82 vec_float4 p2 = (vec_float4)spu_splats((unsigned int)0x3E4E0AA8);
83 vec_float4 p3 = (vec_float4)spu_splats((unsigned int)0xBD241146);
84 vec_float4 p4 = (vec_float4)spu_splats((unsigned int)0x3A4F7F04);
85 vec_float4 p5 = (vec_float4)spu_splats((unsigned int)0x3811EF08);
86
87 vec_float4 q1 = (vec_float4)spu_splats((unsigned int)0xC019D139);
88 vec_float4 q2 = (vec_float4)spu_splats((unsigned int)0x4001572D);
89 vec_float4 q3 = (vec_float4)spu_splats((unsigned int)0xBF303361);
90 vec_float4 q4 = (vec_float4)spu_splats((unsigned int)0x3D9DC62E);
91
92
93 vec_uint4 x_abs = spu_andc((vec_uint4)x,sign_mask);
94 vec_uint4 x_pos = spu_cmpgt(sign_mask,(vec_uint4)x);
95
96
97 vec_uint4 almost_half = spu_splats((unsigned int)0x3EFFFFFF);
98 vec_uint4 sel0 = spu_nand(spu_splats((unsigned int)0xFFFFFFFF),spu_cmpgt(x_abs,almost_half));
99 vec_uint4 sel1 = spu_andc(x_pos,sel0); // pos
100
101 vec_float4 za = spu_sel(spu_sel(spu_add(one,x),spu_sub(one,x),sel1) ,x,sel0);
102 vec_float4 zb = spu_sel(half,x,sel0);
103
104 vec_float4 z = spu_mul(za,zb);
105
106 vec_float4 p;
107 p = spu_madd(z,p5,p4);
108 p = spu_madd(p,z,p3);
109 p = spu_madd(p,z,p2);
110 p = spu_madd(p,z,p1);
111 p = spu_madd(p,z,p0);
112 p = spu_mul(p,z);
113
114 vec_float4 q;
115 q = spu_madd(z,q4,q3);
116 q = spu_madd(q,z,q2);
117 q = spu_madd(q,z,q1);
118 q = spu_madd(q,z,one);
119
120 // Only used by secondaries
121 vec_float4 s = _sqrtf4(z);
122
123 vec_float4 r = _divf4(p,q);
124
125 vec_float4 w1 = spu_msub(r,s,pio2_lo);
126
127
128 vec_float4 df = (vec_float4)spu_and((vec_uint4)s,0xFFFFF000);
129 vec_float4 c = _divf4(spu_nmsub(df,df,z),spu_add(s,df));
130 vec_float4 w2 = spu_madd(r,s,c);
131
132 vec_float4 result0 = spu_sub(pio2_hi,spu_sub(x,spu_nmsub(x,r,pio2_lo)));
133
134
135 vec_float4 result1 = spu_mul(two,spu_add(df,w2));
136 vec_float4 result2 = spu_nmsub(two,spu_add(s,w1),pi);
137
138
139 vec_float4 result;
140
141 result = spu_sel(result2,result1,sel1);
142 result = spu_sel(result,result0,sel0);
143
144 // If |x|==1 then:
145 // if x == 1, return 0
146 // else return pi
147
148 vec_uint4 abs_one = spu_cmpeq(x_abs,(vec_uint4)one);
149 vec_uint4 out_of_bounds = spu_cmpgt(x_abs,(vec_uint4)one);
150 vec_uint4 underflow = spu_cmpgt(denorm_threshold,x_abs);
151
152
153
154 result = spu_sel(result,spu_sel(pi,zero,x_pos),abs_one);
155
156 // If 1 < |x| then return sNaN
157 result = spu_sel(result,snan,out_of_bounds);
158
159 // If |x| < 2**-57, then return pi/2 (OFF BY 1 ULP)
160 result = spu_sel(result,spu_add(pio2_hi,pio2_lo),underflow);
161
162 return result;
163 }
164
165 #endif /* _ACOSF4_H_ */
166 #endif /* __SPU__ */
167
168
169