From ba0dd676aa7e1f57cdec1c02d653e9eb58b9325e Mon Sep 17 00:00:00 2001
From: Erwin Coumans <erwincoumans@replicant.com>
Date: Sun, 2 Nov 2014 14:12:54 -0800
Subject: [PATCH 1/2] comment-out the memset in DEBUG mode in btVector3.cpp and
 b3Vector3.cpp

---
 src/Bullet3Common/b3Vector3.cpp | 488 ++++++++++++++++----------------
 src/LinearMath/btVector3.cpp    |   4 +-
 2 files changed, 246 insertions(+), 246 deletions(-)

diff --git a/src/Bullet3Common/b3Vector3.cpp b/src/Bullet3Common/b3Vector3.cpp
index e43f13897..5f5ac4ac0 100644
--- a/src/Bullet3Common/b3Vector3.cpp
+++ b/src/Bullet3Common/b3Vector3.cpp
@@ -1,16 +1,16 @@
 /*
  Copyright (c) 2011-213 Apple Inc. http://bulletphysics.org
- 
+
  This software is provided 'as-is', without any express or implied warranty.
  In no event will the authors be held liable for any damages arising from the use of this software.
- Permission is granted to anyone to use this software for any purpose, 
- including commercial applications, and to alter it and redistribute it freely, 
+ Permission is granted to anyone to use this software for any purpose,
+ including commercial applications, and to alter it and redistribute it freely,
  subject to the following restrictions:
- 
+
  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.
  2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
  3. This notice may not be removed or altered from any source distribution.
- 
+
  This source version has been altered.
  */
 
@@ -47,35 +47,35 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
     float4 vvec = _mm_loadu_ps( vec );
     float4 vHi = b3CastiTo128f(_mm_shuffle_epi32( b3CastfTo128i( vvec), 0xaa ));          /// zzzz
     float4 vLo = _mm_movelh_ps( vvec, vvec );                               /// xyxy
-    
+
     long maxIndex = -1L;
-    
+
     size_t segment = 0;
     float4 stack_array[ STACK_ARRAY_COUNT ];
-    
+
 #if DEBUG
-    memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
+    // memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
 #endif
-    
+
     size_t index;
     float4 max;
     // Faster loop without cleanup code for full tiles
-    for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 ) 
+    for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 )
     {
         max = dotMax;
-        
-        for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )   
+
+        for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )
         { // do four dot products at a time. Carefully avoid touching the w element.
             float4 v0 = vertices[0];
             float4 v1 = vertices[1];
             float4 v2 = vertices[2];
             float4 v3 = vertices[3];            vertices += 4;
-            
+
             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -86,17 +86,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -107,17 +107,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+1] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -128,17 +128,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+2] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -149,20 +149,20 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+3] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             // 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.
         }
-        
+
         // If we found a new max
         if( 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(max, dotMax)))
-        { 
+        {
             // copy the new max across all lanes of our max accumulator
             max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0x4e));
             max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0xb1));
-            
+
             dotMax = max;
-            
-            // find first occurrence of that max  
+
+            // find first occurrence of that max
             size_t test;
             for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], max))); index++ )   // local_count must be a multiple of 4
             {}
@@ -170,29 +170,29 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             maxIndex = 4*index + segment + indexTable[test];
         }
     }
-    
+
     // account for work we've already done
     count -= segment;
-    
+
     // Deal with the last < STACK_ARRAY_COUNT vectors
     max = dotMax;
     index = 0;
-    
-    
+
+
     if( b3Unlikely( count > 16) )
     {
-        for( ; index + 4 <= count / 4; index+=4 )   
+        for( ; index + 4 <= count / 4; index+=4 )
         { // do four dot products at a time. Carefully avoid touching the w element.
             float4 v0 = vertices[0];
             float4 v1 = vertices[1];
             float4 v2 = vertices[2];
             float4 v3 = vertices[3];            vertices += 4;
-            
+
             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -203,17 +203,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -224,17 +224,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+1] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -245,17 +245,17 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+2] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -266,11 +266,11 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+3] = x;
             max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
-            
+
             // 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.
         }
     }
-    
+
     size_t localCount = (count & -4L) - 4*index;
     if( localCount )
     {
@@ -314,19 +314,19 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
         index += localCount/4;
 #else
         {
-            for( unsigned int i=0; i<localCount/4; i++,index++)   
+            for( unsigned int i=0; i<localCount/4; i++,index++)
             { // do four dot products at a time. Carefully avoid touching the w element.
                 float4 v0 = vertices[0];
                 float4 v1 = vertices[1];
                 float4 v2 = vertices[2];
-                float4 v3 = vertices[3];            
+                float4 v3 = vertices[3];
                 vertices += 4;
-                
+
                 float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
                 float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
                 float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
                 float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-                
+
                 lo0 = lo0*vLo;
                 lo1 = lo1*vLo;
                 float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -353,7 +353,7 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
                 v0 = vertices[0];
                 v1 = vertices[1];
                 v2 = vertices[2];
-                
+
                 // Calculate 3 dot products, transpose, duplicate v2
                 float4 lo0 = _mm_movelh_ps( v0, v1);        // xyxy.lo
                 float4 hi0 = _mm_movehl_ps( v1, v0);        // z?z?.lo
@@ -396,28 +396,28 @@ long b3_maxdot_large( const float *vv, const float *vec, unsigned long count, fl
         max = _mm_max_ps( x, max );         // control the order here so that max is never NaN even if x is nan
         index++;
     }
-    
-    // if we found a new max. 
+
+    // if we found a new max.
     if( 0 == segment || 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(max, dotMax)))
     { // we found a new max. Search for it
       // find max across the max vector, place in all elements of max -- big latency hit here
         max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0x4e));
         max = _mm_max_ps(max, (float4) _mm_shuffle_ps( max, max, 0xb1));
-        
+
         // It is slightly faster to do this part in scalar code when count < 8. However, the common case for
-        // this where it actually makes a difference is handled in the early out at the top of the function, 
-        // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced 
+        // this where it actually makes a difference is handled in the early out at the top of the function,
+        // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced
         // complexity, and removed it.
-        
+
         dotMax = max;
-        
-        // scan for the first occurence of max in the array  
+
+        // scan for the first occurence of max in the array
         size_t test;
         for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], max))); index++ )   // local_count must be a multiple of 4
         {}
         maxIndex = 4*index + segment + indexTable[test];
     }
-    
+
     _mm_store_ss( dotResult, dotMax);
     return maxIndex;
 }
@@ -433,35 +433,35 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
     float4 vvec = _mm_loadu_ps( vec );
     float4 vHi = b3CastiTo128f(_mm_shuffle_epi32( b3CastfTo128i( vvec), 0xaa ));          /// zzzz
     float4 vLo = _mm_movelh_ps( vvec, vvec );                               /// xyxy
-    
+
     long minIndex = -1L;
 
     size_t segment = 0;
     float4 stack_array[ STACK_ARRAY_COUNT ];
-    
+
 #if DEBUG
-    memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
+    // memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
 #endif
-    
+
     size_t index;
     float4 min;
     // Faster loop without cleanup code for full tiles
-    for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 ) 
+    for ( segment = 0; segment + STACK_ARRAY_COUNT*4 <= count; segment += STACK_ARRAY_COUNT*4 )
     {
         min = dotmin;
-        
-        for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )   
+
+        for( index = 0; index < STACK_ARRAY_COUNT; index+= 4 )
         { // do four dot products at a time. Carefully avoid touching the w element.
             float4 v0 = vertices[0];
             float4 v1 = vertices[1];
             float4 v2 = vertices[2];
             float4 v3 = vertices[3];            vertices += 4;
-            
+
             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -472,17 +472,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -493,17 +493,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+1] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -514,17 +514,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+2] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -535,20 +535,20 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+3] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             // 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.
         }
