1 /* -------------------------------------------------------------- */
2 /* (C)Copyright 2006,2008, */
3 /* International Business Machines Corporation */
4 /* All Rights Reserved. */
5 /* */
6 /* Redistribution and use in source and binary forms, with or */
7 /* without modification, are permitted provided that the */
8 /* following conditions are met: */
9 /* */
10 /* - Redistributions of source code must retain the above copyright*/
11 /* notice, this list of conditions and the following disclaimer. */
12 /* */
13 /* - Redistributions in binary form must reproduce the above */
14 /* copyright notice, this list of conditions and the following */
15 /* disclaimer in the documentation and/or other materials */
16 /* provided with the distribution. */
17 /* */
18 /* - Neither the name of IBM Corporation nor the names of its */
19 /* contributors may be used to endorse or promote products */
20 /* derived from this software without specific prior written */
21 /* permission. */
22 /* */
23 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND */
24 /* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, */
25 /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
26 /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
27 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR */
28 /* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, */
29 /* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT */
30 /* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
31 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) */
32 /* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN */
33 /* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR */
34 /* OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, */
35 /* EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
36 /* -------------------------------------------------------------- */
37 /* PROLOG END TAG zYx */
38 #ifdef __SPU__
39 #ifndef _COSHD2_H_
40 #define _COSHD2_H_ 1
41
42 #include <spu_intrinsics.h>
43
44 #include "expd2.h"
45 #include "recipd2.h"
46
47 /*
48 * FUNCTION
49 * vector float _coshd2(vector double angle)
50 *
51 * DESCRIPTION
52 * _coshd2 computes the hyperbolic cosines of a vector of angles
53 * (expressed in radians) to an accuracy of a double precision
54 * floating point.
55 */
_coshd2(vector double x)56 static __inline vector double _coshd2(vector double x)
57 {
58
59 // Coefficents for the power series
60 vec_double2 f02 = spu_splats(5.00000000000000000000E-1); // 1/(2!)
61 vec_double2 f04 = spu_splats(4.16666666666666666667E-2); // 1/(4!)
62 vec_double2 f06 = spu_splats(1.38888888888888888889E-3); // 1/(6!)
63 vec_double2 f08 = spu_splats(2.48015873015873015873E-5); // 1/(8!)
64 vec_double2 f10 = spu_splats(2.75573192239858906526E-7); // 1/(10!)
65 vec_double2 f12 = spu_splats(2.08767569878680989792E-9); // 1/(12!)
66 vec_double2 f14 = spu_splats(1.14707455977297247139E-11); // 1/(14!)
67 vec_double2 f16 = spu_splats(4.77947733238738529744E-14); // 1/(16!)
68 vec_double2 f18 = spu_splats(1.56192069685862264622E-16); // 1/(18!)
69 vec_double2 f20 = spu_splats(4.11031762331216485848E-19); // 1/(20!)
70 vec_double2 f22 = spu_splats(8.89679139245057328675E-22); // 1/(22!)
71
72 // Check if the input is within the range [ -1.0 ... 1.0 ]
73 // If it is, we want to use the power series, otherwise
74 // we want to use the 0.5 * (e^x + e^-x)
75
76 // round to float, check if within range. Results will be in
77 // slots 0 and 2, so we rotate right 4 bytes, and "or" with ourself
78 // to produce 64 bits of all 1's or 0's.
79 vec_uint4 use_exp = spu_cmpabsgt(spu_roundtf(x),spu_splats(1.0f));
80 use_exp = spu_or(use_exp,spu_rlmaskqwbyte(use_exp,-4));
81
82
83 // Perform the calculation of the power series using Horner's method
84 vec_double2 result;
85 vec_double2 x2 = spu_mul(x,x);
86 result = spu_madd(x2,f22,f20);
87 result = spu_madd(x2,result,f18);
88 result = spu_madd(x2,result,f16);
89 result = spu_madd(x2,result,f14);
90 result = spu_madd(x2,result,f12);
91 result = spu_madd(x2,result,f10);
92 result = spu_madd(x2,result,f08);
93 result = spu_madd(x2,result,f06);
94 result = spu_madd(x2,result,f04);
95 result = spu_madd(x2,result,f02);
96 result = spu_madd(x2,result,spu_splats(1.0));
97
98
99 // Perform calculation as a function of 0.5 * (e^x + e^-x)
100 vec_double2 ex = _expd2(x);
101 vec_double2 ex_inv = _recipd2(ex);
102
103 vec_double2 r2= spu_add(ex,ex_inv);
104 r2 = spu_mul(r2,f02); // we can reuse f02 here
105
106
107 // Select either the power series or exp version
108 result = spu_sel(result,r2,(vec_ullong2)use_exp);
109
110 return result;
111
112 }
113
114 #endif /* _COSHD2_H_ */
115 #endif /* __SPU__ */
116
