Merge remote-tracking branch 'origin/master' into modern-harvey3

.travis.yml

@@ -35,6 +35,15 @@ addons:
       - gtk2-engines-pixbuf
 before_install:
   - uname -s
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export HOMEBREW_NO_AUTO_UPDATE=1 ; fi # Don't let homebrew upgrade itself
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then rm -rf ~/.gnupg/; fi  # https://github.com/rvm/rvm/issues/3110
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://rvm.io/mpapis.asc | gpg --import - ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://rvm.io/pkuczynski.asc | gpg --import - ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://get.rvm.io | bash -s stable; fi  # https://github.com/travis-ci/travis-ci/issues/6307
+  #- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then rvm reload ; fi # for homebrew to have 2.6.3, which takes too long. instead we just use HOMEBREW_NO_AUTO_UPDATE=1
+  #- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then rvm install ruby-2.6.3 ; fi # for homebrew
+  #- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then rvm use 2.6 ; fi # for homebrew
+
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then ./build.sh deps-macosx ; else ./build.sh deps-apt-get ; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew ls --versions snappy > /dev/null || brew install snappy; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew ls --versions cmake > /dev/null || brew install cmake; fi
@@ -47,10 +56,6 @@ before_install:
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew services start redis; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew services start postgresql; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew services start memcached; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then rm -rf ~/.gnupg/; fi  # https://github.com/rvm/rvm/issues/3110
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://rvm.io/mpapis.asc | gpg --import - ; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://rvm.io/pkuczynski.asc | gpg --import - ; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then curl -sSL https://get.rvm.io | bash -s stable; fi  # https://github.com/travis-ci/travis-ci/issues/6307
   - if [[ "$TRAVIS_OS_NAME" != "windows" ]]; then
       wget https://github.com/vmt/udis86/archive/v1.7.2.tar.gz && tar xzvf v1.7.2.tar.gz &&
       ( cd udis86-1.7.2/ && ./autogen.sh && ./configure --enable-shared=yes && make && sudo make install ) &&
@@ -69,7 +74,7 @@ script:
   - export CI_BRANCH="${TRAVIS_PULL_REQUEST_BRANCH:-$TRAVIS_BRANCH}"
   - echo "CI_BRANCH=${CI_BRANCH}"
   - DEBUG=1 ./build.sh net-bootstrap < /dev/null
-  - "./factor -e='USING: memory vocabs.hierarchy ; \"zealot\" load save'"
+  - "./factor -e='USING: memory vocabs.hierarchy tools.test namespaces ; \"zealot\" load f long-unit-tests-enabled? set-global save'"
   - './factor -run=zealot.cli-changed-vocabs'
   - './factor -run=tools.test `./factor -run=zealot.cli-changed-vocabs | paste -s -d " " -`'
   - './factor -run=zealot.help-lint `./factor -run=zealot.cli-changed-vocabs | paste -s -d " " -`'

basis/compression/run-length/run-length.factor

@@ -13,7 +13,7 @@ IN: compression.run-length
 
 :: run-length-uncompress-bitmap4 ( byte-array m n -- byte-array' )
     byte-array <sequence-parser> :> sp
-    m  1 + n zero-matrix :> matrix
+    m  1 + n <zero-matrix> :> matrix
     n 4 mod n + :> stride
     0 :> i!
     0 :> j!
@@ -45,7 +45,7 @@ IN: compression.run-length
 
 :: run-length-uncompress-bitmap8 ( byte-array m n -- byte-array' )
     byte-array <sequence-parser> :> sp
-    m  1 + n zero-matrix :> matrix
+    m  1 + n <zero-matrix> :> matrix
     n 4 mod n + :> stride
     0 :> i!
     0 :> j!

basis/io/encodings/string/string.factor

@@ -13,7 +13,7 @@ IN: io.encodings.string
         ] [
             byte-array encoding <byte-reader> :> reader
             byte-array length encoding guess-decoded-length <sbuf> :> buf
-            [ reader stream-read1 dup ] [ buf push ] while drop
+            [ reader stream-read1 ] [ buf push ] while*
             buf "" like
         ] if
     ] if ; inline

basis/json/reader/reader.factor

@@ -126,7 +126,7 @@ DEFER: (read-json-string)
     } case ;
 
 : json-read-input ( stream -- objects )
-    V{ } clone over '[ _ stream-read1 dup ] [ scan ] while drop nip ;
+    V{ } clone over '[ _ stream-read1 ] [ scan ] while* nip ;
 
 ! If there are no json objects, return an empty hashtable
 ! This happens for empty files.

basis/math/matrices/authors.txt

@@ -1 +1,4 @@
 Slava Pestov
+Joe Groff
+Doug Coleman
+Cat Stevens

basis/math/matrices/matrices.factor

@@ -1,197 +1,159 @@
-! Copyright (C) 2005, 2010 Slava Pestov, Joe Groff.
+! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: accessors arrays columns kernel locals math math.bits
+USING: accessors arrays classes.singleton columns combinators
-math.functions math.order math.vectors sequences
+combinators.short-circuit combinators.smart formatting fry
-sequences.private fry math.statistics grouping
+grouping kernel locals math math.bits math.functions math.order
-combinators.short-circuit math.ranges combinators.smart ;
+math.private math.ranges math.statistics math.vectors
+math.vectors.private sequences sequences.deep sequences.private
+slots.private summary ;
 IN: math.matrices
 
-! Matrices
+! defined here because of issue #1943
-: make-matrix ( m n quot -- matrix )
+DEFER: regular-matrix?
-    '[ _ _ replicate ] replicate ; inline
+: regular-matrix? ( object -- ? )
+    [ t ] [
+        dup first-unsafe length
+        '[ length _ = ] all?
+    ] if-empty ;
 
+! the MRO (class linearization) is performed in the order the predicates appear here
+! except that null-matrix is last (but it is relied upon by zero-matrix)
+! in other words:
+! sequence > matrix > zero-matrix > square-matrix > zero-square-matrix > null-matrix
+
+! Factor bug that's hard to repro: using `bi and` in these predicates
+! instead of 1&& will cause spirious no-method and bounds-error errors in <square-cols>
+! and the tests/docs for no apparent reason
+PREDICATE: matrix < sequence
+    { [ [ sequence? ] all? ] [ regular-matrix? ] } 1&& ;
+
+PREDICATE: irregular-matrix < sequence
+    { [ [ sequence? ] all? ] [ regular-matrix? not ] } 1&& ;
+
+DEFER: dimension
+! can't define dim using this predicate for this reason,
+! unless we are going to write two versions of dim, one of which is generic
+PREDICATE: square-matrix < matrix
+    dimension first2-unsafe = ;
+
+PREDICATE: null-matrix < matrix
+    flatten empty? ;
+
+PREDICATE: zero-matrix < matrix
+    dup null-matrix? [ drop f ] [ flatten [ zero? ] all? ] if ;
+
+PREDICATE: zero-square-matrix < square-matrix
+    { [ zero-matrix? ] [ square-matrix? ] } 1&& ;
+
+! Benign matrix constructors
 : <matrix> ( m n element -- matrix )
     '[ _ _ <array> ] replicate ; inline
 
-: zero-matrix ( m n -- matrix )
+: <matrix-by> ( m n quot: ( ... -- elt ) -- matrix )
+    '[ _ _ replicate ] replicate ; inline
+
+: <matrix-by-indices> ( ... m n quot: ( ... m' n' -- ... elt ) -- ... matrix )
+    [ [ <iota> ] bi@ ] dip cartesian-map ; inline
+
+: <zero-matrix> ( m n -- matrix )
     0 <matrix> ; inline
 
