]> gitweb.factorcode.org Git - factor.git/blobdiff - extra/math/matrices/simd/simd.factor
factor: trim using lists
[factor.git] / extra / math / matrices / simd / simd.factor
index d65d1c4103bb326129aedde0c0dce07db7c82d1f..4ee38a3d003caf8b8fba343d460e55658fe8f2a2 100644 (file)
@@ -1,10 +1,13 @@
-! (c)Joe Groff bsd license
-USING: accessors classes.struct fry generalizations kernel locals
-math math.combinatorics math.functions math.matrices.simd math.vectors
-math.vectors.simd sequences sequences.private specialized-arrays
+! Copyright (C) 2009, 2011 Joe Groff.
+! See http://factorcode.org/license.txt for BSD license.
+USING: accessors classes.struct generalizations kernel
+math math.functions math.matrices.simd math.vectors
+math.vectors.simd math.quaternions sequences
+sequences.generalizations sequences.private specialized-arrays
 typed ;
+FROM: sequences.private => nth-unsafe ;
+FROM: math.quaternions.private => (q*sign) ;
 QUALIFIED-WITH: alien.c-types c
-SIMD: c:float
 SPECIALIZED-ARRAY: float-4
 IN: math.matrices.simd
 
@@ -24,19 +27,16 @@ M: matrix4 new-sequence 2drop matrix4 (struct) ; inline
 
 :: set-columns ( c1 c2 c3 c4 c -- c )
     c columns>> :> columns
-    c1 columns set-first
-    c2 columns set-second
-    c3 columns set-third
-    c4 columns set-fourth
+    c1 c2 c3 c4 columns 4 set-firstn-unsafe
     c ; inline
 
-: make-matrix4 ( quot: ( -- c1 c2 c3 c4 ) -- c )
+: make-matrix4 ( ..a quot: ( ..a -- ..b c1 c2 c3 c4 ) -- ..b c )
     matrix4 (struct) swap dip set-columns ; inline
 
 :: 2map-columns ( a b quot -- c )
     [
-        a columns :> a4 :> a3 :> a2 :> a1
-        b columns :> b4 :> b3 :> b2 :> b1
+        a columns :> ( a1 a2 a3 a4 )
+        b columns :> ( b1 b2 b3 b4 )
 
         a1 b1 quot call
         a2 b2 quot call
@@ -44,9 +44,9 @@ M: matrix4 new-sequence 2drop matrix4 (struct) ; inline
         a4 b4 quot call
     ] make-matrix4 ; inline
 
-: map-columns ( a quot -- c )
+: map-columns ( ... a quot: ( ... col -- ... newcol ) -- ... c )
     '[ columns _ 4 napply ] make-matrix4 ; inline
-    
+
 PRIVATE>
 
 TYPED: m4+ ( a: matrix4 b: matrix4 -- c: matrix4 ) [ v+ ] 2map-columns ;