-        
+
         // If we found a new min
         if( 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(min, dotmin)))
-        { 
+        {
             // copy the new min across all lanes of our min accumulator
             min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0x4e));
             min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0xb1));
-            
+
             dotmin = min;
-            
-            // find first occurrence of that min  
+
+            // find first occurrence of that min
             size_t test;
             for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], min))); index++ )   // local_count must be a multiple of 4
             {}
@@ -556,29 +556,29 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             minIndex = 4*index + segment + indexTable[test];
         }
     }
-    
+
     // account for work we've already done
     count -= segment;
-    
+
     // Deal with the last < STACK_ARRAY_COUNT vectors
     min = dotmin;
     index = 0;
-    
-    
+
+
     if(b3Unlikely( count > 16) )
     {
-        for( ; index + 4 <= count / 4; index+=4 )   
+        for( ; index + 4 <= count / 4; index+=4 )
         { // do four dot products at a time. Carefully avoid touching the w element.
             float4 v0 = vertices[0];
             float4 v1 = vertices[1];
             float4 v2 = vertices[2];
             float4 v3 = vertices[3];            vertices += 4;
-            
+
             float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -589,17 +589,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -610,17 +610,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+1] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -631,17 +631,17 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+2] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             v0 = vertices[0];
             v1 = vertices[1];
             v2 = vertices[2];
             v3 = vertices[3];            vertices += 4;
-            
+
             lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
             hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
             lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
             hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-            
+
             lo0 = lo0*vLo;
             lo1 = lo1*vLo;
             z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -652,22 +652,22 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
             x = x+z;
             stack_array[index+3] = x;
             min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
-            
+
             // 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.
         }
     }
-    
+
     size_t localCount = (count & -4L) - 4*index;
     if( localCount )
     {
-        
-        
+
+
 #ifdef __APPLE__
         vertices += localCount;      // counter the offset
         float4 t0, t1, t2, t3, t4;
         size_t byteIndex = -(localCount) * sizeof(float);
         float4 * sap = &stack_array[index + localCount / 4];
-        
+
         asm volatile
         (   ".align 4                                                                   \n\
              0: movaps  %[min], %[t2]                            // move min out of the way to avoid propagating NaNs in min \n\
@@ -702,19 +702,19 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
         index += localCount/4;
 #else
         {
-            for( unsigned int i=0; i<localCount/4; i++,index++)   
+            for( unsigned int i=0; i<localCount/4; i++,index++)
             { // do four dot products at a time. Carefully avoid touching the w element.
                 float4 v0 = vertices[0];
                 float4 v1 = vertices[1];
                 float4 v2 = vertices[2];
-                float4 v3 = vertices[3];            
+                float4 v3 = vertices[3];
                 vertices += 4;
-                
+
                 float4 lo0 = _mm_movelh_ps( v0, v1);    // x0y0x1y1
                 float4 hi0 = _mm_movehl_ps( v1, v0);    // z0?0z1?1
                 float4 lo1 = _mm_movelh_ps( v2, v3);    // x2y2x3y3
                 float4 hi1 = _mm_movehl_ps( v3, v2);    // z2?2z3?3
-                
+
                 lo0 = lo0*vLo;
                 lo1 = lo1*vLo;
                 float4 z = _mm_shuffle_ps(hi0, hi1, 0x88);
@@ -730,7 +730,7 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
 
 #endif
     }
-    
+
     // process the last few points
     if( count & 3 )
     {
@@ -742,7 +742,7 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
                 v0 = vertices[0];
                 v1 = vertices[1];
                 v2 = vertices[2];
-                
+
                 // Calculate 3 dot products, transpose, duplicate v2
                 float4 lo0 = _mm_movelh_ps( v0, v1);        // xyxy.lo
                 float4 hi0 = _mm_movehl_ps( v1, v0);        // z?z?.lo
@@ -785,28 +785,28 @@ long b3_mindot_large( const float *vv, const float *vec, unsigned long count, fl
         min = _mm_min_ps( x, min );         // control the order here so that min is never NaN even if x is nan
         index++;
     }
-    
-    // if we found a new min. 
+
+    // if we found a new min.
     if( 0 == segment || 0xf != _mm_movemask_ps( (float4) _mm_cmpeq_ps(min, dotmin)))
     { // we found a new min. Search for it
       // find min across the min vector, place in all elements of min -- big latency hit here
         min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0x4e));
         min = _mm_min_ps(min, (float4) _mm_shuffle_ps( min, min, 0xb1));
-        
+
         // It is slightly faster to do this part in scalar code when count < 8. However, the common case for
-        // this where it actually makes a difference is handled in the early out at the top of the function, 
-        // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced 
+        // this where it actually makes a difference is handled in the early out at the top of the function,
+        // so it is less than a 1% difference here. I opted for improved code size, fewer branches and reduced
         // complexity, and removed it.
-        
+
         dotmin = min;
-        
-        // scan for the first occurence of min in the array  
+
+        // scan for the first occurence of min in the array
         size_t test;
         for( index = 0; 0 == (test=_mm_movemask_ps( _mm_cmpeq_ps( stack_array[index], min))); index++ )   // local_count must be a multiple of 4
         {}
         minIndex = 4*index + segment + indexTable[test];
     }
-    
+
     _mm_store_ss( dotResult, dotmin);
     return minIndex;
 }
@@ -835,7 +835,7 @@ static long b3_maxdot_large_sel( const float *vv, const float *vec, unsigned lon
         b3_maxdot_large = _maxdot_large_v1;
     else
         b3_maxdot_large = _maxdot_large_v0;
-    
+
     return b3_maxdot_large(vv, vec, count, dotResult);
 }
 
@@ -845,7 +845,7 @@ static long b3_mindot_large_sel( const float *vv, const float *vec, unsigned lon
         b3_mindot_large = _mindot_large_v1;
     else
         b3_mindot_large = _mindot_large_v0;
-    
+
     return b3_mindot_large(vv, vec, count, dotResult);
 }
 
@@ -874,51 +874,51 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
-        
+
         float32x2_t rLo = vpadd_f32( xy0, xy1);
         float32x2_t rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
         uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
         dotMaxHi = vbsl_f32( maskHi, rHi, dotMaxHi);
         iLo = vbsl_u32(maskLo, indexLo, iLo);
         iHi = vbsl_u32(maskHi, indexHi, iHi);
-        indexLo = vadd_u32(indexLo, four); 
+        indexLo = vadd_u32(indexLo, four);
         indexHi = vadd_u32(indexHi, four);
 
         v0 = vld1q_f32_aligned_postincrement( vv );
         v1 = vld1q_f32_aligned_postincrement( vv );
         v2 = vld1q_f32_aligned_postincrement( vv );
         v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         xy0 = vmul_f32( vget_low_f32(v0), vLo);
         xy1 = vmul_f32( vget_low_f32(v1), vLo);
         xy2 = vmul_f32( vget_low_f32(v2), vLo);
         xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         zLo = vmul_f32( z0.val[0], vHi);
         zHi = vmul_f32( z1.val[0], vHi);
-        
+
         rLo = vpadd_f32( xy0, xy1);
         rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         maskLo = vcgt_f32( rLo, dotMaxLo );
         maskHi = vcgt_f32( rHi, dotMaxHi );
         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
@@ -935,22 +935,22 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
-        
+
         float32x2_t rLo = vpadd_f32( xy0, xy1);
         float32x2_t rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
         uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
         dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
@@ -960,7 +960,7 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
         indexLo = vadd_u32(indexLo, four);
         indexHi = vadd_u32(indexHi, four);
     }
