commit 9f4c44fe7e123060fec7553cffc56de6e68f37e8
Author: Leszek Koltunski <leszek@koltunski.pl>
Date:   Fri Mar 6 22:27:38 2020 +0000

    Progress with the Pyraminx - computing all legal quaternions!

diff --git a/src/main/java/org/distorted/object/Cubit.java b/src/main/java/org/distorted/object/Cubit.java
index 6d853d17..f6570ae4 100644
--- a/src/main/java/org/distorted/object/Cubit.java
+++ b/src/main/java/org/distorted/object/Cubit.java
@@ -65,6 +65,8 @@ class Cubit
 
   private void normalizeScrambleQuat(Static4D quat)
     {
+    final float MAX_ERROR = 0.0001f;
+
     float x = quat.get0();
     float y = quat.get1();
     float z = quat.get2();
@@ -74,13 +76,13 @@ class Cubit
     for(float legal: mParent.LEGAL_QUATS)
       {
       diff = x-legal;
-      if( diff*diff<0.01f ) x = legal;
+      if( diff*diff<MAX_ERROR ) x = legal;
       diff = y-legal;
-      if( diff*diff<0.01f ) y = legal;
+      if( diff*diff<MAX_ERROR ) y = legal;
       diff = z-legal;
-      if( diff*diff<0.01f ) z = legal;
+      if( diff*diff<MAX_ERROR ) z = legal;
       diff = w-legal;
-      if( diff*diff<0.01f ) w = legal;
+      if( diff*diff<MAX_ERROR ) w = legal;
       }
 
     if( w<0 )
diff --git a/src/main/java/org/distorted/object/RubikCube.java b/src/main/java/org/distorted/object/RubikCube.java
index 654e141f..712c94a3 100644
--- a/src/main/java/org/distorted/object/RubikCube.java
+++ b/src/main/java/org/distorted/object/RubikCube.java
@@ -62,6 +62,8 @@ class RubikCube extends RubikObject
   // multiplying them and eventually you'll find all (24) legal rotations.
   // 3) linear scan through those shows that the only floats in those 24 quats are those 7 given
   // below.
+  //
+  // Example program in C, res/raw/compute_quats.c , is included.
   private static final float[] LEGALQUATS = new float[]
          {
            0.0f ,
diff --git a/src/main/java/org/distorted/object/RubikPyraminx.java b/src/main/java/org/distorted/object/RubikPyraminx.java
index 094fa95c..e0d01b76 100644
--- a/src/main/java/org/distorted/object/RubikPyraminx.java
+++ b/src/main/java/org/distorted/object/RubikPyraminx.java
@@ -40,12 +40,15 @@ import org.distorted.library.type.Static4D;
 
 public class RubikPyraminx extends RubikObject
 {
+  private static final float SQ2 = (float)Math.sqrt(2);
+  private static final float SQ3 = (float)Math.sqrt(3);
+
   private static final Static3D[] AXIS = new Static3D[]
          {
-           new Static3D(                     0,        1,                       0 ),
-           new Static3D( (float)Math.sqrt(6)/3,  -1.0f/3,  -(float)Math.sqrt(3)/3 ),
-           new Static3D(-(float)Math.sqrt(6)/3,  -1.0f/3,  -(float)Math.sqrt(3)/3 ),
-           new Static3D(                     0,  -1.0f/3, 2*(float)Math.sqrt(2)/3 )
+           new Static3D(         0,        1,       0 ),
+           new Static3D( SQ2*SQ3/3,  -1.0f/3,  -SQ2/3 ),
+           new Static3D(-SQ2*SQ3/3,  -1.0f/3,  -SQ2/3 ),
+           new Static3D(         0,  -1.0f/3, 2*SQ2/3 )
          };
 
   private static final int[] FACE_COLORS = new int[]
@@ -54,9 +57,11 @@ public class RubikPyraminx extends RubikObject
            0xff0000ff, 0xffff0000   // AXIS[2]right (BLUE  ) AXIS[3]right (RED   )
          };
 
+  // computed with res/raw/compute_quats.c
   private static final float[] LEGALQUATS = new float[]
          {
-         // TODO;
+           0.0f, 1.0f, -1.0f, 0.5f, -0.5f, SQ2/2, -SQ2/2, SQ3/2, -SQ3/2,
+           SQ3/3, -SQ3/3, SQ3/6, -SQ3/6, SQ2*SQ3/3, -SQ2*SQ3/3, SQ2*SQ3/6, -SQ2*SQ3/6
          };
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -106,7 +111,6 @@ public class RubikPyraminx extends RubikObject
 
   MeshBase createCubitMesh(int vertices)
     {
-    final float SQ3 = (float)Math.sqrt(3);
     final float angleFaces = (float)((180/Math.PI)*(2*Math.asin(SQ3/3))); // angle between two faces of a tetrahedron
     final int MESHES=4;
 
diff --git a/src/main/res/raw/compute_quats.c b/src/main/res/raw/compute_quats.c
new file mode 100644
index 00000000..0936c6b8
--- /dev/null
+++ b/src/main/res/raw/compute_quats.c
@@ -0,0 +1,141 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#define PYRAMIX
+
+#define SQ2 1.41421356237f
+#define SQ3 1.73205080757f
+#define PI  3.14159265358f
+#define NUM_QUATS  100
+
+#ifdef PYRAMIX 
+#define NUM_AXIS    4
+#define BASIC_ANGLE 3
+
+float axis[NUM_AXIS][3] ={ {         0,       1,       0 } ,
+                           { SQ2*SQ3/3, -1.0f/3,  -SQ2/3 } ,
+                           {-SQ2*SQ3/3, -1.0f/3,  -SQ2/3 } ,
+                           {         0, -1.0f/3, 2*SQ2/3 } };
+#endif
+
+#ifdef CUBE
+#define NUM_AXIS 3
+#define BASIC_ANGLE 4
+float axis[NUM_AXIS][3] = { { 1,0,0 }, {0,1,0}, {0,0,1} };
+#endif
+
+float* quats;
+float* table;
+int inserted=0;
+
+///////////////////////////////////////////////////////////////////
+// q1*q2
+
+void multiply_quats( float* q1, float* q2, float* output)
+  {
+  output[0] = q2[3]*q1[0] - q2[2]*q1[1] + q2[1]*q1[2] + q2[0]*q1[3];
+  output[1] = q2[3]*q1[1] + q2[2]*q1[0] + q2[1]*q1[3] - q2[0]*q1[2];
+  output[2] = q2[3]*q1[2] + q2[2]*q1[3] - q2[1]*q1[0] + q2[0]*q1[1];
+  output[3] = q2[3]*q1[3] - q2[2]*q1[2] - q2[1]*q1[1] - q2[0]*q1[0];
+  }
+
+///////////////////////////////////////////////////////////////////
+// sin(A/2)*x, sin(A/2)*y, sin(A/2)*z, cos(A/2)
+
+void create_quat(float* axis, float angle, float* output)
+  {
+  float cosAngle = cos(angle/2);
+  float sinAngle = sin(angle/2);
+
+  output[0] = sinAngle*axis[0];
+  output[1] = sinAngle*axis[1];
+  output[2] = sinAngle*axis[2];
+  output[3] = cosAngle;
+  }
+
+///////////////////////////////////////////////////////////////////
+// double cover, so q == -q
+
+int is_the_same(float* q1, float* q2)
+  {
+  const float MAX = 0.01f;
+
+  float d0 = q1[0]-q2[0];
+  float d1 = q1[1]-q2[1];
+  float d2 = q1[2]-q2[2];
+  float d3 = q1[3]-q2[3];
+
+  if( d0<MAX && d0>-MAX &&
+      d1<MAX && d1>-MAX &&
+      d2<MAX && d2>-MAX &&
+      d3<MAX && d3>-MAX  ) return 1;
+
+  d0 = q1[0]+q2[0];
+  d1 = q1[1]+q2[1];
+  d2 = q1[2]+q2[2];
+  d3 = q1[3]+q2[3];
+
+  if( d0<MAX && d0>-MAX &&
+      d1<MAX && d1>-MAX &&
+      d2<MAX && d2>-MAX &&
+      d3<MAX && d3>-MAX  ) return 1;
+
+  return 0;
+  }
+
+///////////////////////////////////////////////////////////////////
+
+void insert(float* quat, float* to)
+  {
+  for(int i=0; i<inserted; i++)
+    {
+    if( is_the_same(quat,to+4*i)==1 ) return;
+    }
+
+  to[4*inserted+0] = quat[0];
+  to[4*inserted+1] = quat[1];
+  to[4*inserted+2] = quat[2];
+  to[4*inserted+3] = quat[3];
+
+  inserted++;
+  }
+
+///////////////////////////////////////////////////////////////////
+
+int main(int argc, char** argv)
+  {
+  float tmp[4];
+  int num = 1+NUM_AXIS*(BASIC_ANGLE-1);
+
+  quats = (float*) malloc(4*sizeof(float)*num      );
+  table = (float*) malloc(4*sizeof(float)*NUM_QUATS);
+
+  tmp[0] = 0.0f; tmp[1] = 0.0f; tmp[2] = 0.0f; tmp[3] = 1.0f;
+  insert(tmp,quats);
+
+  for(int angle=1; angle<BASIC_ANGLE; angle++)
+    for( int ax=0; ax<NUM_AXIS; ax++)
+      {
+      create_quat(axis[ax], 2*PI*angle/BASIC_ANGLE, tmp);
+      insert(tmp,quats); 
+      }
+
+  inserted=0;
+
+  for(int i=0; i<num; i++)
+    for(int j=0; j<num; j++)
+      {
+      multiply_quats( quats+4*i, quats+4*j, tmp);
+      insert(tmp,table);
+      }
+
+  printf("inserted: %d\n", inserted);
+
+  for(int i=0; i<inserted; i++)
+    {
+    printf( "%d %7.4f %7.4f %7.4f %7.4f\n", i, table[4*i], table[4*i+1], table[4*i+2], table[4*i+3] );
+    }
+
+return 0;
+}