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 _SINHD2_H_
40 #define _SINHD2_H_ 1
41
42 #include <spu_intrinsics.h>
43
44 #include "expd2.h"
45 #include "recipd2.h"
46
47 /*
48 * FUNCTION
49 * vector double _sinhd2(vector double angle)
50 *
51 * DESCRIPTION
52 * The _sinhd2 function computes the hyperbolic sine of a vector of
53 * angles (expressed in radians) to an accuracy of a double precision
54 * floating point.
55 */
_sinhd2(vector double x)56 static __inline vector double _sinhd2(vector double x)
57 {
58 // Coefficents for the power series
59 vec_double2 f03 = spu_splats(1.66666666666666666667E-01); // 1/(3!)
60 vec_double2 f05 = spu_splats(8.33333333333333333333E-03); // 1/(5!)
61 vec_double2 f07 = spu_splats(1.98412698412698412698E-04); // 1/(7!)
62 vec_double2 f09 = spu_splats(2.75573192239858906526E-06); // 1/(9!)
63 vec_double2 f11 = spu_splats(2.50521083854417187751E-08); // 1/(11!)
64 vec_double2 f13 = spu_splats(1.60590438368216145994E-10); // 1/(13!)
65 vec_double2 f15 = spu_splats(7.64716373181981647590E-13); // 1/(15!)
66 vec_double2 f17 = spu_splats(2.81145725434552076320E-15); // 1/(17!)
67 vec_double2 f19 = spu_splats(8.22063524662432971696E-18); // 1/(19!)
68 vec_double2 f21 = spu_splats(1.95729410633912612308E-20); // 1/(21!)
69 vec_double2 f23 = spu_splats(3.86817017063068403772E-23); // 1/(23!)
70
71
72
73 // Check if the input is within the range [ -1.0 ... 1.0 ]
74 // If it is, we want to use the power series, otherwise
75 // we want to use the 0.5 * (e^x - e^-x)
76
77 // round to float, check if within range. Results will be in
78 // slots 0 and 2, so we rotate right 4 bytes, and "or" with ourself
79 // to produce 64 bits of all 1's or 0's.
80 vec_uint4 use_exp = spu_cmpabsgt(spu_roundtf(x),spu_splats(1.0f));
81 use_exp = spu_or(use_exp,spu_rlmaskqwbyte(use_exp,-4));
82
83
84
85
86 // Perform the calculation of the power series using Horner's method
87 vec_double2 result;
88 vec_double2 x2 = spu_mul(x,x);
89 result = spu_madd(x2,f23,f21);
90 result = spu_madd(x2,result,f19);
91 result = spu_madd(x2,result,f17);
92 result = spu_madd(x2,result,f15);
93 result = spu_madd(x2,result,f13);
94 result = spu_madd(x2,result,f11);
95 result = spu_madd(x2,result,f09);
96 result = spu_madd(x2,result,f07);
97 result = spu_madd(x2,result,f05);
98 result = spu_madd(x2,result,f03);
99 result = spu_madd(x2,result,spu_splats(1.0));
100 result = spu_mul(x,result);
101
102
103 // Perform calculation as a function of 0.5 * (e^x - e^-x)
104 vec_double2 ex = _expd2(x);
105 vec_double2 ex_inv = _recipd2(ex);
106
107 vec_double2 r2= spu_sub(ex,ex_inv);
108 r2 = spu_mul(r2,spu_splats(0.5));
109
110
111 // Select either the power series or exp version
112 result = spu_sel(result,r2,(vec_ullong2)use_exp);
113
114
115
116 return result;
117
118 }
119
120 #endif /* _SINHD2_H_ */
121 #endif /* __SPU__ */
122