-    
+
     switch( count & 3 )
     {
         case 3:
@@ -968,20 +968,20 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
-            
+
             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
             float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
-            
+
             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
             float32x2_t zHi = vmul_f32( vdup_lane_f32(vget_high_f32(v2), 0), vHi);
-            
+
             float32x2_t rLo = vpadd_f32( xy0, xy1);
             float32x2_t rHi = vpadd_f32( xy2, xy2);
             rLo = vadd_f32(rLo, zLo);
             rHi = vadd_f32(rHi, zHi);
-            
+
             uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
             uint32x2_t maskHi = vcgt_f32( rHi, dotMaxHi );
             dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
@@ -994,16 +994,16 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
-            
+
             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
-            
+
             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
-            
+
             float32x2_t rLo = vpadd_f32( xy0, xy1);
             rLo = vadd_f32(rLo, zLo);
-            
+
             uint32x2_t maskLo = vcgt_f32( rLo, dotMaxLo );
             dotMaxLo = vbsl_f32( maskLo, rLo, dotMaxLo);
             iLo = vbsl_u32(maskLo, indexLo, iLo);
@@ -1022,23 +1022,23 @@ long b3_maxdot_large_v0( const float *vv, const float *vec, unsigned long count,
             iLo = vbsl_u32(maskLo, indexLo, iLo);
         }
             break;
-        
+
         default:
             break;
     }
-    
+
     // select best answer between hi and lo results
     uint32x2_t mask = vcgt_f32( dotMaxHi, dotMaxLo );
     dotMaxLo = vbsl_f32(mask, dotMaxHi, dotMaxLo);
     iLo = vbsl_u32(mask, iHi, iLo);
-    
+
     // select best answer between even and odd results
     dotMaxHi = vdup_lane_f32(dotMaxLo, 1);
     iHi = vdup_lane_u32(iLo, 1);
     mask = vcgt_f32( dotMaxHi, dotMaxLo );
     dotMaxLo = vbsl_f32(mask, dotMaxHi, dotMaxLo);
     iLo = vbsl_u32(mask, iHi, iLo);
-    
+
     *dotResult = vget_lane_f32( dotMaxLo, 0);
     return vget_lane_u32(iLo, 0);
 }
@@ -1053,7 +1053,7 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
     uint32x4_t local_index = (uint32x4_t) {0, 1, 2, 3};
     uint32x4_t index = (uint32x4_t) { -1, -1, -1, -1 };
     float32x4_t maxDot = (float32x4_t) { -B3_INFINITY, -B3_INFINITY, -B3_INFINITY, -B3_INFINITY };
-    
+
     unsigned long i = 0;
     for( ; i + 8 <= count; i += 8 )
     {
@@ -1061,23 +1061,23 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         // the next two lines should resolve to a single vswp d, d
         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
         // the next two lines should resolve to a single vswp d, d
         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         float32x4x2_t zb = vuzpq_f32( z0, z1);
         float32x4_t z = vmulq_f32( zb.val[0], vHi);
         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         uint32x4_t mask = vcgtq_f32(x, maxDot);
         maxDot = vbslq_f32( mask, x, maxDot);
         index = vbslq_u32(mask, local_index, index);
@@ -1087,23 +1087,23 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         v1 = vld1q_f32_aligned_postincrement( vv );
         v2 = vld1q_f32_aligned_postincrement( vv );
         v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         // the next two lines should resolve to a single vswp d, d
         xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
         xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
         // the next two lines should resolve to a single vswp d, d
         z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         zb = vuzpq_f32( z0, z1);
         z = vmulq_f32( zb.val[0], vHi);
         xy = vuzpq_f32( xy0, xy1);
         x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         mask = vcgtq_f32(x, maxDot);
         maxDot = vbslq_f32( mask, x, maxDot);
         index = vbslq_u32(mask, local_index, index);
@@ -1123,45 +1123,45 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         // the next two lines should resolve to a single vswp d, d
         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         float32x4x2_t zb = vuzpq_f32( z0, z1);
         float32x4_t z = vmulq_f32( zb.val[0], vHi);
         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         uint32x4_t mask = vcgtq_f32(x, maxDot);
         maxDot = vbslq_f32( mask, x, maxDot);
         index = vbslq_u32(mask, local_index, index);
         local_index = vaddq_u32(local_index, four);
     }
-    
+
     switch (count & 3) {
         case 3:
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
             float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v2));
             // the next two lines should resolve to a single vswp d, d
             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v2));
-            
+
             xy0 = vmulq_f32(xy0, vLo);
             xy1 = vmulq_f32(xy1, vLo);
-            
+
             float32x4x2_t zb = vuzpq_f32( z0, z1);
             float32x4_t z = vmulq_f32( zb.val[0], vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy1);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcgtq_f32(x, maxDot);
             maxDot = vbslq_f32( mask, x, maxDot);
             index = vbslq_u32(mask, local_index, index);
@@ -1173,20 +1173,20 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
             // the next two lines should resolve to a single vswp d, d
             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
-            
+
             xy0 = vmulq_f32(xy0, vLo);
-            
+
             float32x4x2_t zb = vuzpq_f32( z0, z0);
             float32x4_t z = vmulq_f32( zb.val[0], vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcgtq_f32(x, maxDot);
             maxDot = vbslq_f32( mask, x, maxDot);
             index = vbslq_u32(mask, local_index, index);
@@ -1197,19 +1197,19 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         case 1:
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v0));
             // the next two lines should resolve to a single vswp d, d
-            float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0); 
-            
+            float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0);
+
             xy0 = vmulq_f32(xy0, vLo);
-            
+
             z = vmulq_f32( z, vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcgtq_f32(x, maxDot);
             maxDot = vbslq_f32( mask, x, maxDot);
             index = vbslq_u32(mask, local_index, index);
@@ -1220,23 +1220,23 @@ long b3_maxdot_large_v1( const float *vv, const float *vec, unsigned long count,
         default:
             break;
     }
-    
-    
+
+
     // select best answer between hi and lo results
     uint32x2_t mask = vcgt_f32( vget_high_f32(maxDot), vget_low_f32(maxDot));
     float32x2_t maxDot2 = vbsl_f32(mask, vget_high_f32(maxDot), vget_low_f32(maxDot));
     uint32x2_t index2 = vbsl_u32(mask, vget_high_u32(index), vget_low_u32(index));
-    
+
     // select best answer between even and odd results
     float32x2_t maxDotO = vdup_lane_f32(maxDot2, 1);
     uint32x2_t indexHi = vdup_lane_u32(index2, 1);
     mask = vcgt_f32( maxDotO, maxDot2 );
     maxDot2 = vbsl_f32(mask, maxDotO, maxDot2);
     index2 = vbsl_u32(mask, indexHi, index2);
-    
+
     *dotResult = vget_lane_f32( maxDot2, 0);
     return vget_lane_u32(index2, 0);
-    
+
 }
 
 long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count, float *dotResult )
@@ -1252,29 +1252,29 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
     uint32x2_t iLo = (uint32x2_t) {-1, -1};
     uint32x2_t iHi = (uint32x2_t) {-1, -1};
     const uint32x2_t four = (uint32x2_t) {4,4};
