--- /dev/null
+IN: homology
+USING: kernel sequences arrays math words namespaces
+hashtables prettyprint io ;
+
+! Utilities
+: S{ [ [ dup ] map>hash ] [ ] ; parsing
+
+: (lengthen) ( seq n -- seq )
+ over length - f <array> append ;
+
+: lengthen ( sim sim -- sim sim )
+ 2dup max-length tuck (lengthen) >r (lengthen) r> ;
+
+: unswons* 1 over tail swap first ;
+
+: swons* 1array swap append ;
+
+: rot-seq ( seq -- seq ) unswons* add ;
+
+: <point> ( -- sim ) gensym 1array ;
+
+: (C) ( point sim -- sim )
+ [ [ append natural-sort ] map-with ] map-with ;
+
+: (\/) ( sim sim -- sim ) lengthen [ append natural-sort ] 2map ;
+
+: <range> ( from to -- seq ) dup <slice> ;
+
+! Simplicial complexes
+SYMBOL: basepoint
+
+: {*} ( -- sim )
+ #! Initial object in category
+ { { { basepoint } } } ;
+
+: \/ ( sim sim -- sim )
+ #! Glue two complexes at base point
+ (\/) [ prune ] map ;
+
+: +point ( sim -- sim )
+ #! Adjoint an isolated point
+ unswons* <point> add swons* ;
+
+: C ( sim -- sim )
+ #! Cone on a space
+ <point> over first over add >r swap (C) r> swons* ;
+
+: S ( sim -- sim )
+ #! Suspension
+ [
+ <point> <point> 2dup 2array >r
+ pick (C) >r swap (C) r> (\/) r> swons*
+ ] keep (\/) ;
+
+: S^0 ( -- sim )
+ #! Degenerate sphere -- two points
+ {*} +point ;
+
+: S^ ( n -- sim )
+ #! Sphere
+ S^0 swap [ S ] times ;
+
+: D^ ( n -- sim )
+ #! Disc
+ 1- S^ C ;
+
+! Mod 2 matrix algebra
+: remove-1 ( n seq -- seq )
+ >r { } swap dup 1+ r> replace-slice ;
+
+: symmetric-diff ( hash hash -- hash )
+ clone swap [
+ drop dup pick hash [
+ over remove-hash
+ ] [
+ dup pick set-hash
+ ] if
+ ] hash-each ;
+
+SYMBOL: row-basis
+SYMBOL: matrix
+SYMBOL: current-row
+
+: rows ( -- n ) matrix get length ;
+
+: exchange-rows ( m n -- )
+ 2dup = [ 2drop ] [ matrix get exchange ] if ;
+
+: row ( n -- row ) matrix get nth ;
+
+: set-row ( row n -- ) matrix get set-nth ;
+
+: add-row ( src# dst# -- )
+ [ [ row ] 2apply symmetric-diff ] keep set-row ;
+
+: pivot-row ( basis-elt -- n )
+ current-row get rows <range>
+ [ row hash-member? ] find-with nip ;
+
+: kill-column ( basis-elt pivot -- )
+ dup 1+ rows <range> [
+ pick over row hash-member? [ dupd add-row ] [ drop ] if
+ ] each 2drop ;
+
+: with-matrix ( matrix basis quot -- matrix )
+ [
+ >r row-basis set matrix set r> call matrix get
+ ] with-scope ; inline
+
+: (row-reduce)
+ 0 current-row set
+ row-basis get [
+ dup pivot-row dup [
+ current-row get exchange-rows
+ current-row get kill-column
+ current-row inc
+ ] [
+ 2drop
+ ] if
+ ] each ;
+
+: ker/im ( -- ker im )
+ matrix get [ hash-empty? ] subset length
+ row-basis get [
+ matrix get [ hash-member? ] contains-with?
+ ] subset length ;
+
+: row-reduce ( matrix basis -- rowsp colsp matrix )
+ [ (row-reduce) ker/im ] with-matrix ;
+
+! Mod 2 homology
+: (boundary) ( seq -- chain )
+ dup length 1 <= [
+ H{ }
+ ] [
+ dup length [ over remove-1 dup ] map>hash
+ ] if nip ;
+
+: boundary ( chain -- chain )
+ H{ } swap [ drop (boundary) symmetric-diff ] hash-each ;
+
+: homology ( sim -- seq )
+ dup [ [ (boundary) ] map ] map rot-seq
+ [ row-reduce drop 2array ] 2map ;
+
+: print-matrix ( matrix basis -- )
+ swap [
+ swap [
+ ( row basis-elt )
+ swap hash-member? 1 0 ? pprint bl
+ ] each-with terpri
+ ] each-with ;
+
+2 S^ [ [ [ (boundary) ] map ] map unswons* drop ] keep
+[ [ row-reduce 2nip ] 2map ] keep
+[ print-matrix terpri ] 2each