GNU Radio Manual and C++ API Reference  3.7.6.1
The Free & Open Software Radio Ecosystem
volk_32f_sin_32f.h
Go to the documentation of this file.
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2014 Free Software Foundation, Inc.
4  *
5  * This file is part of GNU Radio
6  *
7  * GNU Radio is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 3, or (at your option)
10  * any later version.
11  *
12  * GNU Radio is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with GNU Radio; see the file COPYING. If not, write to
19  * the Free Software Foundation, Inc., 51 Franklin Street,
20  * Boston, MA 02110-1301, USA.
21  */
22 
23 #include <stdio.h>
24 #include <math.h>
25 #include <inttypes.h>
26 
27 #ifndef INCLUDED_volk_32f_sin_32f_a_H
28 #define INCLUDED_volk_32f_sin_32f_a_H
29 
30 #ifdef LV_HAVE_SSE4_1
31 #include <smmintrin.h>
32 /*!
33  \brief Computes sine of input vector and stores results in output vector
34  \param bVector The vector where results will be stored
35  \param aVector The input vector of floats
36  \param num_points Number of points for which sine is to be computed
37 */
38 static inline void volk_32f_sin_32f_a_sse4_1(float* bVector, const float* aVector, unsigned int num_points){
39 
40  float* bPtr = bVector;
41  const float* aPtr = aVector;
42 
43  unsigned int number = 0;
44  unsigned int quarterPoints = num_points / 4;
45  unsigned int i = 0;
46 
47  __m128 aVal, s, m4pi, pio4A, pio4B, cp1, cp2, cp3, cp4, cp5, ffours, ftwos, fones, fzeroes;
48  __m128 sine, cosine, condition1, condition2;
49  __m128i q, r, ones, twos, fours;
50 
51  m4pi = _mm_set1_ps(1.273239545);
52  pio4A = _mm_set1_ps(0.78515625);
53  pio4B = _mm_set1_ps(0.241876e-3);
54  ffours = _mm_set1_ps(4.0);
55  ftwos = _mm_set1_ps(2.0);
56  fones = _mm_set1_ps(1.0);
57  fzeroes = _mm_setzero_ps();
58  ones = _mm_set1_epi32(1);
59  twos = _mm_set1_epi32(2);
60  fours = _mm_set1_epi32(4);
61 
62  cp1 = _mm_set1_ps(1.0);
63  cp2 = _mm_set1_ps(0.83333333e-1);
64  cp3 = _mm_set1_ps(0.2777778e-2);
65  cp4 = _mm_set1_ps(0.49603e-4);
66  cp5 = _mm_set1_ps(0.551e-6);
67 
68  for(;number < quarterPoints; number++){
69  aVal = _mm_load_ps(aPtr);
70  s = _mm_sub_ps(aVal, _mm_and_ps(_mm_mul_ps(aVal, ftwos), _mm_cmplt_ps(aVal, fzeroes)));
71  q = _mm_cvtps_epi32(_mm_mul_ps(s, m4pi));
72  r = _mm_add_epi32(q, _mm_and_si128(q, ones));
73 
74  s = _mm_sub_ps(s, _mm_mul_ps(_mm_cvtepi32_ps(r), pio4A));
75  s = _mm_sub_ps(s, _mm_mul_ps(_mm_cvtepi32_ps(r), pio4B));
76 
77  s = _mm_div_ps(s, _mm_set1_ps(8.0)); // The constant is 2^N, for 3 times argument reduction
78  s = _mm_mul_ps(s, s);
79  // Evaluate Taylor series
80  s = _mm_mul_ps(_mm_add_ps(_mm_mul_ps(_mm_sub_ps(_mm_mul_ps(_mm_add_ps(_mm_mul_ps(_mm_sub_ps(_mm_mul_ps(s, cp5), cp4), s), cp3), s), cp2), s), cp1), s);
81 
82  for(i = 0; i < 3; i++) {
83  s = _mm_mul_ps(s, _mm_sub_ps(ffours, s));
84  }
85  s = _mm_div_ps(s, ftwos);
86 
87  sine = _mm_sqrt_ps(_mm_mul_ps(_mm_sub_ps(ftwos, s), s));
88  cosine = _mm_sub_ps(fones, s);
89 
90  condition1 = _mm_cmpneq_ps(_mm_cvtepi32_ps(_mm_and_si128(_mm_add_epi32(q, ones), twos)), fzeroes);
91  condition2 = _mm_cmpneq_ps(_mm_cmpneq_ps(_mm_cvtepi32_ps(_mm_and_si128(q, fours)), fzeroes), _mm_cmplt_ps(aVal, fzeroes));
92  // Need this condition only for cos
93  //condition3 = _mm_cmpneq_ps(_mm_cvtepi32_ps(_mm_and_si128(_mm_add_epi32(q, twos), fours)), fzeroes);
94 
95  sine = _mm_add_ps(sine, _mm_and_ps(_mm_sub_ps(cosine, sine), condition1));
96  sine = _mm_sub_ps(sine, _mm_and_ps(_mm_mul_ps(sine, _mm_set1_ps(2.0f)), condition2));
97  _mm_store_ps(bPtr, sine);
98  aPtr += 4;
99  bPtr += 4;
100  }
101 
102  number = quarterPoints * 4;
103  for(;number < num_points; number++){
104  *bPtr++ = sin(*aPtr++);
105  }
106 }
107 #endif /* LV_HAVE_SSE4_1 for aligned */
108 
109 #endif /* INCLUDED_volk_32f_sin_32f_a_H */
110 
111 #ifndef INCLUDED_volk_32f_sin_32f_u_H
112 #define INCLUDED_volk_32f_sin_32f_u_H
113 
114 #ifdef LV_HAVE_SSE4_1
115 #include <smmintrin.h>
116 /*!
117  \brief Computes sine of input vector and stores results in output vector
118  \param bVector The vector where results will be stored
119  \param aVector The input vector of floats
120  \param num_points Number of points for which sine is to be computed
121 */
122 static inline void volk_32f_sin_32f_u_sse4_1(float* bVector, const float* aVector, unsigned int num_points){
123 
124  float* bPtr = bVector;
125  const float* aPtr = aVector;
126 
127  unsigned int number = 0;
128  unsigned int quarterPoints = num_points / 4;
129  unsigned int i = 0;
130 
131  __m128 aVal, s, m4pi, pio4A, pio4B, cp1, cp2, cp3, cp4, cp5, ffours, ftwos, fones, fzeroes;
132  __m128 sine, cosine, condition1, condition2;
133  __m128i q, r, ones, twos, fours;
134 
135  m4pi = _mm_set1_ps(1.273239545);
136  pio4A = _mm_set1_ps(0.78515625);
137  pio4B = _mm_set1_ps(0.241876e-3);
138  ffours = _mm_set1_ps(4.0);
139  ftwos = _mm_set1_ps(2.0);
140  fones = _mm_set1_ps(1.0);
141  fzeroes = _mm_setzero_ps();
142  ones = _mm_set1_epi32(1);
143  twos = _mm_set1_epi32(2);
144  fours = _mm_set1_epi32(4);
145 
146  cp1 = _mm_set1_ps(1.0);
147  cp2 = _mm_set1_ps(0.83333333e-1);
148  cp3 = _mm_set1_ps(0.2777778e-2);
149  cp4 = _mm_set1_ps(0.49603e-4);
150  cp5 = _mm_set1_ps(0.551e-6);
151 
152  for(;number < quarterPoints; number++){
153  aVal = _mm_loadu_ps(aPtr);
154  s = _mm_sub_ps(aVal, _mm_and_ps(_mm_mul_ps(aVal, ftwos), _mm_cmplt_ps(aVal, fzeroes)));
155  q = _mm_cvtps_epi32(_mm_mul_ps(s, m4pi));
156  r = _mm_add_epi32(q, _mm_and_si128(q, ones));
157 
158  s = _mm_sub_ps(s, _mm_mul_ps(_mm_cvtepi32_ps(r), pio4A));
159  s = _mm_sub_ps(s, _mm_mul_ps(_mm_cvtepi32_ps(r), pio4B));
160 
161  s = _mm_div_ps(s, _mm_set1_ps(8.0)); // The constant is 2^N, for 3 times argument reduction
162  s = _mm_mul_ps(s, s);
163  // Evaluate Taylor series
164  s = _mm_mul_ps(_mm_add_ps(_mm_mul_ps(_mm_sub_ps(_mm_mul_ps(_mm_add_ps(_mm_mul_ps(_mm_sub_ps(_mm_mul_ps(s, cp5), cp4), s), cp3), s), cp2), s), cp1), s);
165 
166  for(i = 0; i < 3; i++) {
167  s = _mm_mul_ps(s, _mm_sub_ps(ffours, s));
168  }
169  s = _mm_div_ps(s, ftwos);
170 
171  sine = _mm_sqrt_ps(_mm_mul_ps(_mm_sub_ps(ftwos, s), s));
172  cosine = _mm_sub_ps(fones, s);
173 
174  condition1 = _mm_cmpneq_ps(_mm_cvtepi32_ps(_mm_and_si128(_mm_add_epi32(q, ones), twos)), fzeroes);
175  condition2 = _mm_cmpneq_ps(_mm_cmpneq_ps(_mm_cvtepi32_ps(_mm_and_si128(q, fours)), fzeroes), _mm_cmplt_ps(aVal, fzeroes));
176 
177  sine = _mm_add_ps(sine, _mm_and_ps(_mm_sub_ps(cosine, sine), condition1));
178  sine = _mm_sub_ps(sine, _mm_and_ps(_mm_mul_ps(sine, _mm_set1_ps(2.0f)), condition2));
179  _mm_storeu_ps(bPtr, sine);
180  aPtr += 4;
181  bPtr += 4;
182  }
183 
184  number = quarterPoints * 4;
185  for(;number < num_points; number++){
186  *bPtr++ = sin(*aPtr++);
187  }
188 }
189 
190 #endif /* LV_HAVE_SSE4_1 for unaligned */
191 
192 #ifdef LV_HAVE_GENERIC
193 /*!
194  \brief Computes sine of input vector and stores results in output vector
195  \param bVector The vector where results will be stored
196  \param aVector The input vector of floats
197  \param num_points Number of points for which sine is to be computed
198 */
199 static inline void volk_32f_sin_32f_generic(float* bVector, const float* aVector, unsigned int num_points){
200  float* bPtr = bVector;
201  const float* aPtr = aVector;
202  unsigned int number = 0;
203 
204  for(number = 0; number < num_points; number++){
205  *bPtr++ = sin(*aPtr++);
206  }
207 
208 }
209 #endif /* LV_HAVE_GENERIC */
210 
211 #endif /* INCLUDED_volk_32f_sin_32f_u_H */