-    
+
     for( ; i+8 <= count; i+= 8 )
     {
         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
-        
+
         float32x2_t rLo = vpadd_f32( xy0, xy1);
         float32x2_t rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
         uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
@@ -1283,27 +1283,27 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
         iHi = vbsl_u32(maskHi, indexHi, iHi);
         indexLo = vadd_u32(indexLo, four);
         indexHi = vadd_u32(indexHi, four);
-        
+
         v0 = vld1q_f32_aligned_postincrement( vv );
         v1 = vld1q_f32_aligned_postincrement( vv );
         v2 = vld1q_f32_aligned_postincrement( vv );
         v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         xy0 = vmul_f32( vget_low_f32(v0), vLo);
         xy1 = vmul_f32( vget_low_f32(v1), vLo);
         xy2 = vmul_f32( vget_low_f32(v2), vLo);
         xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         zLo = vmul_f32( z0.val[0], vHi);
         zHi = vmul_f32( z1.val[0], vHi);
-        
+
         rLo = vpadd_f32( xy0, xy1);
         rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         maskLo = vclt_f32( rLo, dotMinLo );
         maskHi = vclt_f32( rHi, dotMinHi );
         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
@@ -1320,22 +1320,22 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
         float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
         float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
         float32x2_t xy3 = vmul_f32( vget_low_f32(v3), vLo);
-        
+
         float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x2x2_t z1 = vtrn_f32( vget_high_f32(v2), vget_high_f32(v3));
         float32x2_t zLo = vmul_f32( z0.val[0], vHi);
         float32x2_t zHi = vmul_f32( z1.val[0], vHi);
-        
+
         float32x2_t rLo = vpadd_f32( xy0, xy1);
         float32x2_t rHi = vpadd_f32( xy2, xy3);
         rLo = vadd_f32(rLo, zLo);
         rHi = vadd_f32(rHi, zHi);
-        
+
         uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
         uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
         dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
@@ -1352,20 +1352,20 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
-            
+
             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
             float32x2_t xy2 = vmul_f32( vget_low_f32(v2), vLo);
-            
+
             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
             float32x2_t zHi = vmul_f32( vdup_lane_f32(vget_high_f32(v2), 0), vHi);
-            
+
             float32x2_t rLo = vpadd_f32( xy0, xy1);
             float32x2_t rHi = vpadd_f32( xy2, xy2);
             rLo = vadd_f32(rLo, zLo);
             rHi = vadd_f32(rHi, zHi);
-            
+
             uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
             uint32x2_t maskHi = vclt_f32( rHi, dotMinHi );
             dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
@@ -1378,16 +1378,16 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
-            
+
             float32x2_t xy0 = vmul_f32( vget_low_f32(v0), vLo);
             float32x2_t xy1 = vmul_f32( vget_low_f32(v1), vLo);
-            
+
             float32x2x2_t z0 = vtrn_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x2_t zLo = vmul_f32( z0.val[0], vHi);
-            
+
             float32x2_t rLo = vpadd_f32( xy0, xy1);
             rLo = vadd_f32(rLo, zLo);
-            
+
             uint32x2_t maskLo = vclt_f32( rLo, dotMinLo );
             dotMinLo = vbsl_f32( maskLo, rLo, dotMinLo);
             iLo = vbsl_u32(maskLo, indexLo, iLo);
@@ -1406,23 +1406,23 @@ long b3_mindot_large_v0( const float *vv, const float *vec, unsigned long count,
             iLo = vbsl_u32(maskLo, indexLo, iLo);
         }
             break;
-            
+
         default:
             break;
     }
-    
+
     // select best answer between hi and lo results
     uint32x2_t mask = vclt_f32( dotMinHi, dotMinLo );
     dotMinLo = vbsl_f32(mask, dotMinHi, dotMinLo);
     iLo = vbsl_u32(mask, iHi, iLo);
-    
+
     // select best answer between even and odd results
     dotMinHi = vdup_lane_f32(dotMinLo, 1);
     iHi = vdup_lane_u32(iLo, 1);
     mask = vclt_f32( dotMinHi, dotMinLo );
     dotMinLo = vbsl_f32(mask, dotMinHi, dotMinLo);
     iLo = vbsl_u32(mask, iHi, iLo);
-    
+
     *dotResult = vget_lane_f32( dotMinLo, 0);
     return vget_lane_u32(iLo, 0);
 }
@@ -1436,7 +1436,7 @@ long b3_mindot_large_v1( const float *vv, const float *vec, unsigned long count,
     uint32x4_t local_index = (uint32x4_t) {0, 1, 2, 3};
     uint32x4_t index = (uint32x4_t) { -1, -1, -1, -1 };
     float32x4_t minDot = (float32x4_t) { B3_INFINITY, B3_INFINITY, B3_INFINITY, B3_INFINITY };
-    
+
     unsigned long i = 0;
     for( ; i + 8 <= count; i += 8 )
     {
@@ -1444,182 +1444,182 @@ long b3_mindot_large_v1( const float *vv, const float *vec, unsigned long count,
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         // the next two lines should resolve to a single vswp d, d
         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
         // the next two lines should resolve to a single vswp d, d
         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         float32x4x2_t zb = vuzpq_f32( z0, z1);
         float32x4_t z = vmulq_f32( zb.val[0], vHi);
         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         uint32x4_t mask = vcltq_f32(x, minDot);
         minDot = vbslq_f32( mask, x, minDot);
         index = vbslq_u32(mask, local_index, index);
         local_index = vaddq_u32(local_index, four);
-        
+
         v0 = vld1q_f32_aligned_postincrement( vv );
         v1 = vld1q_f32_aligned_postincrement( vv );
         v2 = vld1q_f32_aligned_postincrement( vv );
         v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         // the next two lines should resolve to a single vswp d, d
         xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
         xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
         // the next two lines should resolve to a single vswp d, d
         z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         zb = vuzpq_f32( z0, z1);
         z = vmulq_f32( zb.val[0], vHi);
         xy = vuzpq_f32( xy0, xy1);
         x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         mask = vcltq_f32(x, minDot);
         minDot = vbslq_f32( mask, x, minDot);
         index = vbslq_u32(mask, local_index, index);
         local_index = vaddq_u32(local_index, four);
     }
-    
+
     for( ; i + 4 <= count; i += 4 )
     {
         float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
         float32x4_t v3 = vld1q_f32_aligned_postincrement( vv );
-        
+
         // the next two lines should resolve to a single vswp d, d
         float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
         float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v3));
         // the next two lines should resolve to a single vswp d, d
         float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
         float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v3));
-        
+
         xy0 = vmulq_f32(xy0, vLo);
         xy1 = vmulq_f32(xy1, vLo);
-        
+
         float32x4x2_t zb = vuzpq_f32( z0, z1);
         float32x4_t z = vmulq_f32( zb.val[0], vHi);
         float32x4x2_t xy = vuzpq_f32( xy0, xy1);
         float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
         x = vaddq_f32(x, z);
-        
+
         uint32x4_t mask = vcltq_f32(x, minDot);
         minDot = vbslq_f32( mask, x, minDot);
         index = vbslq_u32(mask, local_index, index);
         local_index = vaddq_u32(local_index, four);
     }