-: diagonal-matrix ( diagonal-seq -- matrix )
+: <zero-square-matrix> ( n -- matrix )
-    dup length dup zero-matrix
+    dup <zero-matrix> ; inline
-    [ '[ dup _ nth set-nth ] each-index ] keep ; inline
 
-: identity-matrix ( n -- matrix )
+<PRIVATE
-    1 <repetition> diagonal-matrix ; inline
+: (nth-from-end) ( n seq -- n )
+    length 1 - swap - ; inline flushable
 
-: eye ( m n k -- matrix )
+: nth-end ( n seq -- elt )
-    [ [ <iota> ] bi@ ] dip neg '[ _ + = 1 0 ? ] cartesian-map ;
+    [ (nth-from-end) ] keep nth ; inline flushable
 
-: hilbert-matrix ( m n -- matrix )
+: nth-end-unsafe ( n seq -- elt )
-    [ <iota> ] bi@ [ + 1 + recip ] cartesian-map ;
+    [ (nth-from-end) ] keep nth-unsafe ; inline flushable
 
-: toeplitz-matrix ( n -- matrix )
+: array-nth-end-unsafe ( n seq -- elt )
-    <iota> dup [ - abs 1 + ] cartesian-map ;
+    [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
 
-: hankel-matrix ( n -- matrix )
+: set-nth-end ( elt n seq -- )
-    [ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
+    [ (nth-from-end) ] keep set-nth ; inline
 
-: box-matrix ( r -- matrix )
+: set-nth-end-unsafe ( elt n seq -- )
-    2 * 1 + dup '[ _ 1 <array> ] replicate ;
+    [ (nth-from-end) ] keep set-nth-unsafe ; inline
+PRIVATE>
 
-: vandermonde-matrix ( u n -- matrix )
+! main-diagonal matrix
-    <iota> [ v^n ] with map reverse flip ;
+: <diagonal-matrix> ( diagonal-seq -- matrix )
+    [ length <zero-square-matrix> ] keep over
+    '[ dup _ nth set-nth-unsafe ] each-index ; inline
 
-:: rotation-matrix3 ( axis theta -- matrix )
+! could also be written slower as: <diagonal-matrix> [ reverse ] map
-    theta cos :> c
+: <anti-diagonal-matrix> ( diagonal-seq -- matrix )
-    theta sin :> s
+    [ length <zero-square-matrix> ] keep over
-    axis first3 :> ( x y z )
+    '[ dup _ nth set-nth-end-unsafe ] each-index ; inline
-    x sq 1.0 x sq - c * +     x y * 1.0 c - * z s * -   x z * 1.0 c - * y s * + 3array
-    x y * 1.0 c - * z s * +   y sq 1.0 y sq - c * +     y z * 1.0 c - * x s * - 3array
-    x z * 1.0 c - * y s * -   y z * 1.0 c - * x s * +   z sq 1.0 z sq - c * +   3array
-    3array ;
 
-:: rotation-matrix4 ( axis theta -- matrix )
+: <identity-matrix> ( n -- matrix )
-    theta cos :> c
+    1 <repetition> <diagonal-matrix> ; inline
-    theta sin :> s
-    axis first3 :> ( x y z )
-    x sq 1.0 x sq - c * +     x y * 1.0 c - * z s * -   x z * 1.0 c - * y s * +   0 4array
-    x y * 1.0 c - * z s * +   y sq 1.0 y sq - c * +     y z * 1.0 c - * x s * -   0 4array
-    x z * 1.0 c - * y s * -   y z * 1.0 c - * x s * +   z sq 1.0 z sq - c * +     0 4array
-    { 0.0 0.0 0.0 1.0 } 4array ;
 
-:: translation-matrix4 ( offset -- matrix )
+: <eye> ( m n k z -- matrix )
-    offset first3 :> ( x y z )
+    [ [ <iota> ] bi@ ] 2dip
-    {
+    '[ _ neg + = _ 0 ? ]
-        { 1.0 0.0 0.0 x   }
+    cartesian-map ; inline
-        { 0.0 1.0 0.0 y   }
-        { 0.0 0.0 1.0 z   }
-        { 0.0 0.0 0.0 1.0 }
-    } ;
 
-: >scale-factors ( number/sequence -- x y z )
+! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
-    dup number? [ dup dup ] [ first3 ] if ;
+:: <simple-eye> ( m n k -- matrix )
+    m n = k 0 = and
+    [ n <identity-matrix> ]
+    [ m n k 1 <eye> ] if ; inline
 
-:: scale-matrix3 ( factors -- matrix )
+: <coordinate-matrix> ( dim -- coordinates )
-    factors >scale-factors :> ( x y z )
+  first2 [ <iota> ] bi@ cartesian-product ; inline
-    {
-        { x   0.0 0.0 }
-        { 0.0 y   0.0 }
-        { 0.0 0.0 z   }
-    } ;
 
-:: scale-matrix4 ( factors -- matrix )
+ALIAS: <cartesian-indices> <coordinate-matrix>
-    factors >scale-factors :> ( x y z )
-    {
-        { x   0.0 0.0 0.0 }
-        { 0.0 y   0.0 0.0 }
-        { 0.0 0.0 z   0.0 }
-        { 0.0 0.0 0.0 1.0 }
-    } ;
 
-: ortho-matrix4 ( dim -- matrix )
+: <cartesian-square-indices> ( n -- matrix )
-    [ recip ] map scale-matrix4 ;
+    dup 2array <cartesian-indices> ; inline
 
-:: frustum-matrix4 ( xy-dim near far -- matrix )
+ALIAS: transpose flip
-    xy-dim first2 :> ( x y )
-    near x /f :> xf
-    near y /f :> yf
-    near far + near far - /f :> zf
-    2 near far * * near far - /f :> wf
 
-    {
+<PRIVATE
-        { xf  0.0  0.0 0.0 }
+: array-matrix? ( matrix -- ? )
-        { 0.0 yf   0.0 0.0 }
+    [ array? ]
-        { 0.0 0.0  zf  wf  }
+    [ [ array? ] all? ] bi and ; inline foldable flushable
-        { 0.0 0.0 -1.0 0.0 }
-    } ;
 
-:: skew-matrix4 ( theta -- matrix )
+: matrix-cols-iota ( matrix -- cols-iota )
-    theta tan :> zf
+  first-unsafe length <iota> ; inline
 
-    {
+: unshaped-cols-iota ( matrix -- cols-iota )
-        { 1.0 0.0 0.0 0.0 }
+  [ first-unsafe length 1 ] keep
-        { 0.0 1.0 0.0 0.0 }
+  [ length min ] (each) (each-integer) <iota> ; inline
-        { 0.0 zf  1.0 0.0 }
-        { 0.0 0.0 0.0 1.0 }
-    } ;
 
-! Matrix operations
+: generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
-: mneg ( m -- m ) [ vneg ] map ;
+    [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
 
-: n+m  ( n m -- m ) [ n+v ] with map ;
+: array-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
-: m+n  ( m n -- m ) [ v+n ] curry map ;
+    [ <reversed> [ array-nth-end-unsafe ] with map ] curry map ; inline
-: n-m  ( n m -- m ) [ n-v ] with map ;
+PRIVATE>
-: m-n  ( m n -- m ) [ v-n ] curry map ;
-: n*m ( n m -- m ) [ n*v ] with map ;
-: m*n ( m n -- m ) [ v*n ] curry map ;
-: n/m ( n m -- m ) [ n/v ] with map ;
-: m/n ( m n -- m ) [ v/n ] curry map ;
 
-: m+   ( m m -- m ) [ v+ ] 2map ;
+! much faster than [ reverse ] map flip [ reverse ] map
-: m-   ( m m -- m ) [ v- ] 2map ;
+: anti-transpose ( matrix -- newmatrix )
-: m*   ( m m -- m ) [ v* ] 2map ;
+    dup empty? [ ] [
-: m/   ( m m -- m ) [ v/ ] 2map ;
+        [ dup regular-matrix?
+            [ matrix-cols-iota ] [ unshaped-cols-iota ] if
+        ] keep
 
-: v.m ( v m -- v ) flip [ v. ] with map ;
+        dup array-matrix? [
-: m.v ( m v -- v ) [ v. ] curry map ;
+            array-anti-transpose-unsafe
-: m.  ( m m -- m ) flip [ swap m.v ] curry map ;
+        ] [
+            generic-anti-transpose-unsafe
+        ] if
+    ] if ;
 
-: m~  ( m m epsilon -- ? ) [ v~ ] curry 2all? ;
+ALIAS: anti-flip anti-transpose
 
-: mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
+: row ( n matrix -- row )
-: mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
-: mnorm ( m -- n ) dup mmax abs m/n ;
-: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
-: m-1norm ( m -- n ) flip m-infinity-norm ;
-: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
-
-: cross ( vec1 vec2 -- vec3 )
-    [ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ]
-    [ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline
-
-:: normal ( vec1 vec2 vec3 -- vec4 )
-    vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
-
-: proj ( v u -- w )
-    [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
-
-: perp ( v u -- w )
-    dupd proj v- ;
-
-: angle-between ( v u -- a )
-    [ normalize ] bi@ h. acos ;
-
-: (gram-schmidt) ( v seq -- newseq )
-    [ dupd proj v- ] each ;
-
-: gram-schmidt ( seq -- orthogonal )
-    V{ } clone [ over (gram-schmidt) suffix! ] reduce ;
-
-: norm-gram-schmidt ( seq -- orthonormal )
-    gram-schmidt [ normalize ] map ;
-
-ERROR: negative-power-matrix m n ;
-
-: (m^n) ( m n -- n )
-    make-bits over first length identity-matrix
-    [ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
-
-: m^n ( m n -- n )
-    dup 0 >= [ (m^n) ] [ negative-power-matrix ] if ;
-
-: stitch ( m -- m' )
-    [ ] [ [ append ] 2map ] map-reduce ;
-
-: kron ( m1 m2 -- m )
-    '[ [ _ n*m  ] map ] map stitch stitch ;
-
-: outer ( u v -- m )
-    [ n*v ] curry map ;
-
-: row ( n matrix -- col )
     nth ; inline
 
-: rows ( seq matrix -- cols )
+: rows ( seq matrix -- rows )
     '[ _ row ] map ; inline
 
 : col ( n matrix -- col )
@@ -200,77 +162,151 @@ ERROR: negative-power-matrix m n ;
 : cols ( seq matrix -- cols )
     '[ _ col ] map ; inline
 
-: set-index ( object pair matrix -- )
+:: >square-matrix ( m -- subset )
-    [ first2 swap ] dip nth set-nth ; inline
+    m dimension first2 :> ( x y ) {
+        { [ x y = ] [ m ] }
+        { [ x y < ] [ x <iota> m cols transpose ] }
+        { [ x y > ] [ y <iota> m rows ] }
+    } cond ;
 
-: set-indices ( object sequence matrix -- )
+GENERIC: <square-rows> ( desc -- matrix )
-    '[ _ set-index ] with each ; inline
+M: integer <square-rows>
-
+    <iota> <square-rows> ;
-: matrix-map ( matrix quot -- )
+M: sequence <square-rows>
-    '[ _ map ] map ; inline
-
-: column-map ( matrix quot -- seq )
-    [ [ first length <iota> ] keep ] dip '[ _ col @ ] map ; inline
-
-: cartesian-square-indices ( n -- matrix )
-    <iota> dup cartesian-product ; inline
-
-: cartesian-matrix-map ( matrix quot -- matrix' )
-    [ [ first length cartesian-square-indices ] keep ] dip
-    '[ _ @ ] matrix-map ; inline
-
-: cartesian-matrix-column-map ( matrix quot -- matrix' )
-    [ cols first2 ] prepose cartesian-matrix-map ; inline
-
-: cov-matrix-ddof ( matrix ddof -- cov )
-    '[ _ cov-ddof ] cartesian-matrix-column-map ; inline
-
-: population-cov-matrix ( matrix -- cov ) 0 cov-matrix-ddof ; inline
-
-: sample-cov-matrix ( matrix -- cov ) 1 cov-matrix-ddof ; inline
-
-GENERIC: square-rows ( object -- matrix )
-M: integer square-rows <iota> square-rows ;
-M: sequence square-rows
     [ length ] keep >array '[ _ clone ] { } replicate-as ;
 
-GENERIC: square-cols ( object -- matrix )
+M: square-matrix <square-rows> ;
-M: integer square-cols <iota> square-cols ;
+M: matrix <square-rows> >square-matrix ; ! coercing to square is more useful than no-method
-M: sequence square-cols
-    [ length ] keep [ <array> ] with { } map-as ;
 
-: make-matrix-with-indices ( m n quot -- matrix )
+GENERIC: <square-cols> ( desc -- matrix )
-    [ [ <iota> ] bi@ ] dip cartesian-map ; inline
+M: integer <square-cols>
+    <iota> <square-cols> ;
+M: sequence <square-cols>
+    <square-rows> flip ;
 
-: null-matrix? ( matrix -- ? ) empty? ; inline
+M: square-matrix <square-cols> ;
-
+M: matrix <square-cols>
-: well-formed-matrix? ( matrix -- ? )
+    >square-matrix ;
-    [ t ] [
-        [ ] [ first length ] bi
-        '[ length _ = ] all?
-    ] if-empty ;
-
-: dim ( matrix -- pair/f )
-    [ 2 0 <array> ]
-    [ [ length ] [ first length ] bi 2array ] if-empty ;
-
-: square-matrix? ( matrix -- ? )
-    { [ well-formed-matrix? ] [ dim all-eq? ] } 1&& ;
-
-: matrix-coordinates ( dim -- coordinates )
-    first2 [ <iota> ] bi@ cartesian-product ; inline
 
+<PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
 : dimension-range ( matrix -- dim range )
-    dim [ matrix-coordinates ] [ first [1,b] ] bi ;
+    dimension [ <coordinate-matrix> ] [ first [1,b] ] bi ;
 
 : upper-matrix-indices ( matrix -- matrix' )
     dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
 
 : lower-matrix-indices ( matrix -- matrix' )
     dimension-range [ head-slice >array ] 2map concat ;
+PRIVATE>
 
-: make-lower-matrix ( object m n -- matrix )
+! triangulars
-    zero-matrix [ lower-matrix-indices ] [ set-indices ] [ ] tri ;
+DEFER: matrix-set-nths
+: <lower-matrix> ( object m n -- matrix )
+    <zero-matrix> [ lower-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
 
-: make-upper-matrix ( object m n -- matrix )
+: <upper-matrix> ( object m n -- matrix )
-    zero-matrix [ upper-matrix-indices ] [ set-indices ] [ ] tri ;
+    <zero-matrix> [ upper-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
+
+! element- and sequence-wise operations, getters and setters
+: stitch ( m -- m' )
+    [ ] [ [ append ] 2map ] map-reduce ;
+
+: matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
+    '[ _ map ] map ; inline
+
+: column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
+    [ transpose ] dip map transpose ; inline
+
+: matrix-nth ( pair matrix -- elt )
+    [ first2 swap ] dip nth nth ; inline
+
+: matrix-nths ( pairs matrix -- elts )
+    '[ _ matrix-nth ] map ; inline
+
+: matrix-set-nth ( obj pair matrix -- )
+    [ first2 swap ] dip nth set-nth ; inline
+
+: matrix-set-nths ( obj pairs matrix -- )
+    '[ _ matrix-set-nth ] with each ; inline
+
+! -------------------------------------------
+! simple math of matrices follows
+: mneg ( m -- m' ) [ vneg ] map ;
+: mabs ( m -- m' ) [ vabs ] map ;
+
+: n+m ( n m -- m ) [ n+v ] with map ;
+: m+n ( m n -- m ) [ v+n ] curry map ;
+: n-m ( n m -- m ) [ n-v ] with map ;
+: m-n ( m n -- m ) [ v-n ] curry map ;
+: n*m ( n m -- m ) [ n*v ] with map ;
+: m*n ( m n -- m ) [ v*n ] curry map ;
+: n/m ( n m -- m ) [ n/v ] with map ;
+: m/n ( m n -- m ) [ v/n ] curry map ;
+
+: m+  ( m1 m2 -- m ) [ v+ ] 2map ;
+: m-  ( m1 m2 -- m ) [ v- ] 2map ;
+: m*  ( m1 m2 -- m ) [ v* ] 2map ;
+: m/  ( m1 m2 -- m ) [ v/ ] 2map ;
+
+: v.m ( v m -- p ) flip [ v. ] with map ;
+: m.v ( m v -- p ) [ v. ] curry map ;
+: m. ( m m -- m ) flip [ swap m.v ] curry map ;
+
+: m~  ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
+
+: mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
+: mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
+: mnorm ( m -- m' ) dup mmax abs m/n ;
+: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
+: m-1norm ( m -- n ) flip m-infinity-norm ;
+: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
+
+! well-defined for square matrices; but works on nonsquare too
+: main-diagonal ( matrix -- seq )
+    >square-matrix [ swap nth-unsafe ] map-index ; inline
+
+! top right to bottom left; reverse the result if you expected it to start in the lower left
+: anti-diagonal ( matrix -- seq )
+    >square-matrix [ swap nth-end-unsafe ] map-index ; inline
+
+<PRIVATE
+: (rows-iota) ( matrix -- rows-iota )
+    dimension first-unsafe <iota> ;
+: (cols-iota) ( matrix -- cols-iota )
+    dimension second-unsafe <iota> ;
+
+: simple-rows-except ( matrix desc quot -- others )
+    curry [ dup (rows-iota) ] dip
+    pick reject-as swap rows ; inline
+
+: simple-cols-except ( matrix desc quot -- others )
+    curry [ dup (cols-iota) ] dip
+    pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
+
+CONSTANT: scalar-except-quot [ = ]
+CONSTANT: sequence-except-quot [ member? ]
+PRIVATE>
+
+GENERIC: rows-except ( matrix desc -- others )
+M: integer rows-except  scalar-except-quot   simple-rows-except ;
+M: sequence rows-except sequence-except-quot simple-rows-except ;
+
+GENERIC: cols-except ( matrix desc -- others )
+M: integer cols-except  scalar-except-quot   simple-cols-except ;
+M: sequence cols-except sequence-except-quot simple-cols-except ;
+
+! well-defined for any regular matrix
+: matrix-except ( matrix exclude-pair -- submatrix )
+    first2 [ rows-except ] dip cols-except ;
+
+ALIAS: submatrix-excluding matrix-except
+
+:: matrix-except-all ( matrix -- submatrices )
+    matrix dimension [ <iota> ] map first2-unsafe cartesian-product
+    [ [ matrix swap matrix-except ] map ] map ;
+
+ALIAS: all-submatrices matrix-except-all
+
+: dimension ( matrix -- dimension )
+    [ { 0 0 } ]
+    [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;

basis/math/vectors/vectors-docs.factor

@@ -315,11 +315,11 @@ HELP: vclamp
 } ;
 
 HELP: v.
-{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" "a real number" } }
+{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" real } }
 { $description "Computes the dot product of two vectors." } ;
 
 HELP: h.
-{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" "a real number" } }
+{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" real } }
 { $description "Computes the Hermitian inner product of two vectors." } ;
 
 HELP: vs+

basis/math/vectors/vectors-tests.factor

@@ -2,6 +2,7 @@ USING: math.vectors tools.test kernel specialized-arrays compiler
 kernel.private alien.c-types math.functions ;
 SPECIALIZED-ARRAY: int
 
+{ { 10 20 30 } } [ 10 { 1 2 3 } n*v ] unit-test
 { { 1 2 3 } } [ 1/2 { 2 4 6 } n*v ] unit-test
 { { 1 2 3 } } [ { 2 4 6 } 1/2 v*n ] unit-test
 { { 1 2 3 } } [ { 2 4 6 } 2 v/n ] unit-test
@@ -49,3 +50,10 @@ SPECIALIZED-ARRAY: int
 { { 0 30 100 } } [
     { -10 30 120 } { 0 0 0 } { 100 100 100 } vclamp
 ] unit-test
+
+{ { 0 0 1 } } [ { 1 0 0 } { 0 1 0 } cross ] unit-test
+{ { 1 0 0 } } [ { 0 1 0 } { 0 0 1 } cross ] unit-test
+{ { 0 1 0 } } [ { 0 0 1 } { 1 0 0 } cross ] unit-test
+{ { 0.0 -0.707 0.707 } } [ { 1.0 0.0 0.0 } { 0.0 0.707 0.707 } cross ] unit-test
+{ { 0 -2 2 } } [ { -1 -1 -1 } { 1 -1 -1 } cross ] unit-test
+{ { 1 0 0 } } [ { 1 1 0 } { 1 0 0 } proj ] unit-test

basis/math/vectors/vectors.factor

@@ -279,3 +279,19 @@ PRIVATE>
 
 : vclamp ( v min max -- w )
     rot vmin vmax ; inline
+
+: cross ( vec1 vec2 -- vec3 )
+    [ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ]
+    [ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline
+
+:: normal ( vec1 vec2 vec3 -- vec4 )
+    vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
+
+: proj ( v u -- w )
+    [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
+
+: perp ( v u -- w )
+    dupd proj v- ;
+
+: angle-between ( v u -- a )
+    [ normalize ] bi@ h. acos ;

basis/sequences/product/product-docs.factor

@@ -58,7 +58,12 @@ HELP: product-each
 { $description "Calls " { $snippet "quot" } " for every element of the cartesian product of " { $snippet "sequences" } "." }
 { $notes { $snippet "[ ... ] product-each" } " is equivalent to, but more efficient than, " { $snippet "<product-sequence> [ ... ] each" } "." } ;
 
-{ product-map product-each } related-words
+HELP: product-find
+{ $values { "sequences" sequence } { "quot" { $quotation ( ... seq -- ... ? ) } } { "sequence" sequence } }
+{ $description "Calls " { $snippet "quot" } " for every element of the cartesian product of " { $snippet "sequences" } ", returning the first sequence where the quotation returns a true value." }
+{ $notes { $snippet "[ ... ] product-each" } " is equivalent to, but more efficient than, " { $snippet "<product-sequence> [ ... ] find" } "." } ;
+
+{ product-map product-each product-find } related-words
 
 ARTICLE: "sequences.product" "Product sequences"
 "The " { $vocab-link "sequences.product" } " vocabulary provides a virtual sequence and combinators for manipulating the cartesian product of a set of sequences."
@@ -69,6 +74,7 @@ ARTICLE: "sequences.product" "Product sequences"
     product-map-as
     product-map>assoc
     product-each
+    product-find
 } ;
 
 ABOUT: "sequences.product"

basis/sequences/product/product-tests.factor

@@ -1,6 +1,6 @@
 ! Copyright (C) 2009 Joe Groff.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: arrays kernel make sequences sequences.product tools.test ;
+USING: arrays kernel make math sequences sequences.product tools.test ;
 
 { { { 0 "a" } { 1 "a" } { 2 "a" } { 0 "b" } { 1 "b" } { 2 "b" } } }
 [ { { 0 1 2 } { "a" "b" } } <product-sequence> >array ] unit-test
@@ -24,3 +24,13 @@ USING: arrays kernel make sequences sequences.product tools.test ;
 
 { { } } [ { { } { 1 } } [ ] product-map ] unit-test
 { } [ { { } { 1 } } [ drop ] product-each ] unit-test
+
+{ { 2 4 8 } } [
+    { { 1 2 3 } { 4 5 6 } { 7 8 9 } }
+    [ [ even? ] all? ] product-find
+] unit-test
+
+{ f } [
+    { { 1 2 3 } { 4 5 6 } { 7 8 9 } }
+    [ [ 10 > ] all? ] product-find
+] unit-test

basis/sequences/product/product.factor

@@ -76,3 +76,10 @@ M: product-sequence nth
     |[ result |
         sequences [ quot call 2array i result set-nth-unsafe i 1 + i! ] product-each
     ] new-like exemplar assoc-like ; inline
+
+:: product-find ( ... sequences quot: ( ... seq -- ... ? ) -- ... sequence )
+    sequences start-product-iter :> ( ns lengths )
+    lengths [ 0 = ] any? [
+        f [ ns lengths end-product-iter? over or ]
+        [ drop ns sequences nths quot keep and ns lengths product-iter ] until
+    ] unless ; inline

basis/timers/timers-tests.factor

@@ -1,6 +1,6 @@
-USING: timers timers.private calendar concurrency.count-downs
+USING: accessors calendar combinators concurrency.count-downs
 concurrency.promises fry kernel math math.order sequences
-threads tools.test tools.time ;
+threads timers tools.test tools.time ;
 
 { } [
     1 <count-down>
@@ -74,3 +74,20 @@ threads tools.test tools.time ;
     dup restart-timer drop
     700 milliseconds sleep
 ] unit-test
+
+
+{ { 1 } t t t t } [
+    { 0 }
+    dup '[ 0 _ [ 1 + ] change-nth ] 300 milliseconds f <timer>
+    dup start-timer [ thread>> ] keep {
+        [ dup restart-timer thread>> eq? ]
+        [ dup restart-timer thread>> eq? ]
+        [ dup restart-timer thread>> eq? ]
+        [ dup restart-timer thread>> eq? ]
+    } 2cleave
+    700 milliseconds sleep
+] unit-test
+
+[
+    [ ] 1 seconds later start-timer
+] [ timer-already-started? ] must-fail-with

basis/timers/timers.factor

@@ -1,17 +1,17 @@
 ! Copyright (C) 2005, 2008 Slava Pestov, Doug Coleman.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: accessors calendar combinators.short-circuit fry kernel
+
-math math.functions quotations system threads typed ;
+USING: accessors calendar fry kernel math quotations system
+threads ;
+
 IN: timers
 
 TUPLE: timer
     { quot callable initial: [ ] }
-    start-nanos
     delay-nanos
     interval-nanos
-    iteration-start-nanos
+    next-nanos
     quotation-running?
-    restart?
     thread ;
 
 <PRIVATE
@@ -21,60 +21,42 @@ M: f >nanoseconds ;
 M: real >nanoseconds >integer ;
 M: duration >nanoseconds duration>nanoseconds >integer ;
 
-TYPED: set-next-timer-time ( timer: timer -- timer )
+: delay-nanos ( timer -- n )
-    ! start + delay + ceiling((now - (start + delay)) / interval) * interval
+    delay-nanos>> 0 or nano-count + ;
-    nano-count
-    over start-nanos>> -
-    over delay-nanos>> [ - ] when*
-    over interval-nanos>> / ceiling
-    over interval-nanos>> *
-    over start-nanos>> +
-    over delay-nanos>> [ + ] when*
-    >>iteration-start-nanos ;
 
-TYPED: stop-timer? ( timer: timer -- ? )
+: interval-nanos ( timer -- n/f )
-    { [ thread>> self eq? not ] [ restart?>> ] } 1|| ; inline
+    [ next-nanos>> nano-count over - ] [ interval-nanos>> ] bi
+    [ dupd [ mod ] [ swap - ] bi + + ] [ 2drop f ] if* ;
 
-DEFER: call-timer-loop
+: next-nanos ( timer -- timer n/f )
+    dup thread>> self eq? [
+        dup next-nanos>> dup t eq? [
+            drop dup delay-nanos [ >>next-nanos ] keep
+        ] when
+    ] [ f ] if ;
 
-TYPED: loop-timer ( timer: timer -- )
+: run-timer ( timer -- timer )
-    nano-count over
+    dup interval-nanos >>next-nanos
-    [ iteration-start-nanos>> - ] [ interval-nanos>> ] bi <
+    t >>quotation-running?
-    [ set-next-timer-time ] dip
+    dup quot>> call( -- )
-    [ dup iteration-start-nanos>> ] [ 0 ] if
+    f >>quotation-running? ;
-    0 or sleep-until call-timer-loop ;
 
-TYPED: maybe-loop-timer ( timer: timer -- )
+: timer-loop ( timer -- )
-    dup { [ stop-timer? ] [ interval-nanos>> not ] } 1||
+    [ next-nanos ] [
-    [ drop ] [ loop-timer ] if ;
+        dup nano-count <= [
+            drop run-timer yield
+        ] [
+            sleep-until
+        ] if
+    ] while* dup thread>> self eq? [ f >>thread ] when drop ;
 
-TYPED: call-timer-loop ( timer: timer -- )
+: ?interrupt ( thread timer -- )
-    dup stop-timer? [
+    quotation-running?>> [ drop ] [ [ interrupt ] when* ] if ;
-        drop
-    ] [
-        [
-            [ t >>quotation-running? drop ]
-            [ quot>> call( -- ) ]
-            [ f >>quotation-running? drop ] tri
-        ] keep
-        maybe-loop-timer
-    ] if ;
-
-TYPED: sleep-delay ( timer: timer -- )
-    dup stop-timer? [
-        drop
-    ] [
-        nano-count >>start-nanos
-        delay-nanos>> [ sleep ] when*
-    ] if ;
-
-TYPED: timer-loop ( timer: timer -- )
-    [ sleep-delay ]
-    [ nano-count >>iteration-start-nanos call-timer-loop ]
-    [ dup restart?>> [ f >>restart? timer-loop ] [ drop ] if ] tri ;
 
 PRIVATE>
 
+ERROR: timer-already-started timer ;
+
 : <timer> ( quot delay-duration/f interval-duration/f -- timer )
     timer new
         swap >nanoseconds >>interval-nanos
@@ -82,20 +64,19 @@ PRIVATE>
         swap >>quot ; inline
 
 : start-timer ( timer -- )
-    [
+    dup thread>> [ timer-already-started ] when
-        '[ _ timer-loop ] "Timer execution" spawn
+    t >>next-nanos
-    ] keep thread<< ;
+    dup '[ _ timer-loop ] "Timer" <thread>
+    [ >>thread drop ] [ (spawn) ] bi ;
 
 : stop-timer ( timer -- )
-    dup quotation-running?>> [
+    [ f ] change-thread ?interrupt ;
-        dup thread>> [ interrupt ] when*
-    ] unless f >>thread drop ;
 
 : restart-timer ( timer -- )
-    dup quotation-running?>> [
+    dup thread>> [
-        t >>restart? drop
+        t >>next-nanos [ thread>> ] [ ?interrupt ] bi
     ] [
-        dup thread>> [ interrupt ] when* start-timer
+        start-timer
     ] if ;
 
 <PRIVATE
@@ -106,13 +87,10 @@ PRIVATE>
 PRIVATE>
 
 : every ( quot interval-duration -- timer )
-    [ f ] dip (start-timer) ;
+    f swap (start-timer) ;
 
 : later ( quot delay-duration -- timer )
     f (start-timer) ;
 
 : delayed-every ( quot duration -- timer )
     dup (start-timer) ;
-
-: nanos-since ( nano-count -- nanos )
-    [ nano-count ] dip - ;

basis/tools/test/test.factor

@@ -1,14 +1,14 @@
 ! Copyright (C) 2003, 2010 Slava Pestov.
 ! See http://factorcode.org/license.txt for BSD license.
 USING: accessors arrays assocs combinators command-line
-compiler.units constructors continuations debugger effects fry
+compiler.units constructors continuations debugger effects
 generalizations io io.files.temp io.files.unique
-io.streams.string kernel lexer locals macros math.functions
+io.streams.string kernel lexer math math.functions math.vectors
-math.vectors namespaces parser prettyprint quotations sequences
+namespaces parser prettyprint quotations sequences
 sequences.generalizations source-files source-files.errors
 source-files.errors.debugger splitting stack-checker summary
-system tools.errors unicode vocabs vocabs.files vocabs.metadata
+system tools.errors tools.time unicode vocabs vocabs.files
-vocabs.parser words ;
+vocabs.metadata vocabs.parser words ;
 FROM: vocabs.hierarchy => load ;
 IN: tools.test
 
@@ -45,6 +45,7 @@ t restartable-tests? set-global
         swap >>error
         error-continuation get >>continuation ;
 
+INITIALIZED-SYMBOL: long-unit-tests-threshold [ 10,000,000,000 ]
 INITIALIZED-SYMBOL: long-unit-tests-enabled? [ t ]
 
 <PRIVATE
@@ -213,16 +214,26 @@ SYMBOL: forget-tests?
     [ [ [ forget-source ] each ] with-compilation-unit ] [ drop ] if ;
 
 PRIVATE>
+: possible-long-unit-tests ( vocab nanos -- )
+    long-unit-tests-threshold get [
+        dupd > long-unit-tests-enabled? get not and [
+            swap
+            "Warning: possible long unit test for " write
+            vocab-name write " - " write
+            1,000,000,000 /f pprint " seconds" print
+        ] [ 2drop ] if
+    ] [ 2drop ] if* ;
 
 : test-vocab ( vocab -- )
-    lookup-vocab dup [
+    lookup-vocab [
         dup source-loaded?>> [
-            vocab-tests
+            dup vocab-tests [
-            [ [ run-test-file ] each ]
+                [ [ run-test-file ] each ]
-            [ forget-tests ]
+                [ forget-tests ]
-            bi
+                bi
+            ] benchmark possible-long-unit-tests
         ] [ drop ] if
-    ] [ drop ] if ;
+    ] when* ;
 
 : test-vocabs ( vocabs -- ) [ test-vocab ] each ;
 

build.sh

@@ -92,7 +92,7 @@ check_ret() {
 set_downloader() {
     test_program_installed wget
     if [[ $? -ne 0 ]] ; then
-        DOWNLOADER=wget
+        DOWNLOADER="wget -nd"
         DOWNLOADER_NAME=wget
         return
     fi
@@ -154,6 +154,23 @@ clang_version_ok() {
 }
 
 set_cc() {
+
+    # on Cygwin we MUST use the MinGW "cross-compiler", therefore check these first
+    # furthermore, we prefer 64 bit over 32 bit versions if both are available
+    test_programs_installed x86_64-w64-mingw32-gcc x86_64-w64-mingw32-g++
+    if [[ $? -ne 0 ]] ; then
+        [ -z "$CC" ] && CC=x86_64-w64-mingw32-gcc
+        [ -z "$CXX" ] && CXX=x86_64-w64-mingw32-g++
+        return
+    fi
+
+    test_programs_installed i686-w64-mingw32-gcc i686-w64-mingw32-g++
+    if [[ $? -ne 0 ]] ; then
+        [ -z "$CC" ] && CC=i686-w64-mingw32-gcc
+        [ -z "$CXX" ] && CXX=i686-w64-mingw32-g++
+        return
+    fi
+
     test_programs_installed clang clang++
     if [[ $? -ne 0 ]] && clang_version_ok ; then
         [ -z "$CC" ] && CC=clang
@@ -161,6 +178,7 @@ set_cc() {
         return
     fi
 
+    # gcc and g++ will fail to correctly build Factor on Cygwin
     test_programs_installed gcc g++
     if [[ $? -ne 0 ]] ; then
         [ -z "$CC" ] && CC=gcc
@@ -187,8 +205,8 @@ check_installed_programs() {
     ensure_program_installed uname
     ensure_program_installed git
     ensure_program_installed wget curl
-    ensure_program_installed clang gcc
+    ensure_program_installed clang x86_64-w64-mingw32-gcc i686-w64-mingw32-gcc gcc
-    ensure_program_installed clang++ g++ cl
+    ensure_program_installed clang++ x86_64-w64-mingw32-g++ i686-w64-mingw32-g++ g++ cl
     ensure_program_installed make gmake
     ensure_program_installed md5sum md5
     ensure_program_installed cut

core/io/io.factor

@@ -105,9 +105,6 @@ SYMBOL: error-stream
 
 : bl ( -- ) output-stream get stream-bl ;
 
-: each-morsel ( ..a handler: ( ..a data -- ..b ) reader: ( ..b -- ..a data ) -- ..a )
-    [ dup ] compose swap while drop ; inline
-
 <PRIVATE
 
 : stream-exemplar ( stream -- exemplar )
@@ -158,7 +155,7 @@ ERROR: invalid-read-buffer buf stream ;
     input-stream get stream-read-partial-into ; inline
 
 : each-stream-line ( ... stream quot: ( ... line -- ... ) -- ... )
-    swap [ stream-readln ] curry each-morsel ; inline
+    [ [ stream-readln ] curry ] dip while* ; inline
 
 : each-line ( ... quot: ( ... line -- ... ) -- ... )
     input-stream get swap each-stream-line ; inline
@@ -174,15 +171,16 @@ ERROR: invalid-read-buffer buf stream ;
 CONSTANT: each-block-size 65536
 
 : (each-stream-block-slice) ( ... stream quot: ( ... block-slice -- ... ) block-size -- ... )
-    [ [ drop ] prepose swap ] dip
+    -rot [
-    [ swap (new-sequence-for-stream) ] keepd
+        [ (new-sequence-for-stream) ] keep
-    [ stream-read-partial-into ] 2curry each-morsel drop ; inline
+        [ stream-read-partial-into ] 2curry
+    ] dip while drop ; inline
 
 : each-stream-block-slice ( ... stream quot: ( ... block-slice -- ... ) -- ... )
     each-block-size (each-stream-block-slice) ; inline
 
 : (each-stream-block) ( ... stream quot: ( ... block -- ... ) block-size -- ... )
-    rot [ stream-read-partial ] 2curry each-morsel ; inline
+    -rot [ [ stream-read-partial ] 2curry ] dip while* ; inline
 
 : each-stream-block ( ... stream quot: ( ... block -- ... ) -- ... )
     each-block-size (each-stream-block) ; inline

core/kernel/kernel-docs.factor

@@ -885,6 +885,10 @@ HELP: while
 { $values { "pred" { $quotation ( ..a -- ..b ? ) } } { "body" { $quotation ( ..b -- ..a ) } } }
 { $description "Calls " { $snippet "body" } " until " { $snippet "pred" } " returns " { $link f } "." } ;
 
+HELP: while*
+{ $values { "pred" { $quotation ( ..a -- ..b ? ) } } { "body" { $quotation ( ..b ? -- ..a ) } } }
+{ $description "Calls " { $snippet "body" } " until " { $snippet "pred" } " returns " { $link f } "." } ;
+
 HELP: until
 { $values { "pred" { $quotation ( ..a -- ..b ? ) } } { "body" { $quotation ( ..b -- ..a ) } } }
 { $description "Calls " { $snippet "body" } " until " { $snippet "pred" } " returns " { $link t } "." } ;

core/kernel/kernel.factor

@@ -310,6 +310,9 @@ UNION: boolean t postpone: f ;
 : while ( ..a pred: ( ..a -- ..b ? ) body: ( ..b -- ..a ) -- ..b )
     swap do compose [ loop ] curry when ; inline
 
+: while* ( ..a pred: ( ..a -- ..b ? ) body: ( ..b ? -- ..a ) -- ..b )
+    [ [ dup ] compose ] dip while drop ; inline
+
 : until ( ..a pred: ( ..a -- ..b ? ) body: ( ..b -- ..a ) -- ..b )
     [ [ not ] compose ] dip while ; inline
 

core/sequences/sequences-docs.factor

@@ -1620,6 +1620,10 @@ HELP: assert-sequence=
   }
 } ;
 
+HELP: cartesian-find
+{ $values { "seq1" sequence } { "seq2" sequence } { "quot" { $quotation ( ... elt1 elt2 -- ... ? ) } } { "elt1" object } { "elt2" object } }
+{ $description "Applies the quotation to every possible pairing of elements from the two sequences, returning the first two elements where the quotation returns a true value." } ;
+
 HELP: cartesian-each
 { $values { "seq1" sequence } { "seq2" sequence } { "quot" { $quotation ( ... elt1 elt2 -- ... ) } } }
 { $description "Applies the quotation to every possible pairing of elements from the two sequences." } ;
@@ -1980,6 +1984,7 @@ $nl
 { $subsections
     cartesian-each
     cartesian-map
+    cartesian-find
 }
 "Computing the cartesian product of two sequences:"
 { $subsections

core/sequences/sequences-tests.factor

@@ -363,6 +363,9 @@ M: bogus-hashcode hashcode* 2drop 0 >bignum ;
 { { { { 1 "a" } { 1 "b" } } { { 2 "a" } { 2 "b" } } } }
 [ { 1 2 } { "a" "b" } cartesian-product ] unit-test
 
+{ 2 4 } [ { 1 2 3 } { 4 5 6 } [ [ even? ] both? ] cartesian-find ] unit-test
+{ f f } [ { 1 2 3 } { 4 5 6 } [ [ 10 > ] both? ] cartesian-find ] unit-test
+
 [ { } [ string>digits sum ] [ + ] map-reduce ] must-infer
 [ { } [ ] [ + ] map-reduce ] must-fail
 { 4 } [ { 1 1 } [ 1 + ] [ + ] map-reduce ] unit-test

core/sequences/sequences.factor

@@ -1122,6 +1122,9 @@ M: repetition sum [ elt>> ] [ length>> ] bi * ; inline
 : cartesian-product ( seq1 seq2 -- newseq )
     [ { } 2sequence ] cartesian-map ;
 
+: cartesian-find ( ... seq1 seq2 quot: ( ... elt1 elt2 -- ... ? ) -- ... elt1 elt2 )
+    [ f ] 3dip [ with find swap ] 2curry [ nip ] prepose find nip swap ; inline
+
 <PRIVATE
 
 : select-by ( ... seq quot: ( ... elt -- ... x ) compare: ( obj1 obj2 -- ? ) -- ... elt )

extra/audio/wav/wav.factor

@@ -38,12 +38,12 @@ STRUCT: wav-data-chunk
 
 :: read-wav-chunks ( -- fmt data )
     f :> fmt! f :> data!
-    [ { [ fmt data and not ] [ read-chunk ] } 0&& dup ]
+    [ { [ fmt data and not ] [ read-chunk ] } 0&& ]
     [ {
         { [ dup FMT-MAGIC  wav-fmt-chunk  check-chunk ] [ wav-fmt-chunk  memory>struct fmt!  ] }
         { [ dup DATA-MAGIC wav-data-chunk check-chunk ] [ wav-data-chunk memory>struct data! ] }
         [ drop ]
-    } cond ] while drop
+    } cond ] while*
     fmt data 2dup and [ invalid-audio-file ] unless ;
 
 : verify-wav ( chunk -- )

extra/benchmark/3d-matrix-scalar/3d-matrix-scalar.factor

@@ -1,15 +1,15 @@
-USING: kernel locals math math.matrices math.order math.vectors
+USING: kernel locals math math.matrices math.matrices.extras
-prettyprint sequences ;
+math.order math.vectors prettyprint sequences ;
 IN: benchmark.3d-matrix-scalar
 
 :: p-matrix ( dim fov near far -- matrix )
     dim dup first2 min v/n fov v*n near v*n
-    near far frustum-matrix4 ;
+    near far <frustum-matrix4> ;
 
 :: mv-matrix ( pitch yaw location -- matrix )
-    { 1.0 0.0 0.0 } pitch rotation-matrix4
+    { 1.0 0.0 0.0 } pitch <rotation-matrix4>
-    { 0.0 1.0 0.0 } yaw   rotation-matrix4
+    { 0.0 1.0 0.0 } yaw   <rotation-matrix4>
-    location vneg translation-matrix4 m. m. ;
+    location vneg <translation-matrix4> m. m. ;
 
 :: 3d-matrix-scalar-benchmark ( -- )
     f :> result!

extra/benchmark/3d-matrix-vector/3d-matrix-vector.factor

@@ -1,5 +1,6 @@
-USING: kernel locals math math.matrices.simd math.order math.vectors
+USING: kernel locals math math.matrices math.matrices.simd
-math.vectors.simd prettyprint sequences typed ;
+math.order math.vectors math.vectors.simd prettyprint sequences
+typed ;
 QUALIFIED-WITH: alien.c-types c
 IN: benchmark.3d-matrix-vector
 

extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor

@@ -1,4 +1,4 @@
-USING: locals math math.combinatorics math.matrices
+USING: locals math math.combinatorics math.matrices math.matrices.extras
 prettyprint sequences typed ;
 IN: benchmark.matrix-exponential-scalar
 
@@ -15,7 +15,7 @@ IN: benchmark.matrix-exponential-scalar
 
 :: matrix-exponential-scalar-benchmark ( -- )
     f :> result!
-    4 identity-matrix :> i4
+    4 <identity-matrix> :> i4
     10000 [
         i4 20 e^m result!
     ] times

extra/cuesheet/cuesheet.factor

@@ -127,7 +127,7 @@ ERROR: unknown-syntax syntax ;
 PRIVATE>
 
 : read-cuesheet ( -- cuesheet )
-    <cuesheet> [ readln dup ] [ parse-line ] while drop ;
+    <cuesheet> [ readln ] [ parse-line ] while* ;
 
 : file>cuesheet ( path -- cuesheet )
     utf8 [ read-cuesheet ] with-file-reader ;

extra/game/debug/tests/tests.factor

@@ -2,8 +2,9 @@
 ! See http://factorcode.org/license.txt for BSD license.
 USING: accessors colors.constants game.debug game.loop
 game.worlds gpu gpu.framebuffers gpu.util.wasd kernel literals
-locals make math math.matrices math.parser math.trig sequences
+locals make math math.matrices math.matrices.extras math.parser
-specialized-arrays ui.gadgets.worlds ui.pixel-formats ;
+math.trig sequences specialized-arrays ui.gadgets.worlds
+ui.pixel-formats ;
 FROM: alien.c-types => float ;
 SPECIALIZED-ARRAY: float
 IN: game.debug.tests
@@ -21,7 +22,7 @@ IN: game.debug.tests
         { 0 0 0 } { 1 0 0 } color: red   debug-line
         { 0 0 0 } { 0 1 0 } color: green debug-line
         { 0 0 0 } { 0 0 1 } color: blue  debug-line
-        { -1.2 0 0 } { 0 1 0 } 0 deg>rad rotation-matrix3 debug-axes
+        { -1.2 0 0 } { 0 1 0 } 0 deg>rad <rotation-matrix3> debug-axes
         { 3 5 -2 } { 3 2 1 } color: white debug-box
         { 0 9 0 } 8 2 color: blue debug-cylinder
     ] float-array{ } make

extra/gml/geometry/geometry.factor

@@ -1,8 +1,8 @@
 ! Copyright (C) 2010 Slava Pestov.
-USING: arrays kernel math.matrices math.vectors.simd.cords
+USING: arrays gml.runtime kernel math.matrices
-math.trig gml.runtime ;
+math.matrices.extras math.trig math.vectors.simd.cords   ;
 IN: gml.geometry
 
 GML: rot_vec ( v n alpha -- v )
     ! Inefficient!
-    deg>rad rotation-matrix4 swap >array m.v >double-4 ;
+    deg>rad <rotation-matrix4> swap >array m.v >double-4 ;

extra/gpu/demos/raytrace/raytrace.factor

@@ -2,9 +2,9 @@
 ! See http://factorcode.org/license.txt for BSD license.
 USING: accessors arrays combinators.tuple game.loop game.worlds
 generalizations gpu gpu.render gpu.shaders gpu.util gpu.util.wasd
-kernel literals math math.libm math.matrices math.order math.vectors
+kernel literals math math.libm math.matrices math.matrices.extras
-method-chains sequences ui ui.gadgets ui.gadgets.worlds
+math.order math.vectors method-chains sequences ui ui.gadgets
-ui.pixel-formats audio.engine audio.loader locals ;
+ui.gadgets.worlds ui.pixel-formats audio.engine audio.loader locals ;
 IN: gpu.demos.raytrace
 
 GLSL-SHADER-FILE: raytrace-vertex-shader vertex-shader "raytrace.v.glsl"
@@ -48,7 +48,7 @@ TUPLE: raytrace-world < wasd-world
     dup dtheta>> [ + ] curry change-theta drop ;
 
 : sphere-center ( sphere -- center )
-    [ [ axis>> ] [ theta>> ] bi rotation-matrix4 ]
+    [ [ axis>> ] [ theta>> ] bi <rotation-matrix4> ]
     [ home>> ] bi m.v ;
 
 M: sphere audio-position sphere-center ; inline

extra/gpu/util/wasd/wasd.factor

@@ -4,8 +4,8 @@ USING: accessors arrays combinators.smart game.input
 game.input.scancodes game.loop game.worlds
 gpu.render gpu.state kernel literals
 locals math math.constants math.functions math.matrices
-math.order math.vectors opengl.gl sequences
+math.matrices.extras math.order math.vectors opengl.gl
-ui ui.gadgets.worlds specialized-arrays audio.engine ;
+sequences ui ui.gadgets.worlds specialized-arrays audio.engine ;
 FROM: alien.c-types => float ;
 SPECIALIZED-ARRAY: float
 IN: gpu.util.wasd
@@ -38,14 +38,14 @@ GENERIC: wasd-fly-vertically? ( world -- ? )
 M: wasd-world wasd-fly-vertically? drop t ;
 
 : wasd-mv-matrix ( world -- matrix )
-    [ { 1.0 0.0 0.0 } swap pitch>> rotation-matrix4 ]
+    [ { 1.0 0.0 0.0 } swap pitch>> <rotation-matrix4> ]
-    [ { 0.0 1.0 0.0 } swap yaw>>   rotation-matrix4 ]
+    [ { 0.0 1.0 0.0 } swap yaw>>   <rotation-matrix4> ]
-    [ location>> vneg translation-matrix4 ] tri m. m. ;
+    [ location>> vneg <translation-matrix4> ] tri m. m. ;
 
 : wasd-mv-inv-matrix ( world -- matrix )
-    [ location>> translation-matrix4 ]
+    [ location>> <translation-matrix4> ]
-    [ {  0.0 -1.0 0.0 } swap yaw>>   rotation-matrix4 ]
+    [ {  0.0 -1.0 0.0 } swap yaw>>   <rotation-matrix4> ]
-    [ { -1.0  0.0 0.0 } swap pitch>> rotation-matrix4 ] tri m. m. ;
+    [ { -1.0  0.0 0.0 } swap pitch>> <rotation-matrix4> ] tri m. m. ;
 
 : wasd-p-matrix ( world -- matrix )
     p-matrix>> ;
@@ -63,7 +63,7 @@ CONSTANT: fov 0.7
     world wasd-far-plane :> far-plane
 
     world wasd-fov-vector near-plane v*n
-    near-plane far-plane frustum-matrix4 ;
+    near-plane far-plane <frustum-matrix4> ;
 
 :: wasd-pixel-ray ( world loc -- direction )
     loc world dim>> [ /f 0.5 - 2.0 * ] 2map

extra/html/entities/entities.factor

@@ -28,7 +28,7 @@ PRIVATE>
 
 : html-escape ( str -- newstr )
     [
-        [ dup next-escape dup ] [ escape, ] while 2drop ,
+        [ dup next-escape ] [ escape, ] while* drop ,
     ] { } make dup length 1 > [ concat ] [ first ] if ;
 
 <PRIVATE

extra/images/atlas/atlas.factor

@@ -47,7 +47,7 @@ ERROR: atlas-image-formats-dont-match images ;
 :: (pack-images) ( images atlas-width sort-quot -- placements )
     images sort-quot inv-sort-with [ f image-placement boa ] map :> image-placements
     0 :> @y!
-    [ image-placements atlas-width @y (pack-stripe) dup ] [ @y + @y! ] while drop
+    [ image-placements atlas-width @y (pack-stripe) ] [ @y + @y! ] while*
     image-placements ; inline
 
 : atlas-image-format ( image-placements -- component-order component-type upside-down? )

extra/math/matrices/elimination/elimination-docs.factor

@@ -8,8 +8,9 @@ HELP: inverse
 { $examples
   "A matrix multiplied by its inverse is the identity matrix."
   { $example
-    "USING: kernel math.matrices math.matrices.elimination prettyprint ;"
+    "USING: kernel math.matrices prettyprint ;"
-    "{ { 3 4 } { 7 9 } } dup inverse m. 2 identity-matrix = ."
+    "FROM: math.matrices.elimination => inverse ;"
+    "{ { 3 4 } { 7 9 } } dup inverse m. 2 <identity-matrix> = ."
     "t"
   }
 } ;

extra/math/matrices/elimination/elimination.factor

@@ -93,7 +93,7 @@ SYMBOL: matrix
 
 : nullspace ( matrix -- seq )
     echelon reduced dup empty? [
-        dup first length identity-matrix [
+        dup first length <identity-matrix> [
             [
                 dup leading drop
                 [ basis-vector ] [ drop ] if*
@@ -109,5 +109,5 @@ SYMBOL: matrix
 
 : inverse ( matrix -- matrix ) ! Assumes an invertible matrix
     dup length
-    [ identity-matrix [ append ] 2map solution ] keep
+    [ <identity-matrix> [ append ] 2map solution ] keep
     [ tail ] curry map ;

extra/math/matrices/extras/authors.txt

@@ -0,0 +1,4 @@
+Slava Pestov
+Joe Groff
+Doug Coleman
+Cat Stevens

extra/math/matrices/extras/extras-docs.factor

@@ -0,0 +1,641 @@
+USING: arrays generic.single help.markup help.syntax kernel math
+math.matrices math.matrices.private math.matrices.extras
+math.order math.ratios math.vectors opengl.gl random sequences
+urls ;
+IN: math.matrices.extras
+
+ABOUT: "math.matrices.extras"
+
+ARTICLE: "math.matrices.extras" "Extra matrix operations"
+
+"These constructions have special mathematical properties:"
+{ $subsections
+    <box-matrix>
+    <hankel-matrix>
+    <hilbert-matrix>
+    <toeplitz-matrix>
+    <vandermonde-matrix>
+}
+
+"Common transformation matrices:"
+{ $subsections
+    <frustum-matrix4>
+    <ortho-matrix4>
+    <rotation-matrix3>
+    <rotation-matrix4>
+    <scale-matrix3>
+    <scale-matrix4>
+    <skew-matrix4>
+    <translation-matrix4>
+
+    <random-integer-matrix>
+    <random-unit-matrix>
+
+}
+
+
+{ $subsections
+    invertible-matrix?
+    linearly-independent-matrix?
+}
+
+"Common algorithms on matrices:"
+{ $subsections
+    gram-schmidt
+    gram-schmidt-normalize
+    kronecker-product
+    outer-product
+}
+
+"Matrix algebra:"
+{ $subsections
+    mmin
+    mmax
+    mnorm
+    rank
+    nullity
+
+} { $subsections
+    determinant 1/det m*1/det
+    >minors >cofactors
+    multiplicative-inverse
+}
+
+"Covariance in matrices:"
+{ $subsections
+    covariance-matrix
+    covariance-matrix-ddof
+    sample-covariance-matrix
+}
+
+"Errors thrown by this vocabulary:"
+{ $subsections negative-power-matrix non-square-determinant undefined-inverse } ;
+
+HELP: invertible-matrix?
+{ $values { "matrix" matrix } { "?" boolean } }
+{ $description "Tests whether the input matrix has a " { $link multiplicative-inverse } ". In order for a matrix to be invertible, it must be a " { $link square-matrix } ", " { $emphasis "or" } ", if it is non-square, it must not be of " { $link +deficient-rank+ } "." }
+{ $examples { $example "USING: math.matrices.extras prettyprint ;" "" } } ;
+
+HELP: linearly-independent-matrix?
+{ $values { "matrix" matrix } { "?" boolean } }
+{ $description "Tests whether the input matrix is linearly independent." }
+{ $examples { $example "USING: math.matrices.extras prettyprint ;" "" } } ;
+
+! SINGLETON RANK TYPES
+HELP: rank-kind
+{ $class-description "The class of matrix rank quantifiers." } ;
+
+HELP: +full-rank+
+{ $class-description "A " { $link rank-kind } " describing a matrix of full rank." } ;
+HELP: +half-rank+
+{ $class-description "A " { $link rank-kind } " describing a matrix of half rank." } ;
+HELP: +zero-rank+
+{ $class-description "A " { $link rank-kind } " describing a matrix of zero rank." } ;
+HELP: +deficient-rank+
+{ $class-description "A " { $link rank-kind } " describing a matrix of deficient rank." } ;
+HELP: +uncalculated-rank+
+{ $class-description "A " { $link rank-kind } " describing a matrix whose rank is not (yet) known." } ;
+
+! ERRORS
+
+HELP: negative-power-matrix
+{ $values { "m" matrix } { "n" integer } }
+{ $description "Throws a " { $link negative-power-matrix } " error." }
+{ $error-description "Given the semantics of " { $link m^n } ", negative exponents are not within the domain of the power matrix function." } ;
+
+HELP: non-square-determinant
+{ $values { "m" integer } { "n" integer } }
+{ $description "Throws a " { $link non-square-determinant } " error." }
+{ $error-description { $link determinant } " was used with a non-square matrix whose dimensions are " { $snippet "m x n" } ". It is not generally possible to find the determinant of a non-square matrix." } ;
+
+HELP: undefined-inverse
+{ $values { "m" integer } { "n" integer } { "r" rank-kind } }
+{ $description "Throws an " { $link undefined-inverse } " error." }
+{ $error-description { $link multiplicative-inverse } " was used with a non-square matrix of rank " { $snippet "rank" } " whose dimensions are " { $snippet "m x n" } ". It is not generally possible to find the inverse of a " { $link +deficient-rank+ } " non-square " { $link matrix } "." } ;
+
+HELP: <random-integer-matrix>
+{ $values { "m" integer } { "n" integer } { "max" integer } { "matrix" matrix } }
+{ $description "Creates a " { $snippet "m x n" } " " { $link matrix } " full of random, possibly signed " { $link integer } "s whose absolute values are less than or equal to " { $snippet "max" } ", as given by " { $link random-integers } "." }
+{ $notelist
+    { "The signedness of the numbers in the resulting matrix will be randomized. Use " { $link mabs } " with this word to generate a matrix of random positive integers." }
+    { $equiv-word-note "integral" <random-unit-matrix> }
+}
+{ $errors { $link no-method } " if " { $snippet "max"} " is not an " { $link integer } "." }
+{ $examples
+    { $unchecked-example
+        "USING: math.matrices.extras prettyprint ;"
+        "2 4 15 <random-integer-matrix> ."
+        "{ { -9 -9 1 3 } { -14 -8 14 10 } }"
+    }
+} ;
+
+HELP: <random-unit-matrix>
+{ $values { "m" integer } { "n" integer } { "max" number } { "matrix" matrix } }
+{ $description "Creates a " { $snippet "m x n" } " " { $link matrix } " full of random, possibly signed " { $link float } "s  as a fraction of " { $snippet "max" } "." }
+{ $notelist
+    { "The signedness of the numbers in the resulting matrix will be randomized. Use " { $link mabs } " with this word to generate a matrix of random positive numbers." }
+    { $equiv-word-note "real" <random-integer-matrix> }
+    { "This word is implemented by generating sub-integral floats through " { $link random-units } " and multiplying by random integers less than or equal to " { $snippet "max" } "." }
+}
+{ $examples
+    { $unchecked-example
+        "USING: math.matrices.extras prettyprint ;"
+        "4 2 15 <random-unit-matrix> ."
+"{
+    { -3.713295909201797 3.815787135075961 }
+    { -2.460506890603817 1.535222788710546 }
+    { 3.692213981267878 -1.462963244399762 }
+    { 13.8967592095433 -6.688509969360172 }
+}"
+    }
+} ;
+
+
+
+HELP: <hankel-matrix>
+{ $values { "n" integer } { "matrix" matrix } }
+{ $description
+    "A Hankel matrix is a symmetric, " { $link square-matrix } " in which each ascending skew-diagonal from left to right is constant. See " { $url URL" https://en.wikipedia.org/wiki/Hankel_matrix" "hankel matrix" } "."
+    $nl
+    "The following is true of any Hankel matrix" { $snippet "A" } ": " { $snippet "A[i][j] = A[j][i] = a[i+j-2]" } "."
+    $nl
+    "The " { $link <toeplitz-matrix> } " is an upside-down Hankel matrix."
+    $nl
+    "The " { $link <hilbert-matrix> } " is a special case of the Hankel matrix."
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "4 <hankel-matrix> ."
+        "{ { 1 2 3 4 } { 2 3 4 0 } { 3 4 0 0 } { 4 0 0 0 } }"
+    }
+} ;
+
+HELP: <hilbert-matrix>
+{ $values { "m" integer } { "n" integer } { "matrix" matrix } }
+{ $description
+    "A Hilbert matrix is a " { $link square-matrix } " " { $snippet "A" } " in which entries are the unit fractions "
+    { $snippet "A[i][j] = 1/(i+j-1)" }
+    ". See " { $url URL" https://en.wikipedia.org/wiki/Hilbert_matrix" "hilbert matrix" } "."
+    $nl
+    "A Hilbert matrix is a special case of the " { $link <hankel-matrix> } "."
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "1 2 <hilbert-matrix> ."
+        "{ { 1 1/2 } }"
+    }
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "3 6 <hilbert-matrix> ."
+"{
+    { 1 1/2 1/3 1/4 1/5 1/6 }
+    { 1/2 1/3 1/4 1/5 1/6 1/7 }
+    { 1/3 1/4 1/5 1/6 1/7 1/8 }
+}"
+    }
+} ;
+
+HELP: <toeplitz-matrix>
+{ $values { "n" integer } { "matrix" matrix } }
+{ $description "A Toeplitz matrix is an upside-down " { $link <hankel-matrix> } ". Unlike the Hankel matrix, a Toeplitz matrix can be non-square. See " { $url URL" https://en.wikipedia.org/wiki/Hankel_matrix" "hankel matrix" } "."
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "4 <toeplitz-matrix> ."
+        "{ { 1 2 3 4 } { 2 1 2 3 } { 3 2 1 2 } { 4 3 2 1 } }"
+    }
+} ;
+
+HELP: <box-matrix>
+{ $values { "r" integer } { "matrix" matrix } }
+{ $description "Create a box matrix (a " { $link square-matrix } ") with the dimensions of " { $snippet "r x r" } ", filled with ones. The number of elements in the output scales linearly (" { $snippet "(r*2)+1" } ") with " { $snippet "r" } "." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "2 <box-matrix> ."
+"{
+    { 1 1 1 1 1 }
+    { 1 1 1 1 1 }
+    { 1 1 1 1 1 }
+    { 1 1 1 1 1 }
+    { 1 1 1 1 1 }
+}"
+    }
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "3 <box-matrix> ."
+"{
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+    { 1 1 1 1 1 1 1 }
+}"
+    }
+
+} ;
+
+HELP: <scale-matrix3>
+{ $values { "factors" sequence } { "matrix" matrix } }
+{ $description "Make a " { $snippet "3 x 3" } " scaling matrix, used to scale an object in 3 dimensions. See " { $url URL" https://en.wikipedia.org/wiki/Scaling_(geometry)#Matrix_representation" "scaling matrix on Wikipedia" } "." }
+{ $notelist
+    { $finite-input-note "three" "factors" }
+    { $equiv-word-note "3-matrix" <scale-matrix4> }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ 22 33 -44 } <scale-matrix4> ."
+"{
+    { 22 0.0 0.0 0.0 }
+    { 0.0 33 0.0 0.0 }
+    { 0.0 0.0 -44 0.0 }
+    { 0.0 0.0 0.0 1.0 }
+}"
+    }
+} ;
+
+HELP: <scale-matrix4>
+{ $values { "factors" sequence } { "matrix" matrix } }
+{ $description "Make a " { $snippet "4 x 4" } " scaling matrix, used to scale an object in 3 or more dimensions. See " { $url URL" https://en.wikipedia.org/wiki/Scaling_(geometry)#Matrix_representation" "scaling matrix on Wikipedia" } "." }
+{ $notelist
+    { $finite-input-note "three" "factors" }
+    { $equiv-word-note "4-matrix" <scale-matrix3> }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ 22 33 -44 } <scale-matrix4> ."
+"{
+    { 22 0.0 0.0 0.0 }
+    { 0.0 33 0.0 0.0 }
+    { 0.0 0.0 -44 0.0 }
+    { 0.0 0.0 0.0 1.0 }
+}"
+    }
+} ;
+
+HELP: <ortho-matrix4>
+{ $values { "factors" sequence } { "matrix" matrix } }
+{ $description "Create a " { $link <scale-matrix4> } ", with the scale factors inverted." }
+{ $notelist
+    { $finite-input-note "three" "factors" }
+    { $equiv-word-note "inverse" <scale-matrix4> }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ -9.3 100 1/2 } <ortho-matrix4> ."
+"{
+    { -0.1075268817204301 0.0 0.0 0.0 }
+    { 0.0 1/100 0.0 0.0 }
+    { 0.0 0.0 2 0.0 }
+    { 0.0 0.0 0.0 1.0 }
+}"
+    }
+} ;
+
+HELP: <frustum-matrix4>
+{ $values { "xy-dim" pair } { "near" number } { "far" number } { "matrix" matrix } }
+{ $description "Make a " { $snippet "4 x 4" } " matrix suitable for representing an occlusion frustum. A viewing or occlusion frustum is the three-dimensional region of a three-dimensional object which is visible on the screen. See " { $url URL" https://en.wikipedia.org/wiki/Frustum" "frustum on Wikipedia" } "." }
+{ $notes { $finite-input-note "two" "xy-dim" } }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ 5 4 } 5 6 <frustum-matrix4> ."
+"{
+    { 1.0 0.0 0.0 0.0 }
+    { 0.0 1.25 0.0 0.0 }
+    { 0.0 0.0 -11.0 -60.0 }
+    { 0.0 0.0 -1.0 0.0 }
+}"
+    }
+} ;
+{ <frustum-matrix4> glFrustum } related-words
+
+HELP: cartesian-matrix-map
+{ $values { "matrix" matrix } { "quot" { $quotation ( ... pair matrix -- ... matrix' ) } } { "matrix-seq" { $sequence matrix } } }
+{ $description "Calls the quotation with the matrix and the coordinate pair of the current element on the stack, with the matrix on the top of the stack." }
+{ $examples
+  { $example
+    "USING: arrays math.matrices.extras prettyprint ;"
+    "{ { 21 22 } { 23 24 } } [ 2array ] cartesian-matrix-map ."
+"{
+    {
+        { { 0 0 } { { 21 22 } { 23 24 } } }
+        { { 0 1 } { { 21 22 } { 23 24 } } }
+    }
+    {
+        { { 1 0 } { { 21 22 } { 23 24 } } }
+        { { 1 1 } { { 21 22 } { 23 24 } } }
+    }
+}"
+  }
+}
+{ $notelist
+  { $equiv-word-note "orthogonal" cartesian-column-map }
+  { $equiv-word-note "two-dimensional" map-index }
+  $2d-only-note
+} ;
+
+HELP: cartesian-column-map
+{ $values { "matrix" matrix } { "quot" { $quotation ( ... pair matrix -- ... matrix' ) } } { "matrix-seq" { $sequence matrix } } }
+{ $notelist
+  { $equiv-word-note "orthogonal" cartesian-matrix-map }
+  $2d-only-note
+} ;
+
+HELP: gram-schmidt
+{ $values { "matrix" matrix } { "orthogonal" matrix } }
+{ $description "Apply a Gram-Schmidt transform on the matrix." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ { 1 2 } { 3 4 } { 5 6 } } gram-schmidt ."
+        "{ { 1 2 } { 4/5 -2/5 } { 0 0 } }"
+    }
+} ;
+
+HELP: gram-schmidt-normalize
+{ $values { "matrix" matrix } { "orthonormal" matrix } }
+{ $description "Apply a Gram-Schmidt transform on the matrix, and " { $link normalize } " each row of the result, resulting in an orthogonal and normalized matrix (orthonormal)." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ { 1 2 } { 3 4 } { 5 6 } } gram-schmidt-normalize ."
+"{
+    { 0.4472135954999579 0.8944271909999159 }
+    { 0.894427190999916 -0.447213595499958 }
+    { NAN: 8000000000000 NAN: 8000000000000 }
+}"
+    }
+} ;
+
+HELP: m^n
+{ $values { "m" matrix } { "n" object } }
+{ $description "Compute the " { $snippet "nth" } " power of the input matrix. If " { $snippet "n" } " is " { $snippet "-1" } ", the inverse of the matrix is calculated (but see " { $link multiplicative-inverse } " for pitfalls)." }
+{ $errors
+    { $link negative-power-matrix } " if " { $snippet "n" } " is a negative number other than " { $snippet "-1" } "."
+    $nl
+    { $link undefined-inverse } " if " { $snippet "n" } " is " { $snippet "-1" } " and the " { $link multiplicative-inverse } " of " { $snippet "m" } " is undefined."
+}
+{ $notelist
+    { $equiv-word-note "swapped" n^m }
+    $2d-only-note
+    { $matrix-scalar-note max abs / }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ { 1 2 } { 3 4 } } 2 m^n ."
+        "{ { 7 10 } { 15 22 } }"
+    }
+} ;
+
+HELP: n^m
+{ $values { "n" object } { "m" matrix } }
+{ $description "Because it is nonsensical to raise a number to the power of a matrix, this word exists to save typing " { $snippet "swap m^n" } ". See " { $link m^n } " for more information." }
+{ $errors
+    { $link negative-power-matrix } " if " { $snippet "n" } " is a negative number other than " { $snippet "-1" } "."
+    $nl
+    { $link undefined-inverse } " if " { $snippet "n" } " is " { $snippet "-1" } " and the " { $link multiplicative-inverse } " of " { $snippet "m" } " is undefined."
+}
+{ $notelist
+    { $equiv-word-note "swapped" m^n }
+    $2d-only-note
+    { $matrix-scalar-note max abs / }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "2 { { 1 2 } { 3 4 } } n^m ."
+        "{ { 7 10 } { 15 22 } }"
+    }
+} ;
+
+HELP: kronecker-product
+{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
+{ $description "Calculates the " { $url URL" http://enwp.org/Kronecker_product" "Kronecker product" } " of two matrices. This product can be described as a generalization of the vector-based " { $link outer-product } " to matrices. The Kronecker product gives the matrix of the tensor product with respect to a standard choice of basis." }
+{ $notelist
+    { $equiv-word-note "matrix" outer-product }
+    $2d-only-note
+    { $matrix-scalar-note * }
+}
+{ $examples
+    { $unchecked-example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { 1 2 }
+    { 3 4 }
+} {
+    { 0 5 }
+    { 6 7 }
+} kronecker-product ."
+"{
+    { 0 5 0 10 }
+    { 6 7 12 14 }
+    { 0 15 0 20 }
+    { 18 21 24 28 }
+}" }
+} ;
+
+HELP: outer-product
+{ $values { "u" sequence } { "v" sequence } { "matrix" matrix } }
+{ $description "Computes the " { $url URL" http://  enwp.org/Outer_product" "outer-product product" } " of " { $snippet "u" } " and " { $snippet "v" } "." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+        "{ 5 6 7 } { 1 2 3 } outer-product ."
+        "{ { 5 10 15 } { 6 12 18 } { 7 14 21 } }" }
+} ;
+
+HELP: rank
+{ $values { "matrix" matrix } { "rank" rank-kind } }
+{ $contract "The " { $emphasis "rank" } " of a " { $link matrix } " is how its number of linearly independent columns compare to the maximal number of linearly independent columns for a matrix with the same dimension." }
+{ $notes "See " { $url "https://en.wikipedia.org/wiki/Rank_(linear_algebra)" } " for more information." } ;
+
+HELP: nullity
+{ $values { "matrix" matrix } { "nullity" rank-kind } }
+;
+
+HELP: determinant
+{ $values { "matrix" square-matrix } { "determinant" number } }
+{ $contract "Compute the determinant of the input matrix. Generally, the determinant of a matrix is a scaling factor of the transformation described by the matrix." }
+{ $notelist
+    $2d-only-note
+    { $matrix-scalar-note max - * }
+}
+{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    {  3  0 -1 }
+    { -3  1  3 }
+    {  2 -5  4 }
+} determinant ."
+        "44"
+    }
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { -8 -8 13 11 10 -5 -14 }
+    { 3 -11 -8 3 -7 -3 4 }
+    { 10 4 -5 3 0 -6 -12 }
+    { -14 0 -3 -8 10 0 10 }
+    { 3 -6 1 -10 -9 10 0 }
+    { 5 -12 -14 6 5 -1 -7 }
+    { -9 -14 -8 5 2 2 -2 }
+} determinant ."
+        "-103488155"
+    }
+} ;
+
+HELP: 1/det
+{ $values { "matrix" square-matrix } { "1/det" number } }
+{ $description "Find the inverse (" { $link recip } ") of the " { $link determinant } " of the input matrix." }
+{ $notelist
+    $2d-only-note
+    { $matrix-scalar-note determinant recip }
+}
+{ $errors
+    { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "."
+    $nl
+    { $link division-by-zero } " if the " { $link determinant } " of the input matrix is " { $snippet "0" } "."
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { 0 10 -12 4 }
+    { -9 6 -11 9 }
+    { -5 -10 0 2 }
+    { -7 -11 10 11 }
+} 1/det ."
+        "-1/9086"
+    }
+} ;
+
+HELP: m*1/det
+{ $values { "matrix" square-matrix } { "matrix'" square-matrix } }
+{ $description "Multiply the input matrix by the inverse (" { $link recip } ") of its " { $link determinant } "." }
+{ $notelist
+    { "This word is used to implement " { $link recip } " for " { $link square-matrix } "." }
+    $2d-only-note
+    { $matrix-scalar-note determinant recip }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { -14 0 -13 7 }
+    { -4 11 7 -12 }
+    { -3 2 9 -14 }
+    { 3 -5 10 -2 }
+} m*1/det ."
+"{
+    { 7/6855 0 13/13710 -7/13710 }
+    { 2/6855 -11/13710 -7/13710 2/2285 }
+    { 1/4570 -1/6855 -3/4570 7/6855 }
+    { -1/4570 1/2742 -1/1371 1/6855 }
+}"
+    }
+}
+;
+
+HELP: >minors
+{ $values { "matrix" square-matrix } { "matrix'" square-matrix } }
+{ $description "Calculate the " { $emphasis "matrix of minors" } " of the input matrix. See " { $url URL" https://en.wikipedia.org/wiki/Minor_(linear_algebra)" "minor on Wikipedia" } "." }
+{ $notelist
+    $keep-shape-note
+    $2d-only-note
+    { $matrix-scalar-note determinant }
+}
+{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." }
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { -8 0 7 -11 }
+    { 15 0 -3 -11 }
+    { 1 -10 -4 6 }
+    { 11 -15 3 -15 }
+} >minors ."
+"{
+    { 1710 -130 2555 -1635 }
+    { -690 -286 -2965 1385 }
+    { 1650 -754 3795 -1215 }
+    { 1100 416 2530 -810 }
+}"
+    }
+} ;
+
+HELP: >cofactors
+{ $values { "matrix" matrix } { "matrix'" matrix } }
+{ $description "Calculate the " { $emphasis "matrix of cofactors" } " of the input matrix. See " { $url URL" https://en.wikipedia.org/wiki/Minor_(linear_algebra)#Inverse_of_a_matrix" "matrix of cofactors on Wikipedia" } ". Alternating elements of the input matrix have their signs inverted." $nl "On odd rows, the even elements have their signs inverted. On even rows, odd elements have their signs inverted." }
+{ $notelist
+    $keep-shape-note
+    $2d-only-note
+    { $matrix-scalar-note neg }
+}
+{ $examples
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { 8 0 7 11 }
+    { 15 0 3 11 }
+    { 1 10 4 6 }
+    { 11 15 3 15 }
+} >cofactors ."
+"{
+    { 8 0 7 -11 }
+    { -15 0 -3 11 }
+    { 1 -10 4 -6 }
+    { -11 15 -3 15 }
+}"
+    }
+    { $example
+        "USING: math.matrices.extras prettyprint ;"
+"{
+    { -8 0 7 -11 }
+    { 15 0 -3 -11 }
+    { 1 -10 -4 6 }
+    { 11 -15 3 -15 }
+} >cofactors ."
+"{
+    { -8 0 7 11 }
+    { -15 0 3 -11 }
+    { 1 10 -4 -6 }
+    { -11 -15 -3 -15 }
+}"
+    }
+} ;
+
+HELP: multiplicative-inverse
+{ $values { "x" matrix } { "y" matrix } }
+{ $description "Calculate the multiplicative inverse of the input." $nl "If the input is a " { $link square-matrix } ", this is done by multiplying the " { $link transpose } " of the " { $link2 >cofactors "cofactors" } " of the " { $link2 >minors "minors" } " of the input matrix by the " { $link2 1/det "inverse of the determinant" } " of the input matrix."  }
+{ $notelist
+    $keep-shape-note
+    $2d-only-note
+    { $matrix-scalar-note determinant >cofactors 1/det }
+}
+{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." } ;
+
+
+HELP: covariance-matrix-ddof
+{ $values { "matrix" matrix } { "ddof" object } { "cov" matrix } }
+;
+HELP: covariance-matrix
+{ $values { "matrix" matrix } { "cov" matrix } }
+;
+HELP: sample-covariance-matrix
+{ $values { "matrix" matrix } { "cov" matrix } }
+;
+HELP: population-covariance-matrix
+{ $values { "matrix" matrix } { "cov" matrix } }
+;

extra/math/matrices/extras/extras-tests.factor

@@ -0,0 +1,362 @@
+! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
+USING: arrays combinators.short-circuit grouping kernel math
+math.matrices math.matrices.extras math.matrices.extras.private
+math.statistics math.vectors sequences sequences.deep sets
+tools.test ;
+IN: math.matrices.extras
+
+<PRIVATE
+: call-eq? ( obj quots -- ? )
+    [ call( x -- x ) ] with map all-eq? ; !  inline
+PRIVATE>
+
+{ {
+    { 4181 6765 }
+    { 6765 10946 }
+} } [
+  { { 0 1 } { 1 1 } } 20 m^n
+] unit-test
+
+[ { { 0 1 } { 1 1 } } -20 m^n ] [ negative-power-matrix? ] must-fail-with
+[ { { 0 1 } { 1 1 } } -8 m^n ] [ negative-power-matrix? ] must-fail-with
+
+{ { 1 -2 3 -4 } } [ { 1 2 3 4 } t alternating-sign ] unit-test
+{ { -1 2 -3 4 } } [ { 1 2 3 4 } f alternating-sign ] unit-test
+
+
+{ t } [ 50 <box-matrix> dup transpose = ] unit-test
+{ t } [ 50 <box-matrix> dup                  anti-transpose = ] unit-test
+{ f } [ 4 <box-matrix> zero-matrix? ] unit-test
+
+{ t } [ 2 4 15 <random-integer-matrix> mabs {
+    [ flatten [ 15 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 2 = ]
+    [ first length 4 = ]
+} 1&& ] unit-test
+
+{ t } [ 4 4 -45 <random-integer-matrix> mabs {
+    [ flatten [ 45 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 4 = ]
+    [ first length 4 = ]
+} 1&& ] unit-test
+
+{ t } [ 2 2 1 <random-integer-matrix> mabs {
+    [ flatten [ 1 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 2 = ]
+    [ first length 2 = ]
+} 1&& ] unit-test
+
+{ t } [ 2 4 .89 <random-unit-matrix> mabs {
+    [ flatten [ .89 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 2 = ]
+    [ first length 4 = ]
+} 1&& ] unit-test
+
+{ t } [ 2 4 -45.89 <random-unit-matrix> mabs {
+    [ flatten [ 45.89 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 2 = ]
+    [ first length 4 = ]
+} 1&& ] unit-test
+
+{ t } [ 4 4 .89 <random-unit-matrix> mabs {
+    [ flatten [ .89 <= ] all? ]
+    [ regular-matrix? ]
+    [ length 4 = ]
+    [ first length 4 = ]
+} 1&& ] unit-test
+
+{ {
+    { 1   1/2 1/3 1/4 }
+    { 1/2 1/3 1/4 1/5 }
+    { 1/3 1/4 1/5 1/6 }
+} } [ 3 4 <hilbert-matrix> ] unit-test
+
+{ {
+    { 1 2 3 4 }
+    { 2 1 2 3 }
+    { 3 2 1 2 }
+    { 4 3 2 1 }
+} } [ 4 <toeplitz-matrix> ] unit-test
+
+{ {
+    { 1 2 3 4 }
+    { 2 3 4 0 }
+    { 3 4 0 0 }
+    { 4 0 0 0 } }
+} [ 4 <hankel-matrix> ] unit-test
+
+{ {
+    { 1 1 1 }
+    { 4 2 1 }
+    { 9 3 1 }
+    { 25 5 1 } }
+} [
+    { 1 2 3 5 } 3 <vandermonde-matrix>
+] unit-test
+
+{ {
+    { 0 5 0 10 }
+    { 6 7 12 14 }
+    { 0 15 0 20 }
+    { 18 21 24 28 }
+} } [ {
+    { 1 2 }
+    { 3 4 }
+} {
+    { 0 5 }
+    { 6 7 }
+} kronecker-product ] unit-test
+
+{ {
+    { 1  1  1  1 }
+    { 1 -1  1 -1 }
+    { 1  1 -1 -1 }
+    { 1 -1 -1  1 }
+} } [ {
+    { 1  1 }
+    { 1 -1 }
+} dup kronecker-product ] unit-test
+
+{ {
+    { 1 1 1 1 1 1 1 1 }
+    { 1 -1 1 -1 1 -1 1 -1 }
+    { 1 1 -1 -1 1 1 -1 -1 }
+    { 1 -1 -1 1 1 -1 -1 1 }
+    { 1 1 1 1 -1 -1 -1 -1 }
+    { 1 -1 1 -1 -1 1 -1 1 }
+    { 1 1 -1 -1 -1 -1 1 1 }
+    { 1 -1 -1 1 -1 1 1 -1 }
+} } [ {
+    { 1 1 }
+    { 1 -1 }
+} dup dup kronecker-product kronecker-product ] unit-test
+
+{ {
+    { 1 1 1 1 1 1 1 1 }
+    { 1 -1 1 -1 1 -1 1 -1 }
+    { 1 1 -1 -1 1 1 -1 -1 }
+    { 1 -1 -1 1 1 -1 -1 1 }
+    { 1 1 1 1 -1 -1 -1 -1 }
+    { 1 -1 1 -1 -1 1 -1 1 }
+    { 1 1 -1 -1 -1 -1 1 1 }
+    { 1 -1 -1 1 -1 1 1 -1 }
+} } [ {
+    { 1 1 }
+    { 1 -1 }
+} dup dup kronecker-product swap kronecker-product ] unit-test
+
+
+! kronecker-product is not generally commutative, make sure we have the right order
+{ {
+    { 1 2 3 4 5 1 2 3 4 5 }
+    { 6 7 8 9 10 6 7 8 9 10 }
+    { 1 2 3 4 5 -1 -2 -3 -4 -5 }
+    { 6 7 8 9 10 -6 -7 -8 -9 -10 }
+} } [ {
+    { 1 1 }
+    { 1 -1 }
+} {
+    { 1 2 3 4 5 }
+    { 6 7 8 9 10 }
+} kronecker-product ] unit-test
+
+{ {
+    { 1 1 2 2 3 3 4 4 5 5 }
+    { 1 -1 2 -2 3 -3 4 -4 5 -5 }
+    { 6 6 7 7 8 8 9 9 10 10 }
+    { 6 -6 7 -7 8 -8 9 -9 10 -10 }
+} } [ {
+    { 1 1 }
+    { 1 -1 }
+} {
+    { 1 2 3 4 5 }
+    { 6 7 8 9 10 }
+} swap kronecker-product ] unit-test
+
+{ {
+    { 5 10 15 }
+    { 6 12 18 }
+    { 7 14 21 }
+} } [
+    { 5 6 7 }
+    { 1 2 3 }
+    outer-product
+] unit-test
+
+
+CONSTANT: test-points {
+    { 80  27  89 } { 80  27  88 } { 75  25  90 }
+    { 62  24  87 } { 62  22  87 } { 62  23  87 }
+    { 62  24  93 } { 62  24  93 } { 58  23  87 }
+    { 58  18  80 } { 58  18  89 } { 58  17  88 }
+    { 58  18  82 } { 58  19  93 } { 50  18  89 }
+    { 50  18  86 } { 50  19  72 } { 50  19  79 }
+    { 50  20  80 } { 56  20  82 } { 70  20  91 }
+}
+
+{ {
+    { 84+2/35 22+23/35 24+4/7 }
+    { 22+23/35 9+104/105 6+87/140 }
+    { 24+4/7 6+87/140 28+5/7 }
+} } [ test-points sample-covariance-matrix ] unit-test
+
+{ {
+    { 80+8/147 21+85/147 23+59/147 }
+    { 21+85/147 9+227/441 6+15/49 }
+    { 23+59/147 6+15/49 27+17/49 }
+} } [ test-points covariance-matrix ] unit-test
+
+{
+    {
+        { 80+8/147 21+85/147 23+59/147 }
+        { 21+85/147 9+227/441 6+15/49 }
+        { 23+59/147 6+15/49 27+17/49 }
+    }
+} [
+    test-points population-covariance-matrix
+] unit-test
+
+{ t } [ { { 1 } }
+    { [ drop 1 ] [ (1determinant) ] [ 1 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ 0 } [ { { 0 } } determinant ] unit-test
+
+{ t } [ {
+    { 4 6 } ! order is significant
+    { 3 8 }
+} { [ drop 14 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 3 8 }
+    { 4 6 }
+} { [ drop -14 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 2 5 }
+    { 1 -3 }
+} { [ drop -11 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 1 -3 }
+    { 2 5 }
+} { [ drop 11 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 3 0 -1 }
+    { 2 -5 4 }
+    { -3 1 3 }
+} { [ drop -44 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 3 0 -1 }
+    { -3 1 3 }
+    { 2 -5 4 }
+} { [ drop 44 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 2 -3 1 }
+    { 4 2 -1 }
+    { -5 3 -2 }
+} { [ drop -19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 2 -3 1 }
+    { -5 3 -2 }
+    { 4 2 -1 }
+} { [ drop 19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 4 2 -1 }
+    { 2 -3 1 }
+    { -5 3 -2 }
+} { [ drop 19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 5 1 -2 }
+    { -1 0 4 }
+    { 2 -3 3 }
+} { [ drop 65 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 6 1 1 }
+    { 4 -2 5 }
+    { 2 8 7 }
+} { [ drop -306 ] [ (3determinant) ]  [ 3 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { -5  4 -3  2 }
+    { -2  1  0 -1 }
+    { -2 -3 -4 -5  }
+    {  0  2  0  4 }
+} { [ drop -24 ] [ 4 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ t } [ {
+    { 2 4 2 2 }
+    { 5 1 -6 10 }
+    { 4 3 -1 7 }
+    { 9 8 7 3 }
+} { [ drop 272 ] [ 4 swap (ndeterminant) ] [ determinant ] }
+    call-eq?
+] unit-test
+
+{ {
+    { 2 2 2 }
+    { -2 3 3 }
+    { 0 -10 0 }
+} } [ {
+    { 3 0 2 }
+    { 2 0 -2 }
+    { 0 1 1 }
+} >minors ] unit-test
+
+! i think this unit test is wrong
+! { {
+!     { 1 -6 -13 }
+!     { 0 0 0 }
+!     { 1 -6 -13 }
+! } } [ {
+!     { 1 2 1 }
+!     { 6 -1 0 }
+!     { 1 -2 -1 }
+! } >minors ] unit-test
+
+{ {
+    { 1 6 -13 }
+    { 0 0 0 }
+    { 1 6 -13 }
+} } [ {
+    { 1 -6 -13 }
+    { 0 0 0 }
+    { 1 -6 -13 }
+} >cofactors ] unit-test

extra/math/matrices/extras/extras.factor

@@ -0,0 +1,338 @@
+! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
+! See http://factorcode.org/license.txt for BSD license.
+USING: accessors arrays combinators formatting fry kernel locals
+math math.bits math.functions math.matrices
+math.matrices.private math.order math.statistics math.text.english
+math.vectors random sequences sequences.deep summary ;
+IN: math.matrices.extras
+
+! this is a questionable implementation
+SINGLETONS:      +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
+UNION: rank-kind +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
+
+ERROR: negative-power-matrix
+    { m matrix } { n integer } ;
+ERROR: non-square-determinant
+    { m integer }  { n integer } ;
+ERROR: undefined-inverse
+    { m integer }  { n integer } { r rank-kind initial: +uncalculated-rank+ } ;
+
+M: negative-power-matrix summary
+    n>> dup ordinal-suffix "%s%s power of a matrix is undefined" sprintf ;
+M: non-square-determinant summary
+    [ m>> ] [ n>> ] bi "non-square %s x %s matrix has no determinant" sprintf ;
+M: undefined-inverse summary
+    [ m>> ] [ n>> ] [ r>> name>> ] tri "%s x %s matrix of rank %s has no inverse" sprintf ;
+
+<PRIVATE
+DEFER: alternating-sign
+: finish-randomizing-matrix ( matrix -- matrix' )
+    [ f alternating-sign randomize ] map randomize ; inline
+PRIVATE>
+
+: <random-integer-matrix> ( m n max -- matrix )
+    '[ _ _ 1 + random-integers ] replicate
+    finish-randomizing-matrix ; inline
+
+: <random-unit-matrix> ( m n max -- matrix )
+    '[ _ random-units [ _ * ] map ] replicate
+    finish-randomizing-matrix ; inline
+
+<PRIVATE
+: (gram-schmidt) ( v seq -- newseq )
+    [ dupd proj v- ] each ;
+PRIVATE>
+
+: gram-schmidt ( matrix -- orthogonal )
+    [ V{ } clone [ over (gram-schmidt) suffix! ] reduce ] keep like ;
+
+: gram-schmidt-normalize ( matrix -- orthonormal )
+    gram-schmidt [ normalize ] map ; inline
+
+: kronecker-product ( m1 m2 -- m )
+    '[ [ _ n*m  ] map ] map stitch stitch ;
+
+: outer-product ( u v -- matrix )
+    '[ _ n*v ] map ;
+
+! Special matrix constructors follow
+: <hankel-matrix> ( n -- matrix )
+  [ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
+
+: <hilbert-matrix> ( m n -- matrix )
+    [ <iota> ] bi@ [ + 1 + recip ] cartesian-map ;
+
+: <toeplitz-matrix> ( n -- matrix )
+    <iota> dup [ - abs 1 + ] cartesian-map ;
+
+: <box-matrix> ( r -- matrix )
+    2 * 1 + dup '[ _ 1 <array> ] replicate ;
+
+: <vandermonde-matrix> ( u n -- matrix )
+    <iota> [ v^n ] with map reverse flip ;
+
+! Transformation matrices
+:: <rotation-matrix3> ( axis theta -- matrix )
+    theta cos :> c
+    theta sin :> s
+    axis first3 :> ( x y z )
+    x sq 1.0 x sq - c * +    x y * 1.0 c - * z s * -  x z * 1.0 c - * y s * + 3array
+    x y * 1.0 c - * z s * +  y sq 1.0 y sq - c * +    y z * 1.0 c - * x s * - 3array
+    x z * 1.0 c - * y s * -  y z * 1.0 c - * x s * +  z sq 1.0 z sq - c * +   3array
+    3array ;
+
+:: <rotation-matrix4> ( axis theta -- matrix )
+    theta cos :> c
+    theta sin :> s
+    axis first3 :> ( x y z )
+    x sq 1.0 x sq - c * +    x y * 1.0 c - * z s * -  x z * 1.0 c - * y s * +  0 4array
+    x y * 1.0 c - * z s * +  y sq 1.0 y sq - c * +    y z * 1.0 c - * x s * -  0 4array
+    x z * 1.0 c - * y s * -  y z * 1.0 c - * x s * +  z sq 1.0 z sq - c * +    0 4array
+    { 0.0 0.0 0.0 1.0 } 4array ;
+
+:: <translation-matrix4> ( offset -- matrix )
+    offset first3 :> ( x y z )
+    {
+        { 1.0 0.0 0.0 x   }
+        { 0.0 1.0 0.0 y   }
+        { 0.0 0.0 1.0 z   }
+        { 0.0 0.0 0.0 1.0 }
+    } ;
+
+<PRIVATE
+GENERIC: >scale-factors ( object -- x y z )
+M: number >scale-factors
+    dup dup ;
+M: sequence >scale-factors
+    first3 ;
+PRIVATE>
+
+:: <scale-matrix3> ( factors -- matrix )
+    factors >scale-factors :> ( x y z )
+    {
+        { x   0.0 0.0 }
+        { 0.0 y   0.0 }
+        { 0.0 0.0 z   }
+    } ;
+
+:: <scale-matrix4> ( factors -- matrix )
+    factors >scale-factors :> ( x y z )
+    {
+        { x   0.0 0.0 0.0 }
+        { 0.0 y   0.0 0.0 }
+        { 0.0 0.0 z   0.0 }
+        { 0.0 0.0 0.0 1.0 }
+    } ;
+
+: <ortho-matrix4> ( factors -- matrix )
+    [ recip ] map <scale-matrix4> ;
+
+:: <frustum-matrix4> ( xy-dim near far -- matrix )
+    xy-dim first2 :> ( x y )
+    near x /f :> xf
+    near y /f :> yf
+    near far + near far - /f :> zf
+    2 near far * * near far - /f :> wf
+
+    {
+        { xf  0.0  0.0 0.0 }
+        { 0.0 yf   0.0 0.0 }
+        { 0.0 0.0  zf  wf  }
+        { 0.0 0.0 -1.0 0.0 }
+    } ;
+
+:: <skew-matrix4> ( theta -- matrix )
+    theta tan :> zf
+    {
+        { 1.0 0.0 0.0 0.0 }
+        { 0.0 1.0 0.0 0.0 }
+        { 0.0 zf  1.0 0.0 }
+        { 0.0 0.0 0.0 1.0 }
+    } ;
+
+! a simpler verison of this, like matrix-map -except, but map-index, should be possible
+: cartesian-matrix-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
+    [ [ first length <cartesian-square-indices> ] keep ] dip
+    '[ _ @ ] matrix-map ; inline
+
+: cartesian-column-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
+    [ cols first2 ] prepose cartesian-matrix-map ; inline
+
+! -------------------------------------------------
+! numerical analysis of matrices follows
+<PRIVATE
+
+: square-rank ( square-matrix -- rank ) ;
+: nonsquare-rank ( matrix -- rank ) ;
+PRIVATE>
+
+GENERIC: rank ( matrix -- rank )
+M: zero-matrix rank
+    drop +zero-rank+ ;
+
+M: square-matrix rank
+    square-rank ;
+
+M: matrix rank
+    nonsquare-rank ;
+
+GENERIC: nullity ( matrix -- nullity )
+
+
+! implementation details of determinant and inverse
+<PRIVATE
+: alternating-sign ( seq odd-elts? -- seq' )
+    '[ even? _ = [ neg ] unless ] map-index ;
+
+! the determinant of a 1x1 matrix is the value itself
+! this works for any-dimensional matrices too
+: (1determinant) ( matrix -- 1det ) flatten first ; inline
+
+! optimized to find the determinant of a 2x2 matrix
+: (2determinant) ( matrix -- 2det )
+    ! multiply the diagonals and subtract
+    [ main-diagonal ] [ anti-diagonal ] bi [ first2 * ] bi@ - ; inline
+
+! optimized for 3x3
+! https://www.mathsisfun.com/algebra/matrix-determinant.html
+:: (3determinant) ( matrix-seq -- 3det )
+    ! first 3 elements of row 1
+    matrix-seq first first3 :> ( a b c )
+    ! last 2 rows, transposed to make the next step easier
+    matrix-seq rest transpose
+    ! get the lower sub-matrices in reverse order of a b c columns
+    [ rest ] [ [ first ] [ third ] bi 2array ] [ 1 head* ] tri 3array
+    ! find determinants
+    [ (2determinant) ] map
+    ! negate odd elements of a b c and multiply by the new determinants
+    { a b c } t alternating-sign v*
+    ! sum the resulting sequence
+    sum ;
+
+DEFER: (ndeterminant)
+: make-determinants ( n matrix -- seq )
+    <repetition> [
+        cols-except [ length ] keep (ndeterminant) ! recurses here
+    ] map-index ;
+
+DEFER: (determinant)
+! generalized to 4 and higher
+: (ndeterminant) ( n matrix -- ndet )
+    ! TODO? recurse for n < 3
+    over 4 < [ (determinant) ] [
+        [ nip first t alternating-sign ] [ rest make-determinants ] 2bi
+        v* sum
+    ] if ;
+
+! switches on dimensions only
+: (determinant) ( n matrix -- determinant )
+    over {
+        { 1 [ nip (1determinant) ] }
+        { 2 [ nip (2determinant) ] }
+        { 3 [ nip (3determinant) ] }
+        [ drop (ndeterminant) ]
+    } case ;
+PRIVATE>
+
+GENERIC: determinant ( matrix -- determinant )
+M: zero-square-matrix determinant
+    drop 0 ;
+
+M: square-matrix determinant
+    [ length ] keep (determinant) ;
+
+! determinant is undefined for m =/= n, unlike inverse
+M: matrix determinant
+    dimension first2 non-square-determinant ;
+
+: 1/det ( matrix -- 1/det )
+    determinant recip ; inline
+
+! -----------------------------------------------------
+! inverse operations and implementations follow
+ALIAS: multiplicative-inverse recip
+
+! per element, find the determinant of all other elements except the element's row / col
+! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
+: >minors ( matrix -- matrix' )
+    matrix-except-all [ [ determinant ] map ] map ;
+
+! alternately invert values of the matrix (see alternating-sign)
+: >cofactors ( matrix -- matrix' )
+    [ even? alternating-sign ] map-index ;
+
+! multiply a matrix by the inverse of its determinant
+: m*1/det ( matrix -- matrix' )
+    [ 1/det ] keep n*m ; inline
+
+! inverse implementation
+<PRIVATE
+! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
+: (square-inverse) ( square-matrix -- inverted )
+    ! inverse of the determinant of the input matrix
+    [ 1/det ]
+    ! adjugate of the cofactors of the matrix of minors
+    [ >minors >cofactors transpose ]
+    ! adjugate * 1/det
+    bi n*m ;
+
+! TODO
+: (left-inverse) ( matrix -- left-invert )   ;
+: (right-inverse) ( matrix -- right-invert ) ;
+
+! TODO update this when rank works properly
+! only defined for rank(A) = rows(A) OR rank(A) = cols(M)
+! https://en.wikipedia.org/wiki/Invertible_matrix
+: (specialized-inverse) ( rect-matrix -- inverted )
+    dup [ rank ] [ dimension ] bi [ = ] with map {
+        { { t f } [ (left-inverse) ] }
+        { { f t } [ (right-inverse) ] }
+        [ no-case ]
+    } case ;
+PRIVATE>
+
+M: zero-square-matrix recip
+    ; inline
+
+M: square-matrix recip
+    (square-inverse) ; inline
+
+M: zero-matrix recip
+    transpose ; inline ! TODO: error based on rankiness
+
+M: matrix recip
+    (specialized-inverse) ; inline
+
+! TODO: use the faster algorithm: [ determinant zero? ]
+: invertible-matrix? ( matrix -- ? )
+    [ dimension first2 max <identity-matrix> ] keep
+    dup recip m. = ;
+
+: linearly-independent-matrix? ( matrix -- ? ) ;
+
+<PRIVATE
+! this is the original definition of m^n as committed in 2012; it has not been lost
+: (m^n) ( m n -- n )
+    make-bits over first length <identity-matrix>
+    [ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
+PRIVATE>
+
+! A^-1 is the inverse but other negative powers are nonsense
+: m^n ( m n -- n ) {
+        { [ dup -1 = ] [ drop recip ] }
+        { [ dup 0 >= ] [ (m^n) ] }
+        [ negative-power-matrix ]
+    } cond ;
+
+: n^m ( n m -- n ) swap m^n ; inline
+
+: covariance-matrix-ddof ( matrix ddof -- cov )
+    '[ _ cov-ddof ] cartesian-column-map ; inline
+
+: covariance-matrix ( matrix -- cov )
+    0 covariance-matrix-ddof ; inline
+
+: sample-covariance-matrix ( matrix -- cov )
+    1 covariance-matrix-ddof ; inline
+
+: population-covariance-matrix ( matrix -- cov ) 0 covariance-matrix-ddof ; inline

extra/math/matrices/extras/summary.txt

@@ -0,0 +1 @@
+Matrix arithmetic - extra and miscellaneous words

extra/rosetta-code/bitmap/bitmap.factor

@@ -36,7 +36,7 @@ IN: rosetta-code.bitmap
 
 ! The storage functions
 : <raster-image> ( width height -- image )
-    zero-matrix [ drop { 0 0 0 } ] mmap ;
+    <zero-matrix> [ drop { 0 0 0 } ] mmap ;
 : fill-image ( {R,G,B} image -- image )
     swap '[ drop _ ] mmap! ;
 : set-pixel ( {R,G,B} {i,j} image -- ) set-Mi,j ; inline

extra/rosetta-code/conjugate-transpose/conjugate-transpose.factor

@@ -38,4 +38,4 @@ IN: rosetta-code.conjugate-transpose
     dup conj-t [ m. ] [ swap m. ] 2bi = ;
 
 : unitary-matrix? ( matrix -- ? )
-    [ dup conj-t m. ] [ length identity-matrix ] bi = ;
+    [ dup conj-t m. ] [ length <identity-matrix> ] bi = ;

extra/tensors/benchmark/benchmark.factor

@@ -0,0 +1,50 @@
+! Copyright (C) 2019 HMC Clinic.
+! See http://factorcode.org/license.txt for BSD license.
+USING: arrays io kernel locals math prettyprint tensors tools.time ;
+IN: tensors.benchmark
+
+<PRIVATE
+
+:: add-tensors ( trials elems -- time )
+    ! Create the arrays to be added
+    elems 1array naturals dup
+    ! Benchmark!
+    [ trials [ 2dup t+ drop ] times ] benchmark
+    ! Normalize
+    trials / >float
+    nip nip ;
+
+:: matmul-tensors ( trials elems -- time )
+    ! Create the arrays to be multiplied
+    elems elems 2array naturals dup
+    ! Benchmark!
+    [ trials [ 2dup matmul drop ] times ] benchmark
+    ! Normalize
+    trials / >float
+    nip nip ;
+
+:: transpose-tensor ( trials elems -- time )
+    ! Create the array to be transposed
+    elems elems 2array naturals
+    ! benchmark
+    [ trials [ dup transpose drop ] times ] benchmark
+    ! Normalize
+    trials / >float
+    nip ;
+
+PRIVATE>
+
+: run-benchmarks ( -- )
+    "Benchmarking the tensors vocabulary" print
+    "Add two 100 element tensors" print
+    1000000 100 add-tensors .
+    "Add two 100,000 element tensors" print
+    10000 100000 add-tensors .
+    "Multiply two 10x10 matrices" print
+    100000 10 matmul-tensors .
+    "Multiply two 100x100 matrices" print
+    1000 100 matmul-tensors .
+    "Transpose a 10x10 matrix" print
+    10000 10 transpose-tensor .
+    "Transpose a 100x100 matrix" print
+    10 100 transpose-tensor . ;

extra/tensors/tensors.factor

@@ -1,9 +1,12 @@
 ! Copyright (C) 2019 HMC Clinic.
 ! See http://factorcode.org/license.txt for BSD license.
+
 USING: accessors alien.c-types alien.data arrays
 concurrency.combinators grouping kernel locals math.functions
 math.ranges math.statistics math multi-methods quotations sequences
-sequences.private specialized-arrays tensors.tensor-slice typed ;
+sequences.extras sequences.private specialized-arrays
+tensors.tensor-slice typed ;
+
 QUALIFIED-WITH: alien.c-types c
 SPECIALIZED-ARRAY: c:float
 IN: tensors
@@ -47,7 +50,7 @@ PRIVATE>
 ! Construct a one-dimensional tensor with values start, start+step,
 ! ..., stop (inclusive)
 : arange ( a b step -- tensor )
-    <range> [ length 1array ] keep >float-array <tensor> ;
+    <range> [ length >fixnum 1array ] keep >float-array <tensor> ;
 
 ! Construct a tensors with vec { 0 1 2 ... } and reshape to the desired shape
 : naturals ( shape -- tensor )
@@ -98,32 +101,32 @@ PRIVATE>
 ! Add a tensor to either another tensor or a scalar
 multi-methods:GENERIC: t+ ( x y -- tensor )
 METHOD: t+ { tensor tensor } [ + ] t-bop ;
-METHOD: t+ { tensor number } [ + ] curry t-uop ;
+METHOD: t+ { tensor number } >float [ + ] curry t-uop ;
-METHOD: t+ { number tensor } swap [ + ] curry t-uop ;
+METHOD: t+ { number tensor } [ >float ] dip [ + ] with t-uop ;
 
 ! Subtraction between two tensors or a tensor and a scalar
 multi-methods:GENERIC: t- ( x y -- tensor )
 METHOD: t- { tensor tensor } [ - ] t-bop ;
-METHOD: t- { tensor number } [ - ] curry t-uop ;
+METHOD: t- { tensor number } >float [ - ] curry t-uop ;
-METHOD: t- { number tensor } swap [ swap - ] curry t-uop ;
+METHOD: t- { number tensor } [ >float ] dip [ - ] with t-uop ;
 
 ! Multiply a tensor with either another tensor or a scalar
 multi-methods:GENERIC: t* ( x y -- tensor )
 METHOD: t* { tensor tensor } [ * ] t-bop ;
 METHOD: t* { tensor number } [ * ] curry t-uop ;
-METHOD: t* { number tensor } swap [ * ] curry t-uop ;
+METHOD: t* { number tensor } [ >float ] dip [ * ] with t-uop ;
 
 ! Divide two tensors or a tensor and a scalar
 multi-methods:GENERIC: t/ ( x y -- tensor )
 METHOD: t/ { tensor tensor } [ / ] t-bop ;
 METHOD: t/ { tensor number } [ / ] curry t-uop ;
-METHOD: t/ { number tensor } swap [ swap / ] curry t-uop ;
+METHOD: t/ { number tensor } [ / ] with t-uop ;
 
 ! Divide two tensors or a tensor and a scalar
 multi-methods:GENERIC: t% ( x y -- tensor )
 METHOD: t% { tensor tensor } [ mod ] t-bop ;
 METHOD: t% { tensor number } [ mod ] curry t-uop ;
-METHOD: t% { number tensor } swap [ swap mod ] curry t-uop ;
+METHOD: t% { number tensor } [ mod ] with t-uop ;
 
 <PRIVATE
 
@@ -148,20 +151,24 @@ METHOD: t% { number tensor } swap [ swap mod ] curry t-uop ;
 
 ! Perform matrix multiplication muliplying an
 ! mxn matrix with a nxp matrix
-TYPED:: 2d-matmul ( vec1: slice vec2: slice res: slice n: number p: number -- )
+TYPED:: 2d-matmul ( vec1: float-array start1: fixnum
+                      vec2: float-array start2: fixnum
+                      res: float-array start3: fixnum
+                      m: fixnum n: fixnum p: fixnum -- )
     ! For each element in the range, we want to compute the dot product of the
     ! corresponding row and column
-    res
+    m [ :> i
-    [   >fixnum
+        p [ :> j
-        ! Get the row
+            0.0 ! This is the sum
-        [ [ vec1 n ] dip p row ]
+            n [ :> k
-        ! Get the column
+                ! Add to the sum
-        ! [ p mod vec2 swap p every ] bi
+                i n * k + start1 + vec1 nth-unsafe
-        [ p mod f p vec2 <step-slice> ] bi
+                k p * j + start2 + vec2 nth-unsafe
-        ! Take the dot product
+                * +
-        [ * ] [ + ] 2map-reduce
+            ] each-integer
-    ]
+            i p * j + start3 + res set-nth-unsafe
-    map! drop ;
+        ] each-integer
+    ] each-integer ;
 
 PRIVATE>
 
@@ -176,70 +183,52 @@ TYPED:: matmul ( tensor1: tensor tensor2: tensor -- tensor3: tensor )
     tensor1 shape>> unclip-last-slice :> n
     unclip-last-slice :> m :> top-shape
     tensor2 shape>> last :> p
-    top-shape product :> rest
+    top-shape product :> top-prod
 
-    ! Now create the new tensor with { 0 ... m*p-1 } repeating
+    ! Create the shape of the resulting tensor
-    top-shape { m p } append naturals m p * t% :> tensor3
+    top-shape { m p } append
+
+    ! Now create the new float array to store the underlying result
+    dup product c:float (c-array) :> vec3
 
     ! Now update the tensor3 to contain the multiplied matricies
-    rest [0,b)
+    top-prod [
-    [
         :> i
-        ! First make vec1
+        ! Compute vec1 and start1
-        m n * i * dup m n * + tensor1 vec>> <slice>
+        tensor1 vec>> m n * i *
-        ! Now make vec2
+        ! Compute vec2 and start2
-        n p * i * dup n p * + tensor2 vec>> <slice>
+        tensor2 vec>> n p * i *
-        ! Now make the resulting vector
+        ! Compute the result
-        m p * i * dup m p * + tensor3 vec>> <slice>
+        vec3 m p * i *
-        ! Push n and p and multiply the clices
+        ! Push m, n, and p and multiply the arrays
-        n p 2d-matmul
+        m n p 2d-matmul
-        0
+    ] each-integer
-    ] map drop
+    vec3 <tensor> ;
-    tensor3 ;
+
 
 <PRIVATE
-! helper for transpose: gets the turns a shape into a list of things
+! helper for transpose: turns a shape into a list of things
 ! by which to multiply indices to get a full index
 : ind-mults ( shape -- seq )
-    rest-slice <reversed> cum-product { 1 } prepend ;
+    <reversed> 1 swap [ swap [ * ] keep ] map nip ;
 
 ! helper for transpose: given shape, flat index, & mults for the shape, gives nd index
-:: trans-index ( ind shape mults -- seq )
+: transpose-index ( i shape -- seq )
-    ! what we use to divide things
+    <reversed> [ /mod ] map reverse nip ;
-    shape reverse :> S
-    ! accumulator
-    V{ } clone
-    ! loop thru elements & indices of S (mod by elment m)
-    S [| m i |
-        ! we divide by the product of the 1st n elements of S
-        S i head-slice product :> div
-        ! do not mod on the last index
-        i S length 1 - = not :> mod?
-        ! multiply accumulator by mults & sum
-        dup mults [ * ] 2map sum
-        ! subtract from ind & divide
-        ind swap - div /
-        ! mod if necessary
-        mod? [ m mod ] [ ] if
-        ! append to accumulator
-        [ dup ] dip swap push
-    ] each-index
-    reverse ;
 PRIVATE>
 
-! Transpose an n-dimensional tensor
+! Transpose an n-dimensional tensor by flipping the axes
 TYPED:: transpose ( tensor: tensor -- tensor': tensor )
-    ! new shape
+    tensor shape>> :> old-shape
-    tensor shape>> reverse :> newshape
+    tensor vec>> :> vec
-    ! what we multiply by to get indices in the old tensor
+    old-shape reverse :> new-shape
-    tensor shape>> ind-mults :> old-mults
+    ! check that the size is fine
-    ! what we multiply to get indices in new tensor
+    new-shape product vec length assert=
-    newshape ind-mults :> mults
+    old-shape ind-mults reverse :> mults
-    ! new tensor of correct shape
+    ! loop through new tensor
-    newshape naturals dup vec>>
+    new-shape dup product <iota> [
-    [ ! go thru each index
         ! find index in original tensor
-        newshape mults trans-index old-mults [ * ] 2map sum >fixnum
+        old-shape mults [ [ /mod ] dip * ] 2map-sum nip
         ! get that index in original tensor
-        tensor vec>> nth
+        vec nth-unsafe
-    ] map! >>vec ;
+    ] float-array{ } map-as <tensor> ;

extra/trees/trees.factor

@@ -238,10 +238,10 @@ PRIVATE>
     [ root>> (nodepath-at) ] { } make ;
 
 : right-extremity ( node -- node' )
-    [ dup right>> dup ] [ nip ] while drop ;
+    [ dup right>> ] [ nip ] while* ;
 
 : left-extremity ( node -- node' )
-    [ dup left>> dup ] [ nip ] while drop ;
+    [ dup left>> ] [ nip ] while* ;
 
 : lower-node-in-child? ( key node -- ? )
     [ nip left>> ] [ key>> = ] 2bi and ;

misc/vim/syntax/factor.vim

@@ -53,17 +53,17 @@ syn match factorCallQuotation /\<call(\s\+\(\S*\s\+\)*--\(\s\+\S*\)*\s\+)\>/ con
 syn match factorExecute /\<execute(\s\+\(\S*\s\+\)*--\(\s\+\S*\)*\s\+)\>/ contained contains=factorStackEffect
 syn keyword factorCallNextMethod call-next-method
 
-syn keyword factorKeyword (clone) -roll -rot -rotd 2bi 2bi* 2bi@ 2curry 2dip 2drop 2dup 2keep 2keepd 2nip 2nipd 2over 2tri 2tri* 2tri@ 2with 3bi 3curry 3dip 3drop 3dup 3keep 3nip 3nipd 3tri 4dip 4drop 4dup 4keep 4nip 5drop 5nip <wrapper> = >boolean ? ?if and assert assert= assert? bi bi* bi-curry bi-curry* bi-curry@ bi@ boa boolean boolean? both? build call callstack callstack>array callstack? clear clone compose composed composed? curried curried? curry die dip do drop dup dupd either? eq? equal? execute get-callstack get-datastack get-retainstack hashcode hashcode* identity-hashcode identity-tuple identity-tuple? if if* keep keepd keepdd loop most new nip nipd not null object or over overd pick pickd prepose reach roll rot rotd same? spin swap swapd throw tri tri* tri-curry tri-curry* tri-curry@ tri@ tuck tuple tuple? unless unless* until when when* while with wrapper wrapper? xor
+syn keyword factorKeyword (clone) -roll -rot -rotd 2bi 2bi* 2bi@ 2curry 2dip 2drop 2dup 2keep 2keepd 2nip 2nipd 2over 2tri 2tri* 2tri@ 2with 3bi 3curry 3dip 3drop 3dup 3keep 3nip 3nipd 3tri 4dip 4drop 4dup 4keep 4nip 5drop 5nip <wrapper> = >boolean ? ?if and assert assert= assert? bi bi* bi-curry bi-curry* bi-curry@ bi@ boa boolean boolean? both? build call callstack callstack>array callstack? clear clone compose composed composed? curried curried? curry die dip do drop dup dupd either? eq? equal? execute get-callstack get-datastack get-retainstack hashcode hashcode* identity-hashcode identity-tuple identity-tuple? if if* keep keepd keepdd loop most new nip nipd not null object or over overd pick pickd prepose reach roll rot rotd same? spin swap swapd throw tri tri* tri-curry tri-curry* tri-curry@ tri@ tuck tuple tuple? unless unless* until when when* while while* with wrapper wrapper? xor
-syn keyword factorKeyword 2cache <enumerated> >alist ?at ?of assoc assoc-all? assoc-any? assoc-clone-like assoc-combine assoc-diff assoc-diff! assoc-differ assoc-each assoc-empty? assoc-filter assoc-filter! assoc-filter-as assoc-find assoc-hashcode assoc-intersect assoc-like assoc-map assoc-map-as assoc-partition assoc-refine assoc-reject assoc-reject! assoc-reject-as assoc-size assoc-stack assoc-subset? assoc-union assoc-union! assoc-union-as assoc= assoc>map assoc? at at* at+ cache change-at clear-assoc collect-by delete-at delete-at* enumerated enumerated? extract-keys harvest-keys harvest-values inc-at key? keys map>alist map>assoc maybe-set-at new-assoc of push-at rename-at set-at sift-keys sift-values substitute unzip value-at value-at* value? values zip zip-as zip-index zip-index-as
+syn keyword factorKeyword 2cache <enumerated> >alist ?at ?delete-at ?of assoc assoc-all? assoc-any? assoc-clone-like assoc-combine assoc-diff assoc-diff! assoc-differ assoc-each assoc-empty? assoc-filter assoc-filter! assoc-filter-as assoc-find assoc-hashcode assoc-intersect assoc-like assoc-map assoc-map-as assoc-partition assoc-refine assoc-reject assoc-reject! assoc-reject-as assoc-size assoc-stack assoc-subset? assoc-union assoc-union! assoc-union-as assoc= assoc>map assoc? at at* at+ cache change-at clear-assoc collect-by delete-at delete-at* enumerated enumerated? extract-keys harvest-keys harvest-values inc-at key? keys map>alist map>assoc maybe-set-at new-assoc of push-at rename-at set-at sift-keys sift-values substitute unzip value-at value-at* value? values zip zip-as zip-index zip-index-as
 syn keyword factorKeyword 2cleave 2cleave>quot 3cleave 3cleave>quot 4cleave 4cleave>quot alist>quot call-effect case case-find case>quot cleave cleave>quot cond cond>quot deep-spread>quot execute-effect linear-case-quot no-case no-case? no-cond no-cond? recursive-hashcode shallow-spread>quot spread to-fixed-point wrong-values wrong-values?
 syn keyword factorKeyword (all-integers?) iterate-upto find-upto * + - / /f /i /mod 2/ 2^ < <= <fp-nan> > >= >bignum >fixnum >float >fraction >integer >rect ?1+ abs align all-integers? bignum bignum? bit? bitand bitnot bitor bits>double bits>float bitxor complex complex? denominator double>bits each-integer even? find-integer find-last-integer fixnum fixnum? float float>bits float? fp-bitwise= fp-infinity? fp-nan-payload fp-nan? fp-qnan? fp-sign fp-snan? fp-special? gcd if-zero imaginary-part integer integer>fixnum integer>fixnum-strict integer? log2 log2-expects-positive log2-expects-positive? mod neg neg? next-float next-power-of-2 number number= number? numerator odd? power-of-2? prev-float ratio ratio? rational rational? real real-part real? recip rect> rem sgn shift simple-gcd sq times u< u<= u> u>= unless-zero unordered? when-zero zero?
 syn keyword factorKeyword 1sequence 2all? 2each 2each-from 2map 2map-as 2map-reduce 2reduce 2selector 2sequence 3append 3append-as 3each 3map 3map-as 3sequence 4sequence <iota> <repetition> <reversed> <slice> ?first ?last ?nth ?second ?set-nth accumulate accumulate! accumulate* accumulate*! accumulate*-as accumulate-as all? any? append append! append-as assert-sequence assert-sequence= assert-sequence? binary-reduce bounds-check bounds-check? bounds-error bounds-error? but-last but-last-slice cartesian-each cartesian-map cartesian-product change-nth check-slice clone-like collapse-slice collector collector-as collector-for collector-for-as concat concat-as copy count cut cut* cut-slice delete-all delete-slice drop-prefix each each-from each-index empty? exchange filter filter! filter-as find find-from find-index find-index-from find-last find-last-from first first2 first3 first4 flip follow fourth glue halves harvest head head* head-slice head-slice* head? if-empty immutable immutable-sequence immutable-sequence? immutable? index index-from indices infimum infimum-by insert-nth interleave iota iota? join join-as last last-index last-index-from length lengthen like longer longer? longest map map! map-as map-find map-find-last map-index map-index-as map-integers map-reduce map-sum max-length member-eq? member? midpoint@ min-length mismatch move new-like new-resizable new-sequence non-negative-integer-expected non-negative-integer-expected? none? nth nths pad-head pad-tail padding partition pop pop* prefix prepend prepend-as produce produce-as product push push-all push-either push-if reduce reduce-index reject reject! reject-as remove remove! remove-eq remove-eq! remove-nth remove-nth! repetition repetition? replace-slice replicate replicate-as rest rest-slice reverse reverse! reversed reversed? second selector selector-as sequence sequence-hashcode sequence= sequence? set-first set-fourth set-last set-length set-nth set-second set-third short shorten shorter shorter? shortest sift slice slice-error slice-error? slice? snip snip-slice subseq subseq-as subseq-start subseq-start-from subseq? suffix suffix! sum sum-lengths supremum supremum-by surround tail tail* tail-slice tail-slice* tail? third trim trim-head trim-head-slice trim-slice trim-tail trim-tail-slice unclip unclip-last unclip-last-slice unclip-slice unless-empty virtual-exemplar virtual-sequence virtual-sequence? virtual@ when-empty
 syn keyword factorKeyword +@ change change-global counter dec get get-global get-namestack global inc init-namespaces initialize namespace off on set set-global set-namestack toggle with-global with-scope with-variable with-variable-off with-variable-on with-variables
 syn keyword factorKeyword 1array 2array 3array 4array <array> >array array array? pair pair? resize-array
-syn keyword factorKeyword (each-stream-block) (each-stream-block-slice) (stream-contents-by-block) (stream-contents-by-element) (stream-contents-by-length) (stream-contents-by-length-or-block) +byte+ +character+ bad-seek-type bad-seek-type? bl contents each-block each-block-size each-block-slice each-line each-morsel each-stream-block each-stream-block-slice each-stream-line error-stream flush input-stream input-stream? invalid-read-buffer invalid-read-buffer? lines nl output-stream output-stream? print read read-into read-partial read-partial-into read-until read1 readln seek-absolute seek-absolute? seek-end seek-end? seek-input seek-output seek-relative seek-relative? stream-bl stream-contents stream-contents* stream-copy stream-copy* stream-element-type stream-flush stream-length stream-lines stream-nl stream-print stream-read stream-read-into stream-read-partial stream-read-partial-into stream-read-partial-unsafe stream-read-unsafe stream-read-until stream-read1 stream-readln stream-seek stream-seekable? stream-tell stream-write stream-write1 tell-input tell-output with-error-stream with-error-stream* with-error>output with-input-output+error-streams with-input-output+error-streams* with-input-stream with-input-stream* with-output+error-stream with-output+error-stream* with-output-stream with-output-stream* with-output>error with-streams with-streams* write write1
+syn keyword factorKeyword (each-stream-block) (each-stream-block-slice) (stream-contents-by-block) (stream-contents-by-element) (stream-contents-by-length) (stream-contents-by-length-or-block) +byte+ +character+ bad-seek-type bad-seek-type? bl contents each-block each-block-size each-block-slice each-line each-stream-block each-stream-block-slice each-stream-line error-stream flush input-stream input-stream? invalid-read-buffer invalid-read-buffer? lines nl output-stream output-stream? print read read-into read-partial read-partial-into read-until read1 readln seek-absolute seek-absolute? seek-end seek-end? seek-input seek-output seek-relative seek-relative? stream-bl stream-contents stream-contents* stream-copy stream-copy* stream-element-type stream-flush stream-length stream-lines stream-nl stream-print stream-read stream-read-into stream-read-partial stream-read-partial-into stream-read-partial-unsafe stream-read-unsafe stream-read-until stream-read1 stream-readln stream-seek stream-seekable? stream-tell stream-write stream-write1 tell-input tell-output with-error-stream with-error-stream* with-error>output with-input-output+error-streams with-input-output+error-streams* with-input-stream with-input-stream* with-output+error-stream with-output+error-stream* with-output-stream with-output-stream* with-output>error with-streams with-streams* write write1
 syn keyword factorKeyword 1string <string> >string resize-string string string?
 syn keyword factorKeyword 1vector <vector> >vector ?push vector vector?
-syn keyword factorKeyword <condition> <continuation> <restart> attempt-all attempt-all-error attempt-all-error? callback-error-hook callcc0 callcc1 cleanup compute-restarts condition condition? continuation continuation? continue continue-restart continue-with current-continuation error error-continuation error-in-thread error-thread ifcc ignore-error ignore-error/f ignore-errors in-callback? original-error recover restart restart? restarts rethrow rethrow-restarts return return-continuation thread-error-hook throw-continue throw-restarts with-datastack with-return
+syn keyword factorKeyword <condition> <continuation> <restart> attempt-all attempt-all-error attempt-all-error? callback-error-hook callcc0 callcc1 cleanup compute-restarts condition condition? continuation continuation? continue continue-restart continue-with current-continuation error error-continuation error-in-thread error-thread finally ifcc ignore-error ignore-error/f ignore-errors in-callback? original-error recover restart restart? restarts rethrow rethrow-restarts return return-continuation thread-error-hook throw-continue throw-restarts with-datastack with-return
 
 
 syn cluster factorReal          contains=factorInt,factorFloat,factorPosRatio,factorNegRatio,factorBinary,factorHex,factorOctal

removed/tokyo/assoc-functor/assoc-functor.factor

@@ -37,10 +37,10 @@ M: TYPE assoc-size handle>> DBRNUM ;
 : DBKEYS ( db -- keys )
     [ assoc-size <vector> ] [ handle>> ] bi
     dup DBITERINIT drop 0 int <ref>
-    [ 2dup DBITERNEXT dup ] [
+    [ 2dup DBITERNEXT ] [
         [ memory>object ] [ tcfree ] bi
         reach push
-    ] while 3drop ;
+    ] while* 2drop ;
 
 M: TYPE >alist
     [ DBKEYS dup ] keep '[ dup _ at 2array ] map! drop ;

vm/Config.windows

@@ -1,4 +1,4 @@
-SITE_CFLAGS += -mno-cygwin -mwindows
+SITE_CFLAGS += -mwindows
 CFLAGS_CONSOLE += -mconsole
 SHARED_FLAG = -shared
 SHARED_DLL_EXTENSION=.dll
@@ -7,7 +7,7 @@ LIBS = -lm
 
 PLAF_DLL_OBJS += vm/os-windows.o vm/mvm-windows.o
 PLAF_EXE_OBJS += vm/resources.o vm/main-windows.o
-PLAF_MASTER_HEADERS += vm/os-windows.hpp vm/mvm-windows.hpp
+PLAF_MASTER_HEADERS += vm/os-windows.hpp
 
 EXE_SUFFIX=
 EXE_EXTENSION=.exe
@@ -15,5 +15,5 @@ DLL_SUFFIX=
 DLL_EXTENSION=.dll
 CONSOLE_EXTENSION=.com
 
-LINKER = $(CPP) -shared -mno-cygwin -o 
+LINKER = $(CPP) -shared -o 
 LINK_WITH_ENGINE = -l$(DLL_PREFIX)factor$(DLL_SUFFIX)

vm/Config.windows.x86.32

@@ -1,6 +1,6 @@
 PLAF_DLL_OBJS += vm/os-windows-x86.32.o
 PLAF_MASTER_HEADERS += vm/os-windows.32.hpp
 DLL_PATH=http://factorcode.org/dlls
-WINDRES=windres
+WINDRES=windres -F pe-i386
 include vm/Config.windows
 include vm/Config.x86.32

vm/master.hpp

@@ -43,6 +43,8 @@
 #elif defined(__INTEL_COMPILER)
 #define FACTOR_COMPILER_VERSION \
   "Intel C Compiler " FACTOR_STRINGIZE(__INTEL_COMPILER)
+#elif defined(__MINGW32__)
+#define FACTOR_COMPILER_VERSION "MinGW (GCC " __VERSION__ ")"
 #elif defined(__GNUC__)
 #define FACTOR_COMPILER_VERSION "GCC " __VERSION__
 #elif defined(_MSC_FULL_VER)
@@ -79,7 +81,7 @@
 #error "Unsupported architecture"
 #endif
 
-#if defined(_MSC_VER)
+#if defined(_MSC_VER) || defined (__MINGW32__)
 #define WINDOWS
 #define WINNT
 #elif defined(WIN32)

vm/os-windows.cpp

@@ -185,12 +185,14 @@ uint64_t nano_count() {
 
 void sleep_nanos(uint64_t nsec) { Sleep((DWORD)(nsec / 1000000)); }
 
+#ifndef EXCEPTION_DISPOSITION
 typedef enum _EXCEPTION_DISPOSITION {
   ExceptionContinueExecution = 0,
   ExceptionContinueSearch = 1,
   ExceptionNestedException = 2,
   ExceptionCollidedUnwind = 3
 } EXCEPTION_DISPOSITION;
+#endif
 
 LONG factor_vm::exception_handler(PEXCEPTION_RECORD e, void* frame, PCONTEXT c,
                                   void* dispatch) {