bullet: Update to current svn, r2636
[blender.git] / extern / bullet2 / src / LinearMath / btVector3.cpp
1 /*
2  Copyright (c) 2011 Apple Inc.
3  http://continuousphysics.com/Bullet/
4  
5  This software is provided 'as-is', without any express or implied warranty.
6  In no event will the authors be held liable for any damages arising from the use of this software.
7  Permission is granted to anyone to use this software for any purpose, 
8  including commercial applications, and to alter it and redistribute it freely, 
9  subject to the following restrictions:
10  
11  1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required.
12  2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
13  3. This notice may not be removed or altered from any source distribution.
14  
15  This source version has been altered.
16  */
17
18 #if defined (_WIN32) || defined (__i386__)
19 #define BT_USE_SSE_IN_API
20 #endif
21
22
23 #include "btVector3.h"
24
25
26
27 #if defined BT_USE_SIMD_VECTOR3
28
29 #if DEBUG
30 #include <string.h>//for memset
31 #endif
32
33
34 #ifdef __APPLE__
35 #include <stdint.h>
36 typedef  float float4 __attribute__ ((vector_size(16)));
37 #else
38 #define float4 __m128
39 #endif
40 //typedef  uint32_t uint4 __attribute__ ((vector_size(16)));
41
42
43 #if defined BT_USE_SSE || defined _WIN32
44
45 #define LOG2_ARRAY_SIZE     6
46 #define STACK_ARRAY_COUNT   (1UL << LOG2_ARRAY_SIZE)
47
48 #include <emmintrin.h>
49
50 long _maxdot_large( const float *vv, const float *vec, unsigned long count, float *dotResult );
51 long _maxdot_large( const float *vv, const float *vec, unsigned long count, float *dotResult )
52 {
53     const float4 *vertices = (const float4*) vv;
54     static const unsigned char indexTable[16] = {-1, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 };
55     float4 dotMax = btAssign128( -BT_INFINITY,  -BT_INFINITY,  -BT_INFINITY,  -BT_INFINITY );
56     float4 vvec = _mm_loadu_ps( vec );
57     float4 vHi = btCastiTo128f(_mm_shuffle_epi32( btCastfTo128i( vvec), 0xaa ));          /// zzzz
58     float4 vLo = _mm_movelh_ps( vvec, vvec );                               /// xyxy
59     
60     long maxIndex = -1L;
61     
62     size_t segment = 0;
63     float4 stack_array[ STACK_ARRAY_COUNT ];
64     
65 #if DEBUG
66     memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
67 #endif
68     
69     size_t index;
70     float4 max;
71     // Faster loop without cleanup code for full tiles
72     for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 ) 
73     {
74         max = dotMax;
75         
76         for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )   
77         { // do four dot products at a time. Carefully avoid touching the w element.
78             float4 v0 = vertices[0];
79             float4 v1 = vertices[1];
80             float4 v2 = vertices[2];
81             float4 v3 = vertices[3];            vertices += 4;
82             
83             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
84             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
85             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
86             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
87             
88             lo0 = lo0*vLo;
89             lo1 = lo1*vLo;
90             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
91             float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
92             float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
93             z = z*vHi;
94             x = x+y;
95             x = x+z;
96             stack_array[index] = x;
97             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
98             
99             v0 = vertices[0];
100             v1 = vertices[1];
101             v2 = vertices[2];
102             v3 = vertices[3];            vertices += 4;
103             
104             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
105             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
106             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
107             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
108             
109             lo0 = lo0*vLo;
110             lo1 = lo1*vLo;
111             z = _mm_shuffle_ps(hi0, hi1, 0x88);
112             x = _mm_shuffle_ps(lo0, lo1, 0x88);
113             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
114             z = z*vHi;
115             x = x+y;
116             x = x+z;
117             stack_array[index+1] = x;
118             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
119             
120             v0 = vertices[0];
121             v1 = vertices[1];
122             v2 = vertices[2];
123             v3 = vertices[3];            vertices += 4;
124             
125             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
126             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
127             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
128             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
129             
130             lo0 = lo0*vLo;
131             lo1 = lo1*vLo;
132             z = _mm_shuffle_ps(hi0, hi1, 0x88);
133             x = _mm_shuffle_ps(lo0, lo1, 0x88);
134             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
135             z = z*vHi;
136             x = x+y;
137             x = x+z;
138             stack_array[index+2] = x;
139             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
140             
141             v0 = vertices[0];
142             v1 = vertices[1];
143             v2 = vertices[2];
144             v3 = vertices[3];            vertices += 4;
145             
146             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
147             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
148             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
149             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
150             
151             lo0 = lo0*vLo;
152             lo1 = lo1*vLo;
153             z = _mm_shuffle_ps(hi0, hi1, 0x88);
154             x = _mm_shuffle_ps(lo0, lo1, 0x88);
155             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
156             z = z*vHi;
157             x = x+y;
158             x = x+z;
159             stack_array[index+3] = x;
160             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
161             
162             // It is too costly to keep the index of the max here. We will look for it again later.  We save a lot of work this way.
163         }
164         
165         // If we found a new max
166         if( 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(max, dotMax)))
167         { 
168             // copy the new max across all lanes of our max accumulator
169             max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0x4e));
170             max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0xb1));
171             
172             dotMax = max;
173             
174             // find first occurrence of that max  
175             size_t test;
176             for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], max))); index++ )   // local_count must be a multiple of 4
177             {}
178             // record where it is.
179             maxIndex = 4*index + segment + indexTable[test];
180         }
181     }
182     
183     // account for work we've already done
184     count -= segment;
185     
186     // Deal with the last < STACK_ARRAY_COUNT vectors
187     max = dotMax;
188     index = 0;
189     
190     
191     if( btUnlikely( count > 16) )
192     {
193         for( ; index + 4 <= count / 4; index+=4 )   
194         { // do four dot products at a time. Carefully avoid touching the w element.
195             float4 v0 = vertices[0];
196             float4 v1 = vertices[1];
197             float4 v2 = vertices[2];
198             float4 v3 = vertices[3];            vertices += 4;
199             
200             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
201             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
202             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
203             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
204             
205             lo0 = lo0*vLo;
206             lo1 = lo1*vLo;
207             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
208             float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
209             float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
210             z = z*vHi;
211             x = x+y;
212             x = x+z;
213             stack_array[index] = x;
214             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
215             
216             v0 = vertices[0];
217             v1 = vertices[1];
218             v2 = vertices[2];
219             v3 = vertices[3];            vertices += 4;
220             
221             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
222             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
223             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
224             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
225             
226             lo0 = lo0*vLo;
227             lo1 = lo1*vLo;
228             z = _mm_shuffle_ps(hi0, hi1, 0x88);
229             x = _mm_shuffle_ps(lo0, lo1, 0x88);
230             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
231             z = z*vHi;
232             x = x+y;
233             x = x+z;
234             stack_array[index+1] = x;
235             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
236             
237             v0 = vertices[0];
238             v1 = vertices[1];
239             v2 = vertices[2];
240             v3 = vertices[3];            vertices += 4;
241             
242             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
243             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
244             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
245             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
246             
247             lo0 = lo0*vLo;
248             lo1 = lo1*vLo;
249             z = _mm_shuffle_ps(hi0, hi1, 0x88);
250             x = _mm_shuffle_ps(lo0, lo1, 0x88);
251             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
252             z = z*vHi;
253             x = x+y;
254             x = x+z;
255             stack_array[index+2] = x;
256             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
257             
258             v0 = vertices[0];
259             v1 = vertices[1];
260             v2 = vertices[2];
261             v3 = vertices[3];            vertices += 4;
262             
263             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
264             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
265             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
266             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
267             
268             lo0 = lo0*vLo;
269             lo1 = lo1*vLo;
270             z = _mm_shuffle_ps(hi0, hi1, 0x88);
271             x = _mm_shuffle_ps(lo0, lo1, 0x88);
272             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
273             z = z*vHi;
274             x = x+y;
275             x = x+z;
276             stack_array[index+3] = x;
277             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
278             
279             // It is too costly to keep the index of the max here. We will look for it again later.  We save a lot of work this way.
280         }
281     }
282     
283     size_t localCount = (count & -4L) - 4*index;
284     if( localCount )
285     {
286 #ifdef __APPLE__
287         float4 t0, t1, t2, t3, t4;
288         float4 * sap = &stack_array[index + localCount / 4];
289           vertices += localCount;      // counter the offset
290          size_t byteIndex = -(localCount) * sizeof(float);
291         //AT&T Code style assembly
292         asm volatile
293         (   ".align 4                                                                   \n\
294              0: movaps  %[max], %[t2]                            // move max out of the way to avoid propagating NaNs in max \n\
295           movaps  (%[vertices], %[byteIndex], 4),    %[t0]    // vertices[0]      \n\
296           movaps  16(%[vertices], %[byteIndex], 4),  %[t1]    // vertices[1]      \n\
297           movaps  %[t0], %[max]                               // vertices[0]      \n\
298           movlhps %[t1], %[max]                               // x0y0x1y1         \n\
299          movaps  32(%[vertices], %[byteIndex], 4),  %[t3]    // vertices[2]      \n\
300          movaps  48(%[vertices], %[byteIndex], 4),  %[t4]    // vertices[3]      \n\
301           mulps   %[vLo], %[max]                              // x0y0x1y1 * vLo   \n\
302          movhlps %[t0], %[t1]                                // z0w0z1w1         \n\
303          movaps  %[t3], %[t0]                                // vertices[2]      \n\
304          movlhps %[t4], %[t0]                                // x2y2x3y3         \n\
305          mulps   %[vLo], %[t0]                               // x2y2x3y3 * vLo   \n\
306           movhlps %[t3], %[t4]                                // z2w2z3w3         \n\
307           shufps  $0x88, %[t4], %[t1]                         // z0z1z2z3         \n\
308           mulps   %[vHi], %[t1]                               // z0z1z2z3 * vHi   \n\
309          movaps  %[max], %[t3]                               // x0y0x1y1 * vLo   \n\
310          shufps  $0x88, %[t0], %[max]                        // x0x1x2x3 * vLo.x \n\
311          shufps  $0xdd, %[t0], %[t3]                         // y0y1y2y3 * vLo.y \n\
312          addps   %[t3], %[max]                               // x + y            \n\
313          addps   %[t1], %[max]                               // x + y + z        \n\
314          movaps  %[max], (%[sap], %[byteIndex])              // record result for later scrutiny \n\
315          maxps   %[t2], %[max]                               // record max, restore max   \n\
316          add     $16, %[byteIndex]                           // advance loop counter\n\
317          jnz     0b                                          \n\
318      "
319          : [max] "+x" (max), [t0] "=&x" (t0), [t1] "=&x" (t1), [t2] "=&x" (t2), [t3] "=&x" (t3), [t4] "=&x" (t4), [byteIndex] "+r" (byteIndex)
320          : [vLo] "x" (vLo), [vHi] "x" (vHi), [vertices] "r" (vertices), [sap] "r" (sap)
321          : "memory", "cc"
322          );
323         index += localCount/4;
324 #else
325         {
326             for( unsigned int i=0; i<localCount/4; i++,index++)   
327             { // do four dot products at a time. Carefully avoid touching the w element.
328                 float4 v0 = vertices[0];
329                 float4 v1 = vertices[1];
330                 float4 v2 = vertices[2];
331                 float4 v3 = vertices[3];            
332                 vertices += 4;
333                 
334                 float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
335                 float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
336                 float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
337                 float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
338                 
339                 lo0 = lo0*vLo;
340                 lo1 = lo1*vLo;
341                 float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
342                 float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
343                 float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
344                 z = z*vHi;
345                 x = x+y;
346                 x = x+z;
347                 stack_array[index] = x;
348                 max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
349             }
350         }
351 #endif //__APPLE__
352     }
353
354     // process the last few points
355     if( count & 3 )
356     {
357         float4 v0, v1, v2, x, y, z;
358         switch( count & 3 )
359         {
360             case 3:
361             {
362                 v0 = vertices[0];
363                 v1 = vertices[1];
364                 v2 = vertices[2];
365                 
366                 // Calculate 3 dot products, transpose, duplicate v2
367                 float4 lo0 = _mm_movelh_ps( v0, v1);        // xyxy.lo
368                 float4 hi0 = _mm_movehl_ps( v1, v0);        // z?z?.lo
369                 lo0 = lo0*vLo;
370                 z = _mm_shuffle_ps(hi0, v2,  0xa8 );           // z0z1z2z2
371                 z = z*vHi;
372                 float4 lo1 = _mm_movelh_ps(v2, v2);          // xyxy
373                 lo1 = lo1*vLo;
374                 x = _mm_shuffle_ps(lo0, lo1, 0x88);
375                 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
376             }
377                 break;
378             case 2:
379             {
380                 v0 = vertices[0];
381                 v1 = vertices[1];
382                 float4 xy = _mm_movelh_ps(v0, v1);
383                 z = _mm_movehl_ps(v1, v0);
384                 xy = xy*vLo;
385                 z = _mm_shuffle_ps( z, z,  0xa8);
386                 x = _mm_shuffle_ps( xy, xy, 0xa8);
387                 y = _mm_shuffle_ps( xy, xy, 0xfd);
388                 z = z*vHi;
389             }
390                 break;
391             case 1:
392             {
393                 float4 xy = vertices[0];
394                 z =  _mm_shuffle_ps( xy, xy, 0xaa);
395                 xy = xy*vLo;
396                 z = z*vHi;
397                 x = _mm_shuffle_ps(xy, xy, 0);
398                 y = _mm_shuffle_ps(xy, xy, 0x55);
399             }
400                 break;
401         }
402         x = x+y;
403         x = x+z;
404         stack_array[index] = x;
405         max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
406         index++;
407     }
408     
409     // if we found a new max. 
410     if( 0 == segment || 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(max, dotMax)))
411     { // we found a new max. Search for it
412       // find max across the max vector, place in all elements of max -- big latency hit here
413         max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0x4e));
414         max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0xb1));
415         
416         // It is slightly faster to do this part in scalar code when count < 8. However, the common case for
417         // this where it actually makes a difference is handled in the early out at the top of the function, 
418         // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced 
419         // complexity, and removed it.
420         
421         dotMax = max;
422         
423         // scan for the first occurence of max in the array  
424         size_t test;
425         for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], max))); index++ )   // local_count must be a multiple of 4
426         {}
427         maxIndex = 4*index + segment + indexTable[test];
428     }
429     
430     _mm_store_ss( dotResult, dotMax);
431     return maxIndex;
432 }
433
434 long _mindot_large( const float *vv, const float *vec, unsigned long count, float *dotResult );
435
436 long _mindot_large( const float *vv, const float *vec, unsigned long count, float *dotResult )
437 {
438     const float4 *vertices = (const float4*) vv;
439     static const unsigned char indexTable[16] = {-1, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 };
440     float4 dotmin = btAssign128( BT_INFINITY,  BT_INFINITY,  BT_INFINITY,  BT_INFINITY );
441     float4 vvec = _mm_loadu_ps( vec );
442     float4 vHi = btCastiTo128f(_mm_shuffle_epi32( btCastfTo128i( vvec), 0xaa ));          /// zzzz
443     float4 vLo = _mm_movelh_ps( vvec, vvec );                               /// xyxy
444     
445     long minIndex = -1L;
446
447     size_t segment = 0;
448     float4 stack_array[ STACK_ARRAY_COUNT ];
449     
450 #if DEBUG
451     memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
452 #endif
453     
454     size_t index;
455     float4 min;
456     // Faster loop without cleanup code for full tiles
457     for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 ) 
458     {
459         min = dotmin;
460         
461         for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )   
462         { // do four dot products at a time. Carefully avoid touching the w element.
463             float4 v0 = vertices[0];
464             float4 v1 = vertices[1];
465             float4 v2 = vertices[2];
466             float4 v3 = vertices[3];            vertices += 4;
467             
468             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
469             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
470             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
471             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
472             
473             lo0 = lo0*vLo;
474             lo1 = lo1*vLo;
475             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
476             float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
477             float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
478             z = z*vHi;
479             x = x+y;
480             x = x+z;
481             stack_array[index] = x;
482             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
483             
484             v0 = vertices[0];
485             v1 = vertices[1];
486             v2 = vertices[2];
487             v3 = vertices[3];            vertices += 4;
488             
489             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
490             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
491             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
492             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
493             
494             lo0 = lo0*vLo;
495             lo1 = lo1*vLo;
496             z = _mm_shuffle_ps(hi0, hi1, 0x88);
497             x = _mm_shuffle_ps(lo0, lo1, 0x88);
498             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
499             z = z*vHi;
500             x = x+y;
501             x = x+z;
502             stack_array[index+1] = x;
503             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
504             
505             v0 = vertices[0];
506             v1 = vertices[1];
507             v2 = vertices[2];
508             v3 = vertices[3];            vertices += 4;
509             
510             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
511             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
512             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
513             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
514             
515             lo0 = lo0*vLo;
516             lo1 = lo1*vLo;
517             z = _mm_shuffle_ps(hi0, hi1, 0x88);
518             x = _mm_shuffle_ps(lo0, lo1, 0x88);
519             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
520             z = z*vHi;
521             x = x+y;
522             x = x+z;
523             stack_array[index+2] = x;
524             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
525             
526             v0 = vertices[0];
527             v1 = vertices[1];
528             v2 = vertices[2];
529             v3 = vertices[3];            vertices += 4;
530             
531             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
532             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
533             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
534             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
535             
536             lo0 = lo0*vLo;
537             lo1 = lo1*vLo;
538             z = _mm_shuffle_ps(hi0, hi1, 0x88);
539             x = _mm_shuffle_ps(lo0, lo1, 0x88);
540             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
541             z = z*vHi;
542             x = x+y;
543             x = x+z;
544             stack_array[index+3] = x;
545             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
546             
547             // It is too costly to keep the index of the min here. We will look for it again later.  We save a lot of work this way.
548         }
549         
550         // If we found a new min
551         if( 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(min, dotmin)))
552         { 
553             // copy the new min across all lanes of our min accumulator
554             min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0x4e));
555             min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0xb1));
556             
557             dotmin = min;
558             
559             // find first occurrence of that min  
560             size_t test;
561             for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], min))); index++ )   // local_count must be a multiple of 4
562             {}
563             // record where it is.
564             minIndex = 4*index + segment + indexTable[test];
565         }
566     }
567     
568     // account for work we've already done
569     count -= segment;
570     
571     // Deal with the last < STACK_ARRAY_COUNT vectors
572     min = dotmin;
573     index = 0;
574     
575     
576     if(btUnlikely( count > 16) )
577     {
578         for( ; index + 4 <= count / 4; index+=4 )   
579         { // do four dot products at a time. Carefully avoid touching the w element.
580             float4 v0 = vertices[0];
581             float4 v1 = vertices[1];
582             float4 v2 = vertices[2];
583             float4 v3 = vertices[3];            vertices += 4;
584             
585             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
586             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
587             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
588             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
589             
590             lo0 = lo0*vLo;
591             lo1 = lo1*vLo;
592             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
593             float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
594             float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
595             z = z*vHi;
596             x = x+y;
597             x = x+z;
598             stack_array[index] = x;
599             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
600             
601             v0 = vertices[0];
602             v1 = vertices[1];
603             v2 = vertices[2];
604             v3 = vertices[3];            vertices += 4;
605             
606             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
607             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
608             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
609             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
610             
611             lo0 = lo0*vLo;
612             lo1 = lo1*vLo;
613             z = _mm_shuffle_ps(hi0, hi1, 0x88);
614             x = _mm_shuffle_ps(lo0, lo1, 0x88);
615             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
616             z = z*vHi;
617             x = x+y;
618             x = x+z;
619             stack_array[index+1] = x;
620             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
621             
622             v0 = vertices[0];
623             v1 = vertices[1];
624             v2 = vertices[2];
625             v3 = vertices[3];            vertices += 4;
626             
627             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
628             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
629             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
630             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
631             
632             lo0 = lo0*vLo;
633             lo1 = lo1*vLo;
634             z = _mm_shuffle_ps(hi0, hi1, 0x88);
635             x = _mm_shuffle_ps(lo0, lo1, 0x88);
636             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
637             z = z*vHi;
638             x = x+y;
639             x = x+z;
640             stack_array[index+2] = x;
641             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
642             
643             v0 = vertices[0];
644             v1 = vertices[1];
645             v2 = vertices[2];
646             v3 = vertices[3];            vertices += 4;
647             
648             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
649             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
650             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
651             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
652             
653             lo0 = lo0*vLo;
654             lo1 = lo1*vLo;
655             z = _mm_shuffle_ps(hi0, hi1, 0x88);
656             x = _mm_shuffle_ps(lo0, lo1, 0x88);
657             y = _mm_shuffle_ps(lo0, lo1, 0xdd);
658             z = z*vHi;
659             x = x+y;
660             x = x+z;
661             stack_array[index+3] = x;
662             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
663             
664             // It is too costly to keep the index of the min here. We will look for it again later.  We save a lot of work this way.
665         }
666     }
667     
668     size_t localCount = (count & -4L) - 4*index;
669     if( localCount )
670     {
671         
672         
673 #ifdef __APPLE__
674         vertices += localCount;      // counter the offset
675         float4 t0, t1, t2, t3, t4;
676         size_t byteIndex = -(localCount) * sizeof(float);
677         float4 * sap = &stack_array[index + localCount / 4];
678         
679         asm volatile
680         (   ".align 4                                                                   \n\
681              0: movaps  %[min], %[t2]                            // move min out of the way to avoid propagating NaNs in min \n\
682              movaps  (%[vertices], %[byteIndex], 4),    %[t0]    // vertices[0]      \n\
683              movaps  16(%[vertices], %[byteIndex], 4),  %[t1]    // vertices[1]      \n\
684              movaps  %[t0], %[min]                               // vertices[0]      \n\
685              movlhps %[t1], %[min]                               // x0y0x1y1         \n\
686              movaps  32(%[vertices], %[byteIndex], 4),  %[t3]    // vertices[2]      \n\
687              movaps  48(%[vertices], %[byteIndex], 4),  %[t4]    // vertices[3]      \n\
688              mulps   %[vLo], %[min]                              // x0y0x1y1 * vLo   \n\
689              movhlps %[t0], %[t1]                                // z0w0z1w1         \n\
690              movaps  %[t3], %[t0]                                // vertices[2]      \n\
691              movlhps %[t4], %[t0]                                // x2y2x3y3         \n\
692              movhlps %[t3], %[t4]                                // z2w2z3w3         \n\
693              mulps   %[vLo], %[t0]                               // x2y2x3y3 * vLo   \n\
694              shufps  $0x88, %[t4], %[t1]                         // z0z1z2z3         \n\
695              mulps   %[vHi], %[t1]                               // z0z1z2z3 * vHi   \n\
696              movaps  %[min], %[t3]                               // x0y0x1y1 * vLo   \n\
697              shufps  $0x88, %[t0], %[min]                        // x0x1x2x3 * vLo.x \n\
698              shufps  $0xdd, %[t0], %[t3]                         // y0y1y2y3 * vLo.y \n\
699              addps   %[t3], %[min]                               // x + y            \n\
700              addps   %[t1], %[min]                               // x + y + z        \n\
701              movaps  %[min], (%[sap], %[byteIndex])              // record result for later scrutiny \n\
702              minps   %[t2], %[min]                               // record min, restore min   \n\
703              add     $16, %[byteIndex]                           // advance loop counter\n\
704              jnz     0b                                          \n\
705              "
706          : [min] "+x" (min), [t0] "=&x" (t0), [t1] "=&x" (t1), [t2] "=&x" (t2), [t3] "=&x" (t3), [t4] "=&x" (t4), [byteIndex] "+r" (byteIndex)
707          : [vLo] "x" (vLo), [vHi] "x" (vHi), [vertices] "r" (vertices), [sap] "r" (sap)
708          : "memory", "cc"
709          );
710         index += localCount/4;
711 #else
712         {
713             for( unsigned int i=0; i<localCount/4; i++,index++)   
714             { // do four dot products at a time. Carefully avoid touching the w element.
715                 float4 v0 = vertices[0];
716                 float4 v1 = vertices[1];
717                 float4 v2 = vertices[2];
718                 float4 v3 = vertices[3];            
719                 vertices += 4;
720                 
721                 float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
722                 float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
723                 float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
724                 float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
725                 
726                 lo0 = lo0*vLo;
727                 lo1 = lo1*vLo;
728                 float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
729                 float4 x = _mm_shuffle_ps(lo0, lo1, 0x88);
730                 float4 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
731                 z = z*vHi;
732                 x = x+y;
733                 x = x+z;
734                 stack_array[index] = x;
735                 min = _mm_min_ps( x, min );         // control the order here so that max is never NaN even if x is nan
736             }
737         }
738
739 #endif
740     }
741     
742     // process the last few points
743     if( count & 3 )
744     {
745         float4 v0, v1, v2, x, y, z;
746         switch( count & 3 )
747         {
748             case 3:
749             {
750                 v0 = vertices[0];
751                 v1 = vertices[1];
752                 v2 = vertices[2];
753                 
754                 // Calculate 3 dot products, transpose, duplicate v2
755                 float4 lo0 = _mm_movelh_ps( v0, v1);        // xyxy.lo
756                 float4 hi0 = _mm_movehl_ps( v1, v0);        // z?z?.lo
757                 lo0 = lo0*vLo;
758                 z = _mm_shuffle_ps(hi0, v2,  0xa8 );           // z0z1z2z2
759                 z = z*vHi;
760                 float4 lo1 = _mm_movelh_ps(v2, v2);          // xyxy
761                 lo1 = lo1*vLo;
762                 x = _mm_shuffle_ps(lo0, lo1, 0x88);
763                 y = _mm_shuffle_ps(lo0, lo1, 0xdd);
764             }
765                 break;
766             case 2:
767             {
768                 v0 = vertices[0];
769                 v1 = vertices[1];
770                 float4 xy = _mm_movelh_ps(v0, v1);
771                 z = _mm_movehl_ps(v1, v0);
772                 xy = xy*vLo;
773                 z = _mm_shuffle_ps( z, z,  0xa8);
774                 x = _mm_shuffle_ps( xy, xy, 0xa8);
775                 y = _mm_shuffle_ps( xy, xy, 0xfd);
776                 z = z*vHi;
777             }
778                 break;
779             case 1:
780             {
781                 float4 xy = vertices[0];
782                 z =  _mm_shuffle_ps( xy, xy, 0xaa);
783                 xy = xy*vLo;
784                 z = z*vHi;
785                 x = _mm_shuffle_ps(xy, xy, 0);
786                 y = _mm_shuffle_ps(xy, xy, 0x55);
787             }
788                 break;
789         }
790         x = x+y;
791         x = x+z;
792         stack_array[index] = x;
793         min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
794         index++;
795     }
796     
797     // if we found a new min. 
798     if( 0 == segment || 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(min, dotmin)))
799     { // we found a new min. Search for it
800       // find min across the min vector, place in all elements of min -- big latency hit here
801         min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0x4e));
802         min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0xb1));
803         
804         // It is slightly faster to do this part in scalar code when count < 8. However, the common case for
805         // this where it actually makes a difference is handled in the early out at the top of the function, 
806         // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced 
807         // complexity, and removed it.
808         
809         dotmin = min;
810         
811         // scan for the first occurence of min in the array  
812         size_t test;
813         for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], min))); index++ )   // local_count must be a multiple of 4
814         {}
815         minIndex = 4*index + segment + indexTable[test];
816     }
817     
818     _mm_store_ss( dotResult, dotmin);
819     return minIndex;
820 }
821
822
823 #elif defined BT_USE_NEON
824 #define ARM_NEON_GCC_COMPATIBILITY  1
825 #include <arm_neon.h>
826
827
828 static long _maxdot_large_v0( const float *vv, const float *vec, unsigned long count, float *dotResult );
829 static long _maxdot_large_v1( const float *vv, const float *vec, unsigned long count, float *dotResult );
830 static long _maxdot_large_sel( const float *vv, const float *vec, unsigned long count, float *dotResult );
831 static long _mindot_large_v0( const float *vv, const float *vec, unsigned long count, float *dotResult );
832 static long _mindot_large_v1( const float *vv, const float *vec, unsigned long count, float *dotResult );
833 static long _mindot_large_sel( const float *vv, const float *vec, unsigned long count, float *dotResult );
834
835 long (*_maxdot_large)( const float *vv, const float *vec, unsigned long count, float *dotResult ) = _maxdot_large_sel;
836 long (*_mindot_large)( const float *vv, const float *vec, unsigned long count, float *dotResult ) = _mindot_large_sel;
837
838 extern "C" {int  _get_cpu_capabilities( void );}
839
840 static long _maxdot_large_sel( const float *vv, const float *vec, unsigned long count, float *dotResult )
841 {
842     if( _get_cpu_capabilities() & 0x2000 )
843         _maxdot_large = _maxdot_large_v1;
844     else
845         _maxdot_large = _maxdot_large_v0;
846     
847     return _maxdot_large(vv, vec, count, dotResult);
848 }
849
850 static long _mindot_large_sel( const float *vv, const float *vec, unsigned long count, float *dotResult )
851 {
852     if( _get_cpu_capabilities() & 0x2000 )
853         _mindot_large = _mindot_large_v1;
854     else
855         _mindot_large = _mindot_large_v0;
856     
857     return _mindot_large(vv, vec, count, dotResult);
858 }
859
860
861
862 #define vld1q_f32_aligned_postincrement( _ptr ) ({ float32x4_t _r; asm( "vld1.f32  {%0}, [%1, :128]!\n" : "=w" (_r), "+r" (_ptr) ); /*return*/ _r; })
863
864
865 long _maxdot_large_v0( const float *vv, const float *vec, unsigned long count, float *dotResult )
866 {
867     unsigned long i = 0;
868     float32x4_t vvec = vld1q_f32_aligned_postincrement( vec );
869     float32x2_t vLo = vget_low_f32(vvec);
870     float32x2_t vHi = vdup_lane_f32(vget_high_f32(vvec), 0);
871     float32x2_t dotMaxLo = (float32x2_t) { -BT_INFINITY, -BT_INFINITY };
872     float32x2_t dotMaxHi = (float32x2_t) { -BT_INFINITY, -BT_INFINITY };
873     uint32x2_t indexLo = (uint32x2_t) {0, 1};
874     uint32x2_t indexHi = (uint32x2_t) {2, 3};
875     uint32x2_t iLo = (uint32x2_t) {-1, -1};
876     uint32x2_t iHi = (uint32x2_t) {-1, -1};
877     const uint32x2_t four = (uint32x2_t) {4,4};
878
879     for( ; i+8 <= count; i+= 8 )
880     {
881         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
882         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
883         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
884         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
885         
886         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
887         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
888         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
889         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
890         
891         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
892         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
893         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
894         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
895         
896         float32x2_t rLo = vpadd_f32( xy0, xy1);
897         float32x2_t rHi = vpadd_f32( xy2, xy3);
898         rLo = vadd_f32(rLo, zLo);
899         rHi = vadd_f32(rHi, zHi);
900         
901         uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
902         uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
903         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
904         dotMaxHi = vbsl_f32( maskHi, rHi, dotMaxHi);
905         iLo = vbsl_u32(maskLo, indexLo, iLo);
906         iHi = vbsl_u32(maskHi, indexHi, iHi);
907         indexLo = vadd_u32(indexLo, four); 
908         indexHi = vadd_u32(indexHi, four);
909
910         v0 = vld1q_f32_aligned_postincrement( vv );
911         v1 = vld1q_f32_aligned_postincrement( vv );
912         v2 = vld1q_f32_aligned_postincrement( vv );
913         v3 = vld1q_f32_aligned_postincrement( vv );
914         
915         xy0 = vmul_f32( vget_low_f32(v0), vLo);
916         xy1 = vmul_f32( vget_low_f32(v1), vLo);
917         xy2 = vmul_f32( vget_low_f32(v2), vLo);
918         xy3 = vmul_f32( vget_low_f32(v3), vLo);
919         
920         z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
921         z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
922         zLo = vmul_f32( z0.val[0], vHi);
923         zHi = vmul_f32( z1.val[0], vHi);
924         
925         rLo = vpadd_f32( xy0, xy1);
926         rHi = vpadd_f32( xy2, xy3);
927         rLo = vadd_f32(rLo, zLo);
928         rHi = vadd_f32(rHi, zHi);
929         
930         maskLo = vcgt_f32( rLo, dotMaxLo );
931         maskHi = vcgt_f32( rHi, dotMaxHi );
932         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
933         dotMaxHi = vbsl_f32( maskHi, rHi, dotMaxHi);
934         iLo = vbsl_u32(maskLo, indexLo, iLo);
935         iHi = vbsl_u32(maskHi, indexHi, iHi);
936         indexLo = vadd_u32(indexLo, four);
937         indexHi = vadd_u32(indexHi, four);
938     }
939
940     for( ; i+4 <= count; i+= 4 )
941     {
942         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
943         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
944         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
945         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
946         
947         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
948         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
949         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
950         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
951         
952         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
953         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
954         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
955         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
956         
957         float32x2_t rLo = vpadd_f32( xy0, xy1);
958         float32x2_t rHi = vpadd_f32( xy2, xy3);
959         rLo = vadd_f32(rLo, zLo);
960         rHi = vadd_f32(rHi, zHi);
961         
962         uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
963         uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
964         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
965         dotMaxHi = vbsl_f32( maskHi, rHi, dotMaxHi);
966         iLo = vbsl_u32(maskLo, indexLo, iLo);
967         iHi = vbsl_u32(maskHi, indexHi, iHi);
968         indexLo = vadd_u32(indexLo, four);
969         indexHi = vadd_u32(indexHi, four);
970     }
971     
972     switch( count & 3 )
973     {
974         case 3:
975         {
976             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
977             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
978             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
979             
980             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
981             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
982             float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
983             
984             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
985             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
986             float32x2_t zHi = vmul_f32( vdup_lane_f32(vget_high_f32(v2), 0), vHi);
987             
988             float32x2_t rLo = vpadd_f32( xy0, xy1);
989             float32x2_t rHi = vpadd_f32( xy2, xy2);
990             rLo = vadd_f32(rLo, zLo);
991             rHi = vadd_f32(rHi, zHi);
992             
993             uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
994             uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
995             dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
996             dotMaxHi = vbsl_f32( maskHi, rHi, dotMaxHi);
997             iLo = vbsl_u32(maskLo, indexLo, iLo);
998             iHi = vbsl_u32(maskHi, indexHi, iHi);
999         }
1000             break;
1001         case 2:
1002         {
1003             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1004             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1005             
1006             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1007             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
1008             
1009             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1010             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
1011             
1012             float32x2_t rLo = vpadd_f32( xy0, xy1);
1013             rLo = vadd_f32(rLo, zLo);
1014             
1015             uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
1016             dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
1017             iLo = vbsl_u32(maskLo, indexLo, iLo);
1018         }
1019             break;
1020         case 1:
1021         {
1022             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1023             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1024             float32x2_t z0 = vdup_lane_f32(vget_high_f32(v0), 0);
1025             float32x2_t zLo = vmul_f32( z0, vHi);
1026             float32x2_t rLo = vpadd_f32( xy0, xy0);
1027             rLo = vadd_f32(rLo, zLo);
1028             uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
1029             dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
1030             iLo = vbsl_u32(maskLo, indexLo, iLo);
1031         }
1032             break;
1033         
1034         default:
1035             break;
1036     }
1037     
1038     // select best answer between hi and lo results
1039     uint32x2_t mask = vcgt_f32( dotMaxHi, dotMaxLo );
1040     dotMaxLo = vbsl_f32(mask, dotMaxHi, dotMaxLo);
1041     iLo = vbsl_u32(mask, iHi, iLo);
1042     
1043     // select best answer between even and odd results
1044     dotMaxHi = vdup_lane_f32(dotMaxLo, 1);
1045     iHi = vdup_lane_u32(iLo, 1);
1046     mask = vcgt_f32( dotMaxHi, dotMaxLo );
1047     dotMaxLo = vbsl_f32(mask, dotMaxHi, dotMaxLo);
1048     iLo = vbsl_u32(mask, iHi, iLo);
1049     
1050     *dotResult = vget_lane_f32( dotMaxLo, 0);
1051     return vget_lane_u32(iLo, 0);
1052 }
1053
1054
1055 long _maxdot_large_v1( const float *vv, const float *vec, unsigned long count, float *dotResult )
1056 {
1057     float32x4_t vvec = vld1q_f32_aligned_postincrement( vec );
1058     float32x4_t vLo = vcombine_f32(vget_low_f32(vvec), vget_low_f32(vvec));
1059     float32x4_t vHi = vdupq_lane_f32(vget_high_f32(vvec), 0);
1060     const uint32x4_t four = (uint32x4_t){ 4, 4, 4, 4 };
1061     uint32x4_t local_index = (uint32x4_t) {0, 1, 2, 3};
1062     uint32x4_t index = (uint32x4_t) { -1, -1, -1, -1 };
1063     float32x4_t maxDot = (float32x4_t) { -BT_INFINITY, -BT_INFINITY, -BT_INFINITY, -BT_INFINITY };
1064     
1065     unsigned long i = 0;
1066     for( ; i + 8 <= count; i += 8 )
1067     {
1068         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1069         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1070         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1071         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1072         
1073         // the next two lines should resolve to a single vswp d, d
1074         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1075         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1076         // the next two lines should resolve to a single vswp d, d
1077         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1078         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1079         
1080         xy0 = vmulq_f32(xy0, vLo);
1081         xy1 = vmulq_f32(xy1, vLo);
1082         
1083         float32x4x2_t zb = vuzpq_f32( z0, z1);
1084         float32x4_t z = vmulq_f32( zb.val[0], vHi);
1085         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1086         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1087         x = vaddq_f32(x, z);
1088         
1089         uint32x4_t mask = vcgtq_f32(x, maxDot);
1090         maxDot = vbslq_f32( mask, x, maxDot);
1091         index = vbslq_u32(mask, local_index, index);
1092         local_index = vaddq_u32(local_index, four);
1093
1094         v0 = vld1q_f32_aligned_postincrement( vv );
1095         v1 = vld1q_f32_aligned_postincrement( vv );
1096         v2 = vld1q_f32_aligned_postincrement( vv );
1097         v3 = vld1q_f32_aligned_postincrement( vv );
1098         
1099         // the next two lines should resolve to a single vswp d, d
1100         xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1101         xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1102         // the next two lines should resolve to a single vswp d, d
1103         z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1104         z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1105         
1106         xy0 = vmulq_f32(xy0, vLo);
1107         xy1 = vmulq_f32(xy1, vLo);
1108         
1109         zb = vuzpq_f32( z0, z1);
1110         z = vmulq_f32( zb.val[0], vHi);
1111         xy = vuzpq_f32( xy0, xy1);
1112         x = vaddq_f32(xy.val[0], xy.val[1]);
1113         x = vaddq_f32(x, z);
1114         
1115         mask = vcgtq_f32(x, maxDot);
1116         maxDot = vbslq_f32( mask, x, maxDot);
1117         index = vbslq_u32(mask, local_index, index);
1118         local_index = vaddq_u32(local_index, four);
1119     }
1120
1121     for( ; i + 4 <= count; i += 4 )
1122     {
1123         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1124         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1125         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1126         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1127
1128         // the next two lines should resolve to a single vswp d, d
1129         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1130         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1131         // the next two lines should resolve to a single vswp d, d
1132         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1133         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1134         
1135         xy0 = vmulq_f32(xy0, vLo);
1136         xy1 = vmulq_f32(xy1, vLo);
1137         
1138         float32x4x2_t zb = vuzpq_f32( z0, z1);
1139         float32x4_t z = vmulq_f32( zb.val[0], vHi);
1140         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1141         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1142         x = vaddq_f32(x, z);
1143         
1144         uint32x4_t mask = vcgtq_f32(x, maxDot);
1145         maxDot = vbslq_f32( mask, x, maxDot);
1146         index = vbslq_u32(mask, local_index, index);
1147         local_index = vaddq_u32(local_index, four);
1148     }
1149     
1150     switch (count & 3) {
1151         case 3:
1152         {
1153             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1154             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1155             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1156             
1157             // the next two lines should resolve to a single vswp d, d
1158             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1159             float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v2));
1160             // the next two lines should resolve to a single vswp d, d
1161             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1162             float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v2));
1163             
1164             xy0 = vmulq_f32(xy0, vLo);
1165             xy1 = vmulq_f32(xy1, vLo);
1166             
1167             float32x4x2_t zb = vuzpq_f32( z0, z1);
1168             float32x4_t z = vmulq_f32( zb.val[0], vHi);
1169             float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1170             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1171             x = vaddq_f32(x, z);
1172             
1173             uint32x4_t mask = vcgtq_f32(x, maxDot);
1174             maxDot = vbslq_f32( mask, x, maxDot);
1175             index = vbslq_u32(mask, local_index, index);
1176             local_index = vaddq_u32(local_index, four);
1177         }
1178             break;
1179
1180         case 2:
1181         {
1182             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1183             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1184             
1185             // the next two lines should resolve to a single vswp d, d
1186             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1187             // the next two lines should resolve to a single vswp d, d
1188             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1189             
1190             xy0 = vmulq_f32(xy0, vLo);
1191             
1192             float32x4x2_t zb = vuzpq_f32( z0, z0);
1193             float32x4_t z = vmulq_f32( zb.val[0], vHi);
1194             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
1195             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1196             x = vaddq_f32(x, z);
1197             
1198             uint32x4_t mask = vcgtq_f32(x, maxDot);
1199             maxDot = vbslq_f32( mask, x, maxDot);
1200             index = vbslq_u32(mask, local_index, index);
1201             local_index = vaddq_u32(local_index, four);
1202         }
1203             break;
1204
1205         case 1:
1206         {
1207             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1208             
1209             // the next two lines should resolve to a single vswp d, d
1210             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v0));
1211             // the next two lines should resolve to a single vswp d, d
1212             float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0); 
1213             
1214             xy0 = vmulq_f32(xy0, vLo);
1215             
1216             z = vmulq_f32( z, vHi);
1217             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
1218             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1219             x = vaddq_f32(x, z);
1220             
1221             uint32x4_t mask = vcgtq_f32(x, maxDot);
1222             maxDot = vbslq_f32( mask, x, maxDot);
1223             index = vbslq_u32(mask, local_index, index);
1224             local_index = vaddq_u32(local_index, four);
1225         }
1226             break;
1227
1228         default:
1229             break;
1230     }
1231     
1232     
1233     // select best answer between hi and lo results
1234     uint32x2_t mask = vcgt_f32( vget_high_f32(maxDot), vget_low_f32(maxDot));
1235     float32x2_t maxDot2 = vbsl_f32(mask, vget_high_f32(maxDot), vget_low_f32(maxDot));
1236     uint32x2_t index2 = vbsl_u32(mask, vget_high_u32(index), vget_low_u32(index));
1237     
1238     // select best answer between even and odd results
1239     float32x2_t maxDotO = vdup_lane_f32(maxDot2, 1);
1240     uint32x2_t indexHi = vdup_lane_u32(index2, 1);
1241     mask = vcgt_f32( maxDotO, maxDot2 );
1242     maxDot2 = vbsl_f32(mask, maxDotO, maxDot2);
1243     index2 = vbsl_u32(mask, indexHi, index2);
1244     
1245     *dotResult = vget_lane_f32( maxDot2, 0);
1246     return vget_lane_u32(index2, 0);
1247     
1248 }
1249
1250 long _mindot_large_v0( const float *vv, const float *vec, unsigned long count, float *dotResult )
1251 {
1252     unsigned long i = 0;
1253     float32x4_t vvec = vld1q_f32_aligned_postincrement( vec );
1254     float32x2_t vLo = vget_low_f32(vvec);
1255     float32x2_t vHi = vdup_lane_f32(vget_high_f32(vvec), 0);
1256     float32x2_t dotMinLo = (float32x2_t) { BT_INFINITY, BT_INFINITY };
1257     float32x2_t dotMinHi = (float32x2_t) { BT_INFINITY, BT_INFINITY };
1258     uint32x2_t indexLo = (uint32x2_t) {0, 1};
1259     uint32x2_t indexHi = (uint32x2_t) {2, 3};
1260     uint32x2_t iLo = (uint32x2_t) {-1, -1};
1261     uint32x2_t iHi = (uint32x2_t) {-1, -1};
1262     const uint32x2_t four = (uint32x2_t) {4,4};
1263     
1264     for( ; i+8 <= count; i+= 8 )
1265     {
1266         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1267         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1268         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1269         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1270         
1271         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1272         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
1273         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
1274         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
1275         
1276         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1277         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
1278         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
1279         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
1280         
1281         float32x2_t rLo = vpadd_f32( xy0, xy1);
1282         float32x2_t rHi = vpadd_f32( xy2, xy3);
1283         rLo = vadd_f32(rLo, zLo);
1284         rHi = vadd_f32(rHi, zHi);
1285         
1286         uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
1287         uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
1288         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1289         dotMinHi = vbsl_f32( maskHi, rHi, dotMinHi);
1290         iLo = vbsl_u32(maskLo, indexLo, iLo);
1291         iHi = vbsl_u32(maskHi, indexHi, iHi);
1292         indexLo = vadd_u32(indexLo, four);
1293         indexHi = vadd_u32(indexHi, four);
1294         
1295         v0 = vld1q_f32_aligned_postincrement( vv );
1296         v1 = vld1q_f32_aligned_postincrement( vv );
1297         v2 = vld1q_f32_aligned_postincrement( vv );
1298         v3 = vld1q_f32_aligned_postincrement( vv );
1299         
1300         xy0 = vmul_f32( vget_low_f32(v0), vLo);
1301         xy1 = vmul_f32( vget_low_f32(v1), vLo);
1302         xy2 = vmul_f32( vget_low_f32(v2), vLo);
1303         xy3 = vmul_f32( vget_low_f32(v3), vLo);
1304         
1305         z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1306         z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
1307         zLo = vmul_f32( z0.val[0], vHi);
1308         zHi = vmul_f32( z1.val[0], vHi);
1309         
1310         rLo = vpadd_f32( xy0, xy1);
1311         rHi = vpadd_f32( xy2, xy3);
1312         rLo = vadd_f32(rLo, zLo);
1313         rHi = vadd_f32(rHi, zHi);
1314         
1315         maskLo = vclt_f32( rLo, dotMinLo );
1316         maskHi = vclt_f32( rHi, dotMinHi );
1317         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1318         dotMinHi = vbsl_f32( maskHi, rHi, dotMinHi);
1319         iLo = vbsl_u32(maskLo, indexLo, iLo);
1320         iHi = vbsl_u32(maskHi, indexHi, iHi);
1321         indexLo = vadd_u32(indexLo, four);
1322         indexHi = vadd_u32(indexHi, four);
1323     }
1324
1325     for( ; i+4 <= count; i+= 4 )
1326     {
1327         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1328         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1329         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1330         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1331         
1332         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1333         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
1334         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
1335         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
1336         
1337         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1338         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
1339         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
1340         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
1341         
1342         float32x2_t rLo = vpadd_f32( xy0, xy1);
1343         float32x2_t rHi = vpadd_f32( xy2, xy3);
1344         rLo = vadd_f32(rLo, zLo);
1345         rHi = vadd_f32(rHi, zHi);
1346         
1347         uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
1348         uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
1349         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1350         dotMinHi = vbsl_f32( maskHi, rHi, dotMinHi);
1351         iLo = vbsl_u32(maskLo, indexLo, iLo);
1352         iHi = vbsl_u32(maskHi, indexHi, iHi);
1353         indexLo = vadd_u32(indexLo, four);
1354         indexHi = vadd_u32(indexHi, four);
1355     }
1356     switch( count & 3 )
1357     {
1358         case 3:
1359         {
1360             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1361             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1362             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1363             
1364             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1365             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
1366             float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
1367             
1368             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1369             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
1370             float32x2_t zHi = vmul_f32( vdup_lane_f32(vget_high_f32(v2), 0), vHi);
1371             
1372             float32x2_t rLo = vpadd_f32( xy0, xy1);
1373             float32x2_t rHi = vpadd_f32( xy2, xy2);
1374             rLo = vadd_f32(rLo, zLo);
1375             rHi = vadd_f32(rHi, zHi);
1376             
1377             uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
1378             uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
1379             dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1380             dotMinHi = vbsl_f32( maskHi, rHi, dotMinHi);
1381             iLo = vbsl_u32(maskLo, indexLo, iLo);
1382             iHi = vbsl_u32(maskHi, indexHi, iHi);
1383         }
1384             break;
1385         case 2:
1386         {
1387             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1388             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1389             
1390             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1391             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
1392             
1393             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
1394             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
1395             
1396             float32x2_t rLo = vpadd_f32( xy0, xy1);
1397             rLo = vadd_f32(rLo, zLo);
1398             
1399             uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
1400             dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1401             iLo = vbsl_u32(maskLo, indexLo, iLo);
1402         }
1403             break;
1404         case 1:
1405         {
1406             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1407             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
1408             float32x2_t z0 = vdup_lane_f32(vget_high_f32(v0), 0);
1409             float32x2_t zLo = vmul_f32( z0, vHi);
1410             float32x2_t rLo = vpadd_f32( xy0, xy0);
1411             rLo = vadd_f32(rLo, zLo);
1412             uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
1413             dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
1414             iLo = vbsl_u32(maskLo, indexLo, iLo);
1415         }
1416             break;
1417             
1418         default:
1419             break;
1420     }
1421     
1422     // select best answer between hi and lo results
1423     uint32x2_t mask = vclt_f32( dotMinHi, dotMinLo );
1424     dotMinLo = vbsl_f32(mask, dotMinHi, dotMinLo);
1425     iLo = vbsl_u32(mask, iHi, iLo);
1426     
1427     // select best answer between even and odd results
1428     dotMinHi = vdup_lane_f32(dotMinLo, 1);
1429     iHi = vdup_lane_u32(iLo, 1);
1430     mask = vclt_f32( dotMinHi, dotMinLo );
1431     dotMinLo = vbsl_f32(mask, dotMinHi, dotMinLo);
1432     iLo = vbsl_u32(mask, iHi, iLo);
1433     
1434     *dotResult = vget_lane_f32( dotMinLo, 0);
1435     return vget_lane_u32(iLo, 0);
1436 }
1437
1438 long _mindot_large_v1( const float *vv, const float *vec, unsigned long count, float *dotResult )
1439 {
1440     float32x4_t vvec = vld1q_f32_aligned_postincrement( vec );
1441     float32x4_t vLo = vcombine_f32(vget_low_f32(vvec), vget_low_f32(vvec));
1442     float32x4_t vHi = vdupq_lane_f32(vget_high_f32(vvec), 0);
1443     const uint32x4_t four = (uint32x4_t){ 4, 4, 4, 4 };
1444     uint32x4_t local_index = (uint32x4_t) {0, 1, 2, 3};
1445     uint32x4_t index = (uint32x4_t) { -1, -1, -1, -1 };
1446     float32x4_t minDot = (float32x4_t) { BT_INFINITY, BT_INFINITY, BT_INFINITY, BT_INFINITY };
1447     
1448     unsigned long i = 0;
1449     for( ; i + 8 <= count; i += 8 )
1450     {
1451         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1452         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1453         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1454         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1455         
1456         // the next two lines should resolve to a single vswp d, d
1457         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1458         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1459         // the next two lines should resolve to a single vswp d, d
1460         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1461         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1462         
1463         xy0 = vmulq_f32(xy0, vLo);
1464         xy1 = vmulq_f32(xy1, vLo);
1465         
1466         float32x4x2_t zb = vuzpq_f32( z0, z1);
1467         float32x4_t z = vmulq_f32( zb.val[0], vHi);
1468         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1469         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1470         x = vaddq_f32(x, z);
1471         
1472         uint32x4_t mask = vcltq_f32(x, minDot);
1473         minDot = vbslq_f32( mask, x, minDot);
1474         index = vbslq_u32(mask, local_index, index);
1475         local_index = vaddq_u32(local_index, four);
1476         
1477         v0 = vld1q_f32_aligned_postincrement( vv );
1478         v1 = vld1q_f32_aligned_postincrement( vv );
1479         v2 = vld1q_f32_aligned_postincrement( vv );
1480         v3 = vld1q_f32_aligned_postincrement( vv );
1481         
1482         // the next two lines should resolve to a single vswp d, d
1483         xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1484         xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1485         // the next two lines should resolve to a single vswp d, d
1486         z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1487         z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1488         
1489         xy0 = vmulq_f32(xy0, vLo);
1490         xy1 = vmulq_f32(xy1, vLo);
1491         
1492         zb = vuzpq_f32( z0, z1);
1493         z = vmulq_f32( zb.val[0], vHi);
1494         xy = vuzpq_f32( xy0, xy1);
1495         x = vaddq_f32(xy.val[0], xy.val[1]);
1496         x = vaddq_f32(x, z);
1497         
1498         mask = vcltq_f32(x, minDot);
1499         minDot = vbslq_f32( mask, x, minDot);
1500         index = vbslq_u32(mask, local_index, index);
1501         local_index = vaddq_u32(local_index, four);
1502     }
1503     
1504     for( ; i + 4 <= count; i += 4 )
1505     {
1506         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1507         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1508         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1509         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
1510         
1511         // the next two lines should resolve to a single vswp d, d
1512         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1513         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
1514         // the next two lines should resolve to a single vswp d, d
1515         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1516         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
1517         
1518         xy0 = vmulq_f32(xy0, vLo);
1519         xy1 = vmulq_f32(xy1, vLo);
1520         
1521         float32x4x2_t zb = vuzpq_f32( z0, z1);
1522         float32x4_t z = vmulq_f32( zb.val[0], vHi);
1523         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1524         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1525         x = vaddq_f32(x, z);
1526         
1527         uint32x4_t mask = vcltq_f32(x, minDot);
1528         minDot = vbslq_f32( mask, x, minDot);
1529         index = vbslq_u32(mask, local_index, index);
1530         local_index = vaddq_u32(local_index, four);
1531     }
1532     
1533     switch (count & 3) {
1534         case 3:
1535         {
1536             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1537             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1538             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
1539             
1540             // the next two lines should resolve to a single vswp d, d
1541             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1542             float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v2));
1543             // the next two lines should resolve to a single vswp d, d
1544             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1545             float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v2));
1546             
1547             xy0 = vmulq_f32(xy0, vLo);
1548             xy1 = vmulq_f32(xy1, vLo);
1549             
1550             float32x4x2_t zb = vuzpq_f32( z0, z1);
1551             float32x4_t z = vmulq_f32( zb.val[0], vHi);
1552             float32x4x2_t xy = vuzpq_f32( xy0, xy1);
1553             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1554             x = vaddq_f32(x, z);
1555             
1556             uint32x4_t mask = vcltq_f32(x, minDot);
1557             minDot = vbslq_f32( mask, x, minDot);
1558             index = vbslq_u32(mask, local_index, index);
1559             local_index = vaddq_u32(local_index, four);
1560         }
1561             break;
1562             
1563         case 2:
1564         {
1565             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1566             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
1567             
1568             // the next two lines should resolve to a single vswp d, d
1569             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
1570             // the next two lines should resolve to a single vswp d, d
1571             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
1572             
1573             xy0 = vmulq_f32(xy0, vLo);
1574             
1575             float32x4x2_t zb = vuzpq_f32( z0, z0);
1576             float32x4_t z = vmulq_f32( zb.val[0], vHi);
1577             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
1578             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1579             x = vaddq_f32(x, z);
1580             
1581             uint32x4_t mask = vcltq_f32(x, minDot);
1582             minDot = vbslq_f32( mask, x, minDot);
1583             index = vbslq_u32(mask, local_index, index);
1584             local_index = vaddq_u32(local_index, four);
1585         }
1586             break;
1587             
1588         case 1:
1589         {
1590             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
1591             
1592             // the next two lines should resolve to a single vswp d, d
1593             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v0));
1594             // the next two lines should resolve to a single vswp d, d
1595             float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0); 
1596             
1597             xy0 = vmulq_f32(xy0, vLo);
1598             
1599             z = vmulq_f32( z, vHi);
1600             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
1601             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
1602             x = vaddq_f32(x, z);
1603             
1604             uint32x4_t mask = vcltq_f32(x, minDot);
1605             minDot = vbslq_f32( mask, x, minDot);
1606             index = vbslq_u32(mask, local_index, index);
1607             local_index = vaddq_u32(local_index, four);
1608         }
1609             break;
1610             
1611         default:
1612             break;
1613     }
1614     
1615     
1616     // select best answer between hi and lo results
1617     uint32x2_t mask = vclt_f32( vget_high_f32(minDot), vget_low_f32(minDot));
1618     float32x2_t minDot2 = vbsl_f32(mask, vget_high_f32(minDot), vget_low_f32(minDot));
1619     uint32x2_t index2 = vbsl_u32(mask, vget_high_u32(index), vget_low_u32(index));
1620     
1621     // select best answer between even and odd results
1622     float32x2_t minDotO = vdup_lane_f32(minDot2, 1);
1623     uint32x2_t indexHi = vdup_lane_u32(index2, 1);
1624     mask = vclt_f32( minDotO, minDot2 );
1625     minDot2 = vbsl_f32(mask, minDotO, minDot2);
1626     index2 = vbsl_u32(mask, indexHi, index2);
1627     
1628     *dotResult = vget_lane_f32( minDot2, 0);
1629     return vget_lane_u32(index2, 0);
1630     
1631 }
1632
1633 #else
1634     #error Unhandled __APPLE__ arch
1635 #endif
1636
1637 #endif  /* __APPLE__ */
1638
1639