-    
+
     switch (count & 3) {
         case 3:
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v2 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
             float32x4_t xy1 = vcombine_f32( vget_low_f32(v2), vget_low_f32(v2));
             // the next two lines should resolve to a single vswp d, d
             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
             float32x4_t z1 = vcombine_f32( vget_high_f32(v2), vget_high_f32(v2));
-            
+
             xy0 = vmulq_f32(xy0, vLo);
             xy1 = vmulq_f32(xy1, vLo);
-            
+
             float32x4x2_t zb = vuzpq_f32( z0, z1);
             float32x4_t z = vmulq_f32( zb.val[0], vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy1);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcltq_f32(x, minDot);
             minDot = vbslq_f32( mask, x, minDot);
             index = vbslq_u32(mask, local_index, index);
             local_index = vaddq_u32(local_index, four);
         }
             break;
-            
+
         case 2:
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
             float32x4_t v1 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v1));
             // the next two lines should resolve to a single vswp d, d
             float32x4_t z0 = vcombine_f32( vget_high_f32(v0), vget_high_f32(v1));
-            
+
             xy0 = vmulq_f32(xy0, vLo);
-            
+
             float32x4x2_t zb = vuzpq_f32( z0, z0);
             float32x4_t z = vmulq_f32( zb.val[0], vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcltq_f32(x, minDot);
             minDot = vbslq_f32( mask, x, minDot);
             index = vbslq_u32(mask, local_index, index);
             local_index = vaddq_u32(local_index, four);
         }
             break;
-            
+
         case 1:
         {
             float32x4_t v0 = vld1q_f32_aligned_postincrement( vv );
-            
+
             // the next two lines should resolve to a single vswp d, d
             float32x4_t xy0 = vcombine_f32( vget_low_f32(v0), vget_low_f32(v0));
             // the next two lines should resolve to a single vswp d, d
-            float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0); 
-            
+            float32x4_t z = vdupq_lane_f32(vget_high_f32(v0), 0);
+
             xy0 = vmulq_f32(xy0, vLo);
-            
+
             z = vmulq_f32( z, vHi);
             float32x4x2_t xy = vuzpq_f32( xy0, xy0);
             float32x4_t x = vaddq_f32(xy.val[0], xy.val[1]);
             x = vaddq_f32(x, z);
-            
+
             uint32x4_t mask = vcltq_f32(x, minDot);
             minDot = vbslq_f32( mask, x, minDot);
             index = vbslq_u32(mask, local_index, index);
             local_index = vaddq_u32(local_index, four);
         }
             break;
-            
+
         default:
             break;
     }
-    
-    
+
+
     // select best answer between hi and lo results
     uint32x2_t mask = vclt_f32( vget_high_f32(minDot), vget_low_f32(minDot));
     float32x2_t minDot2 = vbsl_f32(mask, vget_high_f32(minDot), vget_low_f32(minDot));
     uint32x2_t index2 = vbsl_u32(mask, vget_high_u32(index), vget_low_u32(index));
-    
+
     // select best answer between even and odd results
     float32x2_t minDotO = vdup_lane_f32(minDot2, 1);
     uint32x2_t indexHi = vdup_lane_u32(index2, 1);
     mask = vclt_f32( minDotO, minDot2 );
     minDot2 = vbsl_f32(mask, minDotO, minDot2);
     index2 = vbsl_u32(mask, indexHi, index2);
-    
+
     *dotResult = vget_lane_f32( minDot2, 0);
     return vget_lane_u32(index2, 0);
-    
+
 }
 
 #else
