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 _EXP2D2_H_
44 #define _EXP2D2_H_ 1
45
46 #include <spu_intrinsics.h>
47
48
49 /*
50 * FUNCTION
51 * vector double _exp2d2(vector double x)
52 *
53 * DESCRIPTION
54 * _exp2d2 computes 2 raised to the input x for each
55 * of the double word elements of x. Computation is
56 * performed by observing the 2^(a+b) = 2^a * 2^b.
57 * We decompose x into a and b (above) by letting.
58 * a = ceil(x), b = x - a;
59 *
60 * 2^a is easily computed by placing a into the exponent
61 * or a floating point number whose mantissa is all zeros.
62 *
63 * 2^b is computed using the polynomial approximation.
64 *
65 * __13_
66 * \
67 * \
68 * 2^x = / Ci*x^i
69 * /____
70 * i=0
71 *
72 * for x in the range 0.0 to 1.0.
73 *
74 */
75 #define EXP_C00 1.0
76 #define EXP_C01 6.93147180559945286227e-01
77 #define EXP_C02 2.40226506959100694072e-01
78 #define EXP_C03 5.55041086648215761801e-02
79 #define EXP_C04 9.61812910762847687873e-03
80 #define EXP_C05 1.33335581464284411157e-03
81 #define EXP_C06 1.54035303933816060656e-04
82 #define EXP_C07 1.52527338040598376946e-05
83 #define EXP_C08 1.32154867901443052734e-06
84 #define EXP_C09 1.01780860092396959520e-07
85 #define EXP_C10 7.05491162080112087744e-09
86 #define EXP_C11 4.44553827187081007394e-10
87 #define EXP_C12 2.56784359934881958182e-11
88 #define EXP_C13 1.36914888539041240648e-12
89
_exp2d2(vector double vx)90 static __inline vector double _exp2d2(vector double vx)
91 {
92 vec_int4 ix, exp;
93 vec_uint4 overflow, underflow;
94 vec_float4 vxf;
95 vec_double2 p1, p2, x2, x4, x8;
96 vec_double2 vy, vxw, out_of_range;
97
98 /* Compute: vxw = x - ceil(x)
99 */
100 vxw = spu_add(vx, spu_splats(0.5));
101 vxf = spu_roundtf(vxw);
102 ix = spu_convts(vxf, 0);
103 ix = spu_add(ix, (vec_int4)spu_andc(spu_cmpgt(spu_splats(0.0f), vxf), spu_cmpeq(ix, spu_splats((int)0x80000000))));
104 vxf = spu_convtf(ix, 0);
105 vxw = spu_sub(vx, spu_extend(vxf));
106
107 /* Detect overflow and underflow. If overflow, force the result
108 * to infinity (at the end).
109 */
110 exp = spu_shuffle(ix, ix, ((vec_uchar16) { 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 }));
111
112 overflow = spu_cmpgt(exp, 1023);
113 underflow = spu_cmpgt(exp, -1023);
114 out_of_range = (vec_double2)spu_and(overflow, ((vec_uint4) { 0x7FF00000, 0, 0x7FF00000, 0 }));
115
116 /* Calculate the result by evaluating the 13th order polynomial.
117 * For efficiency, the polynomial is broken into two parts and
118 * evaluate then using nested
119 *
120 * result = (((((c13*x + c12)*x + c11)*x + c10)*x + c9)*x + c8)*x^8 +
121 * ((((((c7*x + c6)*x + c5)*x + c4)*x + c3)*x + c2)*x + c1)*x + c0
122 */
123 p2 = spu_madd(spu_splats(EXP_C07), vxw, spu_splats(EXP_C06));
124 p1 = spu_madd(spu_splats(EXP_C13), vxw, spu_splats(EXP_C12));
125 x2 = spu_mul(vxw, vxw);
126 p2 = spu_madd(vxw, p2, spu_splats(EXP_C05));
127 p1 = spu_madd(vxw, p1, spu_splats(EXP_C11));
128 x4 = spu_mul(x2, x2);
129 p2 = spu_madd(vxw, p2, spu_splats(EXP_C04));
130 p1 = spu_madd(vxw, p1, spu_splats(EXP_C10));
131 p2 = spu_madd(vxw, p2, spu_splats(EXP_C03));
132 p1 = spu_madd(vxw, p1, spu_splats(EXP_C09));
133 x8 = spu_mul(x4, x4);
134 p2 = spu_madd(vxw, p2, spu_splats(EXP_C02));
135 p1 = spu_madd(vxw, p1, spu_splats(EXP_C08));
136 p2 = spu_madd(vxw, p2, spu_splats(EXP_C01));
137 p2 = spu_madd(vxw, p2, spu_splats(EXP_C00));
138 vy = spu_madd(x8, p1, p2);
139
140 /* Align the integer integer portion of x with the exponent.
141 */
142 ix = spu_sl(ix, ((vec_uint4) { 20, 32, 20, 32 }));
143 vy = (vec_double2)spu_add((vec_int4)vy, ix);
144
145 /* Select the result if not overflow or underflow. Otherwise select the
146 * the out of range value.
147 */
148 return (spu_sel(vy, out_of_range, (vec_ullong2)spu_orc(overflow, underflow)));
149 }
150
151 #endif /* _EXP2D2_H_ */
152 #endif /* __SPU__ */
153