@@ -61,20 +61,20 @@ TYPED: n/m4 ( a: float b: matrix4 -- c: matrix4 ) [ n/v ] with map-columns ;
 
 TYPED:: m4. ( a: matrix4 b: matrix4 -- c: matrix4 )
     [
-        a columns :> a4 :> a3 :> a2 :> a1
-        b columns :> b4 :> b3 :> b2 :> b1
+        a columns :> ( a1 a2 a3 a4 )
+        b columns :> ( b1 b2 b3 b4 )
 
         b1 first  a1 n*v :> c1a
         b2 first  a1 n*v :> c2a
         b3 first  a1 n*v :> c3a
         b4 first  a1 n*v :> c4a
 
-        b1 second a2 n*v c1a v+ :> c1b 
+        b1 second a2 n*v c1a v+ :> c1b
         b2 second a2 n*v c2a v+ :> c2b
         b3 second a2 n*v c3a v+ :> c3b
         b4 second a2 n*v c4a v+ :> c4b
 
-        b1 third  a3 n*v c1b v+ :> c1c 
+        b1 third  a3 n*v c1b v+ :> c1c
         b2 third  a3 n*v c2b v+ :> c2c
         b3 third  a3 n*v c3b v+ :> c3c
         b4 third  a3 n*v c4b v+ :> c4c
@@ -86,15 +86,15 @@ TYPED:: m4. ( a: matrix4 b: matrix4 -- c: matrix4 )
     ] make-matrix4 ;
 
 TYPED:: m4.v ( m: matrix4 v: float-4 -- v': float-4 )
-    m columns :> m4 :> m3 :> m2 :> m1
-    
+    m columns :> ( m1 m2 m3 m4 )
+
     v first  m1 n*v
     v second m2 n*v v+
     v third  m3 n*v v+
     v fourth m4 n*v v+ ;
 
 TYPED:: v.m4 ( v: float-4 m: matrix4 -- c: float-4 )
-    m columns [ v v. ] 4 napply float-4-boa ;
+    m columns [ v vdot ] 4 napply float-4-boa ;
 
 CONSTANT: identity-matrix4
     S{ matrix4 f
@@ -145,19 +145,31 @@ TYPED: transpose-matrix4 ( matrix: matrix4 -- matrix: matrix4 )
 TYPED: translation-matrix4 ( offset: float-4 -- matrix: matrix4 )
     [
         linear>homogeneous
-        [ 
+        [
             float-4{ 1.0 0.0 0.0 0.0 }
             float-4{ 0.0 1.0 0.0 0.0 }
             float-4{ 0.0 0.0 1.0 0.0 }
         ] dip
     ] make-matrix4 ;
 
+:: (rotation-matrix4) ( diagonal triangle-hi triangle-lo -- matrix )
+    matrix4 (struct) :> triangle-m
+    diagonal scale-matrix4 :> diagonal-m
+
+    triangle-hi { 3 2 1 3 } vshuffle
+    triangle-hi { 3 3 0 3 } vshuffle triangle-lo { 2 3 3 3 } vshuffle vbitor
+                                     triangle-lo { 1 0 3 3 } vshuffle
+    float-4 new
+
+    triangle-m set-columns drop
+
+    diagonal-m triangle-m m4+ ; inline
+
 TYPED:: rotation-matrix4 ( axis: float-4 theta: float -- matrix: matrix4 )
     !   x*x + c*(1.0 - x*x)   x*y*(1.0 - c) + s*z   x*z*(1.0 - c) - s*y   0
     !   x*y*(1.0 - c) - s*z   y*y + c*(1.0 - y*y)   y*z*(1.0 - c) + s*x   0
     !   x*z*(1.0 - c) + s*y   y*z*(1.0 - c) - s*x   z*z + c*(1.0 - z*z)   0
     !   0                     0                     0                     1
-    matrix4 (struct) :> triangle-m
     theta cos :> c
     theta sin :> s
 
@@ -177,23 +189,14 @@ TYPED:: rotation-matrix4 ( axis: float-4 theta: float -- matrix: matrix4 )
     triangle-a triangle-b v+ :> triangle-lo
     triangle-a triangle-b v- :> triangle-hi
 
-    diagonal scale-matrix4 :> diagonal-m
-
-    triangle-hi { 3 2 1 3 } vshuffle
-    triangle-hi { 3 3 0 3 } vshuffle triangle-lo { 2 3 3 3 } vshuffle v+
-    triangle-lo { 1 0 3 3 } vshuffle
-    float-4 new
-
-    triangle-m set-columns drop
-
-    diagonal-m triangle-m m4+ ;
+    diagonal triangle-hi triangle-lo (rotation-matrix4) ;
 
 TYPED:: frustum-matrix4 ( xy: float-4 near: float far: float -- matrix: matrix4 )
     [
         near near near far + 2 near far * * float-4-boa ! num
         float-4{ t t f f } xy near far - float-4-with v? ! denom
         v/ :> fov
-        
+
         float-4{ 0.0 -1.0 0.0 0.0 } :> negone
 
         fov vmerge-diagonal
@@ -201,3 +204,30 @@ TYPED:: frustum-matrix4 ( xy: float-4 near: float far: float -- matrix: matrix4
         [ negone (vmerge) ] bi*
     ] make-matrix4 ;
 
+! interface with quaternions
+M: float-4 (q*sign)
+    float-4{ -0.0  0.0  0.0  0.0 } vbitxor ; inline
+M: float-4 qconjugate
+    float-4{  0.0 -0.0 -0.0 -0.0 } vbitxor ; inline
+
+: euler4 ( phi theta psi -- q )
+    float-4{ 0 0 0 0 } euler-like ; inline
+
+TYPED:: q>matrix4 ( q: float-4 -- matrix: matrix4 )
+    !   a*a + b*b - c*c - d*d  2*b*c - 2*a*d          2*b*d + 2*a*c          0
+    !   2*b*c + 2*a*d          a*a - b*b + c*c - d*d  2*c*d - 2*a*b          0
+    !   2*b*d - 2*a*c          2*c*d + 2*a*b          a*a - b*b - c*c + d*d  0
+    !   0                      0                      0                      1
+    q { 2 1 1 3 } vshuffle  q { 3 3 2 3 } vshuffle  v*  :> triangle-a
+    q { 0 0 0 3 } vshuffle  q { 1 2 3 3 } vshuffle  v*  :> triangle-b
+
+    triangle-a float-4{ 2.0 2.0 2.0 0.0 } v*  triangle-b float-4{ -2.0 2.0 -2.0 0.0 } v*
+    [ v- ] [ v+ ] 2bi :> ( triangle-hi triangle-lo )
+
+    q q v* first4 {
+        [ [ + ] [ - ] [ - ] tri* ]
+        [ [ - ] [ + ] [ - ] tri* ]
+        [ [ - ] [ - ] [ + ] tri* ]
+    } 4 ncleave 1.0 float-4-boa :> diagonal
+
+    diagonal triangle-hi triangle-lo (rotation-matrix4) ;