diff --git a/src/LinearMath/btVector3.cpp b/src/LinearMath/btVector3.cpp
index 80c7055d5..e05bdccd6 100644
--- a/src/LinearMath/btVector3.cpp
+++ b/src/LinearMath/btVector3.cpp
@@ -63,7 +63,7 @@ long _maxdot_large( const float *vv, const float *vec, unsigned long count, floa
     float4 stack_array[ STACK_ARRAY_COUNT ];
     
 #if DEBUG
-    memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
+    //memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
 #endif
     
     size_t index;
@@ -448,7 +448,7 @@ long _mindot_large( const float *vv, const float *vec, unsigned long count, floa
     float4 stack_array[ STACK_ARRAY_COUNT ];
     
 #if DEBUG
-    memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
+    //memset( stack_array, -1, STACK_ARRAY_COUNT * sizeof(stack_array[0]) );
 #endif
     
     size_t index;

From eafbf817c371bf102583f3e50d9639c4d985affc Mon Sep 17 00:00:00 2001
From: erwin coumans <erwincoumans@google.com>
Date: Wed, 5 Nov 2014 08:58:24 -0800
Subject: [PATCH 2/2] fix MSVC 2013 compilation (std::min needs #include
 <algorithm> add example use of GJK+MPR in VoronoiFractureDemo

---
 Demos/VoronoiFractureDemo/CMakeLists.txt      |   6 +
 .../VoronoiFractureDemo.cpp                   |  47 ++-
 .../VoronoiFractureDemo/VoronoiFractureDemo.h |   2 +
 .../btConvexConvexMprAlgorithm.cpp            | 288 ++++++++++++++++++
 .../btConvexConvexMprAlgorithm.h              |  87 ++++++
 Demos/VoronoiFractureDemo/main.cpp            |   2 +-
 .../dynamics/fem/fem_add_plasticity_force.h   |   4 +-
 7 files changed, 429 insertions(+), 7 deletions(-)
 create mode 100644 Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.cpp
 create mode 100644 Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.h

diff --git a/Demos/VoronoiFractureDemo/CMakeLists.txt b/Demos/VoronoiFractureDemo/CMakeLists.txt
index 8b1a10963..2abff4ece 100644
--- a/Demos/VoronoiFractureDemo/CMakeLists.txt
+++ b/Demos/VoronoiFractureDemo/CMakeLists.txt
@@ -29,6 +29,8 @@ ADD_EXECUTABLE(AppVoronoiFractureDemo
 		main.cpp
 		VoronoiFractureDemo.cpp 
 		VoronoiFractureDemo.h
+		btConvexConvexMprAlgorithm.cpp
+		btConvexConvexMprAlgorithm.h
 		${BULLET_PHYSICS_SOURCE_DIR}/build3/bullet.rc
 	)
 ELSE()
@@ -36,6 +38,8 @@ ELSE()
 		main.cpp
 		VoronoiFractureDemo.cpp 
 		VoronoiFractureDemo.h
+		btConvexConvexMprAlgorithm.cpp
+		btConvexConvexMprAlgorithm.h
 	)
 ENDIF()
 
@@ -57,6 +61,8 @@ ELSE (USE_GLUT)
 		Win32VoronoiFractureDemo.cpp
 		VoronoiFractureDemo.cpp 
 		VoronoiFractureDemo.h
+		btConvexConvexMprAlgorithm.cpp
+		btConvexConvexMprAlgorithm.h
 		${BULLET_PHYSICS_SOURCE_DIR}/build3/bullet.rc
 	)
 	
diff --git a/Demos/VoronoiFractureDemo/VoronoiFractureDemo.cpp b/Demos/VoronoiFractureDemo/VoronoiFractureDemo.cpp
index d920a6dff..19cf7ddd2 100644
--- a/Demos/VoronoiFractureDemo/VoronoiFractureDemo.cpp
+++ b/Demos/VoronoiFractureDemo/VoronoiFractureDemo.cpp
@@ -25,20 +25,23 @@ Voronoi fracture and shatter code and demo copyright (c) 2011 Alain Ducharme
 
 //maximum number of objects (and allow user to shoot additional boxes)
 #define MAX_PROXIES (2048)
-#define BREAKING_THRESHOLD 2
+#define BREAKING_THRESHOLD 3
 #define CONVEX_MARGIN 0.04
+static int useMpr = 0;
 
 #include "VoronoiFractureDemo.h"
 #include "GlutStuff.h"
 ///btBulletDynamicsCommon.h is the main Bullet include file, contains most common include files.
 #include "btBulletDynamicsCommon.h"
-
+#include "GLDebugFont.h"
 #include <stdio.h> //printf debugging
 
 #include "GLDebugDrawer.h"
 static GLDebugDrawer sDebugDraw;
 static bool useGenericConstraint = false;
 
+#include "btConvexConvexMprAlgorithm.h"
+
 
 void VoronoiFractureDemo::attachFixedConstraints()
 {
@@ -493,13 +496,13 @@ void VoronoiFractureDemo::clientMoveAndDisplay()
 	///step the simulation
 	if (m_dynamicsWorld)
 	{
-		m_dynamicsWorld->stepSimulation(ms / 1000000.f);
+		m_dynamicsWorld->stepSimulation(1. / 60., 0);// ms / 1000000.f);
 		//optional but useful: debug drawing
 		m_dynamicsWorld->debugDrawWorld();
 	}
 		
 	renderme(); 
-
+	
 	glFlush();
 
 	swapBuffers();
@@ -520,9 +523,30 @@ void VoronoiFractureDemo::displayCallback(void) {
 	swapBuffers();
 }
 
+void VoronoiFractureDemo::renderme()
+{
+	DemoApplication::renderme();
+	char buf[124];
+
+	int lineWidth = 200;
+	int xStart = m_glutScreenWidth - lineWidth;
+
+	if (useMpr)
+	{
+		sprintf(buf, "Using GJK+MPR");
+	}
+	else
+	{
+		sprintf(buf, "Using GJK+EPA");
+	}
+	GLDebugDrawString(xStart, 20, buf);
+
+}
+
 
 void	VoronoiFractureDemo::initPhysics()
 {
+	srand(13);
 	useGenericConstraint = !useGenericConstraint;
 	printf("useGenericConstraint = %d\n", useGenericConstraint);
 
@@ -538,7 +562,22 @@ void	VoronoiFractureDemo::initPhysics()
 
 	///use the default collision dispatcher. For parallel processing you can use a diffent dispatcher (see Extras/BulletMultiThreaded)
 	m_dispatcher = new	btCollisionDispatcher(m_collisionConfiguration);
+	
+	useMpr = 1 - useMpr;
 
+	if (useMpr)
+	{
+		printf("using GJK+MPR convex-convex collision detection\n");
+		btConvexConvexMprAlgorithm::CreateFunc* cf = new btConvexConvexMprAlgorithm::CreateFunc;
+		m_dispatcher->registerCollisionCreateFunc(CONVEX_HULL_SHAPE_PROXYTYPE, CONVEX_HULL_SHAPE_PROXYTYPE, cf);
+		m_dispatcher->registerCollisionCreateFunc(CONVEX_HULL_SHAPE_PROXYTYPE, BOX_SHAPE_PROXYTYPE, cf);
+		m_dispatcher->registerCollisionCreateFunc(BOX_SHAPE_PROXYTYPE, CONVEX_HULL_SHAPE_PROXYTYPE, cf);
+	}
+	else
+	{
+		printf("using default (GJK+EPA) convex-convex collision detection\n");
+	}
+	
 	m_broadphase = new btDbvtBroadphase();
 
 	///the default constraint solver. For parallel processing you can use a different solver (see Extras/BulletMultiThreaded)
diff --git a/Demos/VoronoiFractureDemo/VoronoiFractureDemo.h b/Demos/VoronoiFractureDemo/VoronoiFractureDemo.h
index a21737426..3c36819be 100644
--- a/Demos/VoronoiFractureDemo/VoronoiFractureDemo.h
+++ b/Demos/VoronoiFractureDemo/VoronoiFractureDemo.h
@@ -68,6 +68,8 @@ class VoronoiFractureDemo : public PlatformDemoApplication
 
 	void	exitPhysics();
 
+	virtual void renderme();
+
 	void getVerticesInsidePlanes(const btAlignedObjectArray<btVector3>& planes, btAlignedObjectArray<btVector3>& verticesOut, std::set<int>& planeIndicesOut);
 	void voronoiBBShatter(const btAlignedObjectArray<btVector3>& points, const btVector3& bbmin, const btVector3& bbmax, const btQuaternion& bbq, const btVector3& bbt, btScalar matDensity);
 	void voronoiConvexHullShatter(const btAlignedObjectArray<btVector3>& points, const btAlignedObjectArray<btVector3>& verts, const btQuaternion& bbq, const btVector3& bbt, btScalar matDensity);
diff --git a/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.cpp b/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.cpp
new file mode 100644
index 000000000..e14a0e67d
--- /dev/null
+++ b/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.cpp
@@ -0,0 +1,288 @@
+/*
+Bullet Continuous Collision Detection and Physics Library
+Copyright (c) 2003-2014 Erwin Coumans  http://continuousphysics.com/Bullet/
+
+This software is provided 'as-is', without any express or implied warranty.
+In no event will the authors be held liable for any damages arising from the use of this software.
+Permission is granted to anyone to use this software for any purpose, 
+including commercial applications, and to alter it and redistribute it freely, 
+subject to the following restrictions:
+
+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.
+2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+3. This notice may not be removed or altered from any source distribution.
+*/
+
+
+#include "btConvexConvexMprAlgorithm.h"
+
+//#include <stdio.h>
+#include "BulletCollision/NarrowPhaseCollision/btDiscreteCollisionDetectorInterface.h"
+#include "BulletCollision/BroadphaseCollision/btBroadphaseInterface.h"
+#include "BulletCollision/CollisionDispatch/btCollisionObject.h"
+#include "BulletCollision/CollisionShapes/btConvexShape.h"
+
+#include "BulletCollision/CollisionShapes/btTriangleShape.h"
+
+
+
+#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
+#include "BulletCollision/BroadphaseCollision/btBroadphaseProxy.h"
+#include "BulletCollision/CollisionDispatch/btCollisionDispatcher.h"
+#include "BulletCollision/CollisionShapes/btBoxShape.h"
+#include "BulletCollision/CollisionDispatch/btManifoldResult.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btConvexPenetrationDepthSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btContinuousConvexCollision.h"
+#include "BulletCollision/NarrowPhaseCollision/btSubSimplexConvexCast.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkConvexCast.h"
+
+
+
+#include "BulletCollision/NarrowPhaseCollision/btVoronoiSimplexSolver.h"
+#include "BulletCollision/CollisionShapes/btSphereShape.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btMinkowskiPenetrationDepthSolver.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpa2.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpaPenetrationDepthSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btPolyhedralContactClipping.h"
+#include "BulletCollision/CollisionDispatch/btCollisionObjectWrapper.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btComputeGjkEpaPenetration.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpa3.h"
+#include "BulletCollision/NarrowPhaseCollision/btMprPenetration.h"
+
+//this is just an internal debug variable to switch between GJK+MPR or GJK+EPA
+bool gUseMprCollisionFunction = true;
+
+btConvexConvexMprAlgorithm::CreateFunc::CreateFunc()
+{
+
+}
+
+btConvexConvexMprAlgorithm::CreateFunc::~CreateFunc() 
+{ 
+}
+
+btConvexConvexMprAlgorithm::btConvexConvexMprAlgorithm(btPersistentManifold* mf,const btCollisionAlgorithmConstructionInfo& ci,const btCollisionObjectWrapper* body0Wrap,const btCollisionObjectWrapper* body1Wrap)
+: btActivatingCollisionAlgorithm(ci,body0Wrap,body1Wrap),
+m_ownManifold (false),
+m_manifoldPtr(mf)
+{
+	(void)body0Wrap;
+	(void)body1Wrap;
+}
+
+
+
+
+btConvexConvexMprAlgorithm::~btConvexConvexMprAlgorithm()
+{
+	if (m_ownManifold)
+	{
+		if (m_manifoldPtr)
+			m_dispatcher->releaseManifold(m_manifoldPtr);
+	}
+}
+
+
+btVector3 btBulletShapeSupportFunc(const void* shapeAptr, const btVector3& dir, bool includeMargin)
+{
+	btConvexShape* shape = (btConvexShape*) shapeAptr;
+	if (includeMargin)
+	{
+		return shape->localGetSupportingVertex(dir);
+	}
+	
+	return shape->localGetSupportingVertexWithoutMargin(dir);
+}
+
+btVector3 btBulletShapeCenterFunc(const void* shapeAptr)
+{
+	return btVector3(0,0,0);
+}
+
+
+struct btMprConvexWrap
+{
+	const btConvexShape* m_convex;
+	btTransform m_worldTrans;
+	inline btScalar getMargin() const
+	{
+		return m_convex->getMargin();
+	}
+	inline btVector3 getObjectCenterInWorld() const
+	{
+		return m_worldTrans.getOrigin();
+	}
+	inline const btTransform& getWorldTransform() const
+	{
+		return m_worldTrans;
+	}
+	inline btVector3 getLocalSupportWithMargin(const btVector3& dir) const
+	{
+		return m_convex->localGetSupportingVertex(dir);
+	}
+	inline btVector3 getLocalSupportWithoutMargin(const btVector3& dir) const
+	{
+		return m_convex->localGetSupportingVertexWithoutMargin(dir);
+	}
+};
+
+struct btMyDistanceInfo
+{
+	btVector3	m_pointOnA;
+	btVector3	m_pointOnB;
+	btVector3	m_normalBtoA;
+	btScalar	m_distance;
+};
+
+//
+// Convex-Convex collision algorithm
+//
+void btConvexConvexMprAlgorithm ::processCollision (const btCollisionObjectWrapper* body0Wrap,const btCollisionObjectWrapper* body1Wrap,const btDispatcherInfo& dispatchInfo,btManifoldResult* resultOut)
+{
+
+	if (!m_manifoldPtr)
+	{
+		//swapped?
+		m_manifoldPtr = m_dispatcher->getNewManifold(body0Wrap->getCollisionObject(),body1Wrap->getCollisionObject());
+		m_ownManifold = true;
+	}
+	resultOut->setPersistentManifold(m_manifoldPtr);
+
+	//comment-out next line to test multi-contact generation
+	//resultOut->getPersistentManifold()->clearManifold();
+	
+
+	const btConvexShape* min0 = static_cast<const btConvexShape*>(body0Wrap->getCollisionShape());
+	const btConvexShape* min1 = static_cast<const btConvexShape*>(body1Wrap->getCollisionShape());
+
+	btVector3  normalOnB;
+	btVector3  pointOnBWorld;
+
+	btGjkPairDetector::ClosestPointInput input;
+
+	btVoronoiSimplexSolver vs;
+	btGjkEpaPenetrationDepthSolver epa;
+	
+	
+	if (gUseMprCollisionFunction)
+	{
+
+		btMprConvexWrap a,b;
+		a.m_worldTrans = body0Wrap->getWorldTransform();
+		b.m_worldTrans = body1Wrap->getWorldTransform();
+		a.m_convex = (const btConvexShape*)body0Wrap->getCollisionShape();
+		b.m_convex = (const btConvexShape*)body1Wrap->getCollisionShape();
+		btVoronoiSimplexSolver simplexSolver;
+		simplexSolver.reset();
+		btGjkCollisionDescription colDesc;
+		btMyDistanceInfo distInfo;
+		int res = btComputeGjkDistance(a,b,colDesc,&distInfo);
+		if (res==0)
+		{
+			   //printf("use GJK results in distance %f\n",distInfo.m_distance);
+		} else
+		{
+			btMprCollisionDescription mprDesc;
+			res = btComputeMprPenetration(a,b,mprDesc, &distInfo);
+
+			//printf("use MPR results in distance %f\n",distInfo.m_distance);
+
+		}
+		if (res == 0)
+		{
+#if 0
+			printf("Dist=%f,normalOnB[%f,%f,%f],pA=[%f,%f,%f],pB[%f,%f,%f]\n",
+				distInfo.m_distance, distInfo.m_normalBtoA[0], distInfo.m_normalBtoA[1], distInfo.m_normalBtoA[2],
+				distInfo.m_pointOnA[0], distInfo.m_pointOnA[1], distInfo.m_pointOnA[2],
+				distInfo.m_pointOnB[0], distInfo.m_pointOnB[1], distInfo.m_pointOnB[2]);
+#endif
+
+			if (distInfo.m_distance<=0)
+			{
+				resultOut->addContactPoint(distInfo.m_normalBtoA, distInfo.m_pointOnB, distInfo.m_distance);
+			}
+			//ASSERT_EQ(0,result);
+			//ASSERT_NEAR(btFabs(btScalar(i-z))-btScalar(j)-ssd.m_radiusB, distInfo.m_distance, abs_error);
+			//btVector3 computedA = distInfo.m_pointOnB+distInfo.m_distance*distInfo.m_normalBtoA;
+			//ASSERT_NEAR(computedA.x(),distInfo.m_pointOnA.x(),abs_error);
+			//ASSERT_NEAR(computedA.y(),distInfo.m_pointOnA.y(),abs_error);
+			//ASSERT_NEAR(computedA.z(),distInfo.m_pointOnA.z(),abs_error);
+		}
+
+
+#if 0
+		btCollisionDescription colDesc;
+		colDesc.m_objA = min0;
+		colDesc.m_objB = min1;
+		colDesc.m_localSupportFuncA = &btBulletShapeSupportFunc;
+		colDesc.m_localSupportFuncB = &btBulletShapeSupportFunc;
+		colDesc.m_localOriginFuncA = &btBulletShapeCenterFunc;
+		colDesc.m_localOriginFuncB = &btBulletShapeCenterFunc;
+
+		colDesc.m_transformA = body0Wrap->getWorldTransform();
+		colDesc.m_transformB = body1Wrap->getWorldTransform();
+		colDesc.m_marginA = body0Wrap->getCollisionShape()->getMargin();
+		colDesc.m_marginB = body1Wrap->getCollisionShape()->getMargin();
+		btDistanceInfo distInfo;
+		//int	result = btComputeGjkEpaPenetration(colDesc, &distInfo);
+		//int	result = btComputeGjkEpaPenetration2(colDesc, &distInfo);
+		int	result = btComputeMprPenetration(colDesc, &distInfo);
+		
+		if (result==0)
+		{
+			resultOut->addContactPoint(distInfo.m_normalBtoA,distInfo.m_pointOnB,distInfo.m_distance);
+		}
+
+		//bool res = b3MprPenetration(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,convexData,collidable2,cpuVertices,sepAxis,hasSepAxis,depthOut,dirOut,posOut);
+	
+		/*btCollisionDescription colDesc;
+		btDistanceInfo distInfo;
+		int	btComputeGjkEpaPenetration(min0, min1, &colDesc, &distInfo);
+		*/
+#endif
+	} else
+	{
+
+		btGjkPairDetector	gjkPairDetector(min0,min1,&vs,&epa);//m_simplexSolver,m_pdSolver);
+		//TODO: if (dispatchInfo.m_useContinuous)
+		gjkPairDetector.setMinkowskiA(min0);
+		gjkPairDetector.setMinkowskiB(min1);
+
+		{
+			//if (dispatchInfo.m_convexMaxDistanceUseCPT)
+			//{
+			//	input.m_maximumDistanceSquared = min0->getMargin() + min1->getMargin() + m_manifoldPtr->getContactProcessingThreshold();
+			//} else
+			//{
+			input.m_maximumDistanceSquared = min0->getMargin() + min1->getMargin() + m_manifoldPtr->getContactBreakingThreshold();
+	//		}
+
+			input.m_maximumDistanceSquared*= input.m_maximumDistanceSquared;
+		}
+
+		input.m_transformA = body0Wrap->getWorldTransform();
+		input.m_transformB = body1Wrap->getWorldTransform();
+
+
+		gjkPairDetector.getClosestPoints(input,*resultOut,dispatchInfo.m_debugDraw);
+	}
+	if (m_ownManifold)
+	{
+		resultOut->refreshContactPoints();
+	}
+
+}
+
+
+btScalar	btConvexConvexMprAlgorithm::calculateTimeOfImpact(btCollisionObject* col0,btCollisionObject* col1,const btDispatcherInfo& dispatchInfo,btManifoldResult* resultOut)
+{
+	(void)resultOut;
+	(void)dispatchInfo;
+	btAssert(0);
+	return 0;
+}
+
diff --git a/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.h b/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.h
new file mode 100644
index 000000000..6edbdccae
--- /dev/null
+++ b/Demos/VoronoiFractureDemo/btConvexConvexMprAlgorithm.h
@@ -0,0 +1,87 @@
+/*
+Bullet Continuous Collision Detection and Physics Library
+Copyright (c) 2003-2006 Erwin Coumans  http://continuousphysics.com/Bullet/
+
+This software is provided 'as-is', without any express or implied warranty.
+In no event will the authors be held liable for any damages arising from the use of this software.
+Permission is granted to anyone to use this software for any purpose, 
+including commercial applications, and to alter it and redistribute it freely, 
+subject to the following restrictions:
+
+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.
+2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+3. This notice may not be removed or altered from any source distribution.
+*/
+
+#ifndef BT_CONVEX_CONVEX_MPR_ALGORITHM_H
+#define BT_CONVEX_CONVEX_MPR_ALGORITHM_H
+
+#include "BulletCollision/CollisionDispatch/btActivatingCollisionAlgorithm.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
+#include "BulletCollision/NarrowPhaseCollision/btPersistentManifold.h"
+#include "BulletCollision/BroadphaseCollision/btBroadphaseProxy.h"
+#include "BulletCollision/NarrowPhaseCollision/btVoronoiSimplexSolver.h"
+#include "BulletCollision/CollisionDispatch/btCollisionCreateFunc.h"
+#include "BulletCollision/CollisionDispatch/btCollisionDispatcher.h"
+#include "LinearMath/btTransformUtil.h" //for btConvexSeparatingDistanceUtil
+
+class btConvexPenetrationDepthSolver;
+
+///Enabling USE_SEPDISTANCE_UTIL2 requires 100% reliable distance computation. However, when using large size ratios GJK can be imprecise
+///so the distance is not conservative. In that case, enabling this USE_SEPDISTANCE_UTIL2 would result in failing/missing collisions.
+///Either improve GJK for large size ratios (testing a 100 units versus a 0.1 unit object) or only enable the util
+///for certain pairs that have a small size ratio
+
+///The convexConvexAlgorithm collision algorithm implements time of impact, convex closest points and penetration depth calculations between two convex objects.
+///Multiple contact points are calculated by perturbing the orientation of the smallest object orthogonal to the separating normal.
+///This idea was described by Gino van den Bergen in this forum topic http://www.bulletphysics.com/Bullet/phpBB3/viewtopic.php?f=4&t=288&p=888#p888
+class btConvexConvexMprAlgorithm : public btActivatingCollisionAlgorithm
+{
+
+	
+	bool	m_ownManifold;
+	btPersistentManifold*	m_manifoldPtr;
+	///cache separating vector to speedup collision detection
+	
+public:
+
+	btConvexConvexMprAlgorithm(btPersistentManifold* mf,const btCollisionAlgorithmConstructionInfo& ci,const btCollisionObjectWrapper* body0Wrap,const btCollisionObjectWrapper* body1Wrap);
+
+	virtual ~btConvexConvexMprAlgorithm();
+
+	virtual void processCollision (const btCollisionObjectWrapper* body0Wrap,const btCollisionObjectWrapper* body1Wrap,const btDispatcherInfo& dispatchInfo,btManifoldResult* resultOut);
+
+	virtual btScalar calculateTimeOfImpact(btCollisionObject* body0,btCollisionObject* body1,const btDispatcherInfo& dispatchInfo,btManifoldResult* resultOut);
+
+	virtual	void	getAllContactManifolds(btManifoldArray&	manifoldArray)
+	{
+		///should we use m_ownManifold to avoid adding duplicates?
+		if (m_manifoldPtr && m_ownManifold)
+			manifoldArray.push_back(m_manifoldPtr);
+	}
+
+const btPersistentManifold*	getManifold()
+	{
+		return m_manifoldPtr;
+	}
+
+	struct CreateFunc :public 	btCollisionAlgorithmCreateFunc
+	{
+
+		
+
+		CreateFunc();
+		
+		virtual ~CreateFunc();
+
+		virtual	btCollisionAlgorithm* CreateCollisionAlgorithm(btCollisionAlgorithmConstructionInfo& ci, const btCollisionObjectWrapper* body0Wrap,const btCollisionObjectWrapper* body1Wrap)
+		{
+			void* mem = ci.m_dispatcher1->allocateCollisionAlgorithm(sizeof(btConvexConvexMprAlgorithm));
+			return new(mem) btConvexConvexMprAlgorithm(ci.m_manifold,ci,body0Wrap,body1Wrap);
+		}
+	};
+
+
+};
+
+#endif //BT_CONVEX_CONVEX_MPR_ALGORITHM_H
diff --git a/Demos/VoronoiFractureDemo/main.cpp b/Demos/VoronoiFractureDemo/main.cpp
index e1f7a9d76..9da8fe59e 100644
--- a/Demos/VoronoiFractureDemo/main.cpp
+++ b/Demos/VoronoiFractureDemo/main.cpp
@@ -25,7 +25,7 @@ subject to the following restrictions:
 int main(int argc,char** argv)
 {
 
-
+	
 	VoronoiFractureDemo ccdDemo;
 	ccdDemo.initPhysics();
 
diff --git a/Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem_add_plasticity_force.h b/Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem_add_plasticity_force.h
index 78ad78e31..12097a692 100644
--- a/Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem_add_plasticity_force.h
+++ b/Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem_add_plasticity_force.h
@@ -8,7 +8,7 @@
 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
 //
 #include <OpenTissue/configuration.h>
-
+#include <algorithm> //std::min
 namespace OpenTissue
 {
   namespace fem
@@ -22,7 +22,7 @@ namespace OpenTissue
 		{
 
 			        using std::min;
-        using std::sqrt;
+					using std::sqrt;
 
         typedef typename tetrahedron_type::vector3_type     vector3_type;