diff --git a/collects/math/array.rkt b/collects/math/array.rkt
index 5d90c239a1..e1159c4c04 100644
--- a/collects/math/array.rkt
+++ b/collects/math/array.rkt
@@ -8,6 +8,7 @@
          "private/array/array-transform.rkt"
          "private/array/array-convert.rkt"
          "private/array/array-fold.rkt"
+         "private/array/array-special-folds.rkt"
          "private/array/array-print.rkt"
          "private/array/array-fft.rkt"
          "private/array/array-syntax.rkt"
@@ -34,6 +35,7 @@
           "private/array/array-transform.rkt"
           "private/array/array-convert.rkt"
           "private/array/array-fold.rkt"
+          "private/array/array-special-folds.rkt"
           "private/array/array-print.rkt"
           "private/array/array-syntax.rkt"
           "private/array/array-fft.rkt"
diff --git a/collects/math/private/array/array-fold.rkt b/collects/math/private/array/array-fold.rkt
index 532b12eb3f..4b5b0caa42 100644
--- a/collects/math/private/array/array-fold.rkt
+++ b/collects/math/private/array/array-fold.rkt
@@ -1,35 +1,32 @@
 #lang racket/base
 
-(require typed/untyped-utils
-         (except-in "typed-array-fold.rkt"
-                    array-axis-sum
-                    array-axis-prod
-                    array-axis-min
-                    array-axis-max
-                    array-all-sum
-                    array-all-prod
-                    array-all-min
-                    array-all-max))
+(require (for-syntax racket/base)
+         "typed-array-fold.rkt")
 
-(require/untyped-contract
- (begin (require "array-struct.rkt"))
- "typed-array-fold.rkt"
- [array-axis-sum  (case-> ((Array Number) Integer -> (Array Number))
-                          ((Array Number) Integer Number -> (Array Number)))]
- [array-axis-prod  (case-> ((Array Number) Integer -> (Array Number))
-                           ((Array Number) Integer Number -> (Array Number)))]
- [array-axis-min  (case-> ((Array Real) Integer -> (Array Real))
-                          ((Array Real) Integer Real -> (Array Real)))]
- [array-axis-max  (case-> ((Array Real) Integer -> (Array Real))
-                          ((Array Real) Integer Real -> (Array Real)))]
- [array-all-sum  (case-> ((Array Number) -> Number)
-                         ((Array Number) Number -> Number))]
- [array-all-prod  (case-> ((Array Number) -> Number)
-                          ((Array Number) Number -> Number))]
- [array-all-min  (case-> ((Array Real) -> Real)
-                         ((Array Real) Real -> Real))]
- [array-all-max  (case-> ((Array Real) -> Real)
-                         ((Array Real) Real -> Real))])
+;; ===================================================================================================
+;; Standard folds
+
+(define-syntax-rule (define-axis-fold name f)
+  (define-syntax (name stx)
+    (syntax-case stx ()
+      [(_ arr k)  (syntax/loc stx (array-axis-fold arr k f))]
+      [(_ arr k init)  (syntax/loc stx (array-axis-fold arr k f init))])))
+
+(define-syntax-rule (define-all-fold name f)
+  (define-syntax (name stx)
+    (syntax-case stx ()
+      [(_ arr)  (syntax/loc stx (array-all-fold arr f))]
+      [(_ arr init)  (syntax/loc stx (array-all-fold arr f init))])))
+
+(define-axis-fold array-axis-sum +)
+(define-axis-fold array-axis-prod *)
+(define-axis-fold array-axis-min min)
+(define-axis-fold array-axis-max max)
+
+(define-all-fold array-all-sum +)
+(define-all-fold array-all-prod *)
+(define-all-fold array-all-min min)
+(define-all-fold array-all-max max)
 
 (provide array-axis-fold
          array-axis-sum
@@ -37,17 +34,15 @@
          array-axis-min
          array-axis-max
          array-axis-count
-         array-axis-andmap
-         array-axis-ormap
+         array-axis-and
+         array-axis-or
          array-fold
+         array-all-fold
          array-all-sum
          array-all-prod
          array-all-min
          array-all-max
-         array-all-count
-         array-all-andmap
-         array-all-ormap
-         array-all-equal?
-         array-all-eqv?
-         array-all-eq?
-         array-all=)
+         array-all-and
+         array-all-or
+         array-axis-reduce
+         unsafe-array-axis-reduce)
diff --git a/collects/math/private/array/array-print.rkt b/collects/math/private/array/array-print.rkt
index 6bfd51de63..727abf4168 100644
--- a/collects/math/private/array/array-print.rkt
+++ b/collects/math/private/array/array-print.rkt
@@ -18,7 +18,10 @@
 (: print-array (All (A) ((Array A) Symbol Output-Port (U #t #f 0 1) -> Any)))
 ;; The logic in `print-array' causes the REPL printer to try printing an array in each layout, and
 ;; keep the first successful one. An overflowing line means failure.
-(define (print-array arr name port mode)
+(define (print-array orig-arr name port mode)
+  ;; Try to compute each element only once
+  (define arr (array-lazy orig-arr))
+  
   ;; Called to print array elements; may recur (e.g. printing arrays of arrays)
   ;; We never have to consider the `mode' argument again after defining `recur-print'
   (define recur-print
diff --git a/collects/math/private/array/array-special-folds.rkt b/collects/math/private/array/array-special-folds.rkt
new file mode 100644
index 0000000000..215d28d850
--- /dev/null
+++ b/collects/math/private/array/array-special-folds.rkt
@@ -0,0 +1,26 @@
+#lang typed/racket/base
+
+(require "array-struct.rkt"
+         "array-fold.rkt"
+         "array-pointwise.rkt")
+
+(provide array-count)
+
+(: array-count
+   (All (A B T ...)
+        (case-> ((A -> Any) (Array A) -> Index)
+                ((A B T ... T -> Any) (Array A) (Array B) (Array T) ... T -> Index))))
+(define array-count
+  (case-lambda:
+    [([f : (A -> Any)] [arr0 : (Array A)])
+     (assert (array-all-sum (inline-array-map (λ: ([a : A]) (if (f a) 1 0)) arr0)) index?)]
+    [([f : (A B -> Any)] [arr0 : (Array A)] [arr1 : (Array B)])
+     (assert
+      (array-all-sum (inline-array-map (λ: ([a : A] [b : B]) (if (f a b) 1 0)) arr0 arr1))
+      index?)]
+    [([f : (A B T ... T -> Any)] [arr0 : (Array A)] [arr1 : (Array B)] . [arrs : (Array T) ... T])
+     (assert
+      (array-all-sum (apply array-map
+                            (λ: ([a : A] [b : B] . [ts : T ... T]) (if (apply f a b ts) 1 0))
+                            arr0 arr1 arrs))
+      index?)]))
diff --git a/collects/math/private/array/typed-array-fold.rkt b/collects/math/private/array/typed-array-fold.rkt
index 4c65e8a93c..ed299b12cc 100644
--- a/collects/math/private/array/typed-array-fold.rkt
+++ b/collects/math/private/array/typed-array-fold.rkt
@@ -1,11 +1,11 @@
 #lang typed/racket/base
 
 (require racket/performance-hint
+         racket/fixnum
          "../unsafe.rkt"
          "array-struct.rkt"
          "array-indexing.rkt"
-         "utils.rkt"
-         "for-each.rkt")
+         "utils.rkt")
 
 (provide (all-defined-out))
 
@@ -16,16 +16,32 @@
 (define (check-array-axis name arr k)
   (define dims (array-dims arr))
   (cond
-    [(= dims 0)  (raise-argument-error 'name "Array with at least one axis" 0 arr k)]
+    [(fx= dims 0)  (raise-argument-error name "Array with at least one axis" 0 arr k)]
     [(or (0 . > . k) (k . >= . dims))
-     (raise-argument-error 'name (format "Index < ~a" dims) 1 arr k)]
+     (raise-argument-error name (format "Index < ~a" dims) 1 arr k)]
     [else  k]))
 
+(: array-axis-reduce (All (A B) ((Array A) Integer (Index (Integer -> A) -> B) -> (Array B))))
+(define (array-axis-reduce arr k f)
+  (let ([k  (check-array-axis 'array-axis-reduce arr k)])
+    (define ds (array-shape arr))
+    (define dk (unsafe-vector-ref ds k))
+    (define new-ds (unsafe-vector-remove ds k))
+    (define proc (unsafe-array-proc arr))
+    (unsafe-build-array
+     new-ds (λ: ([js : Indexes])
+              (define old-js (unsafe-vector-insert js k 0))
+              (f dk (λ: ([jk : Integer])
+                      (cond [(or (jk . < . 0) (jk . >= . dk))
+                             (raise-argument-error 'array-axis-reduce (format "Index < ~a" dk) jk)]
+                            [else
+                             (unsafe-vector-set! old-js k jk)
+                             (proc (vector-copy-all old-js))])))))))
+
 (: unsafe-array-axis-reduce (All (A B) ((Array A) Index (Index (Index -> A) -> B) -> (Array B))))
 (begin-encourage-inline
   (define (unsafe-array-axis-reduce arr k f)
     (define ds (array-shape arr))
-    (define dims (vector-length ds))
     (define dk (unsafe-vector-ref ds k))
     (define new-ds (unsafe-vector-remove ds k))
     (define proc (unsafe-array-proc arr))
@@ -42,18 +58,18 @@
     (unsafe-array-axis-reduce
      arr k (λ: ([dk : Index] [proc : (Index -> A)])
              (let: loop : B ([jk : Nonnegative-Fixnum  0] [acc : B  init])
-               (cond [(jk . < . dk)  (loop (+ jk 1) (f (proc jk) acc))]
+               (cond [(jk . fx< . dk)  (loop (fx+ jk 1) (f (proc jk) acc))]
                      [else  acc]))))))
 
 (: array-axis-fold/no-init (All (A) ((Array A) Integer (A A -> A) -> (Array A))))
 (define (array-axis-fold/no-init arr k f)
   (let ([k  (check-array-axis 'array-axis-fold arr k)])
-    (when (= (unsafe-vector-ref (array-shape arr) k) 0)
+    (when (fx= (unsafe-vector-ref (array-shape arr) k) 0)
       (raise-argument-error 'array-axis-fold "nonzero axis" 0 arr k))
     (unsafe-array-axis-reduce
      arr k (λ: ([dk : Index] [proc : (Index -> A)])
              (let: loop : A ([jk : Nonnegative-Fixnum  1] [acc : A  (proc 0)])
-               (cond [(jk . < . dk)  (loop (+ jk 1) (f (proc jk) acc))]
+               (cond [(jk . fx< . dk)  (loop (fx+ jk 1) (f (proc jk) acc))]
                      [else  acc]))))))
 
 (: array-axis-fold (All (A B) (case-> ((Array A) Integer (A A -> A) -> (Array A))
@@ -63,57 +79,33 @@
     [(arr k f)  (array-axis-fold/no-init arr k f)]
     [(arr k f init)  (array-axis-fold/init arr k f init)]))
 
-;; ---------------------------------------------------------------------------------------------------
-;; Whole-array fold
+;; ===================================================================================================
+;; Whole-array folds
 
-(: array-fold (All (A) ((Array A) ((Array A) Index -> (Array A)) -> (Array A))))
 (begin-encourage-inline
+  
+  (: array-fold (All (A) ((Array A) ((Array A) Index -> (Array A)) -> (Array A))))
   (define (array-fold arr f)
     (define dims (array-dims arr))
     (let loop ([#{k : Index} dims] [arr arr])
-      (cond [(zero? k)  arr]
-            [else  (let ([k  (sub1 k)])
-                     (loop k (f arr k)))]))))
-
-;; ===================================================================================================
-;; Standard axis folds
-
-(define-syntax-rule (define-axis-fold name op T ...)
-  (begin-encourage-inline
-    (: name (case-> ((Array T) Integer -> (Array T)) ...
-                    ((Array T) Integer T -> (Array T)) ...))
-    (define name
-      (case-lambda
-        [(arr k)  (array-axis-fold arr k op)]
-        [(arr k init)  (array-axis-fold arr k op init)]))))
-
-(define-axis-fold array-axis-sum + Float Real Float-Complex Number)
-(define-axis-fold array-axis-prod * Float Real Float-Complex Number)
-(define-axis-fold array-axis-min min Float Real)
-(define-axis-fold array-axis-max max Float Real)
-
-;; ===================================================================================================
-;; Standard whole-array folds
-
-(define-syntax-rule (define-fold name array-axis-op T ...)
-  (begin-encourage-inline
-    (: name (case-> ((Array T) -> T) ...
-                    ((Array T) T -> T) ...))
-    (define name 
-      (case-lambda
-        [(arr)  (array-ref (array-fold arr array-axis-op) #())]
-        [(arr init)
-         (plet: (A) ([arr : (Array A)  arr]
-                     [array-axis-op : ((Array A) Index A -> (Array A))  array-axis-op]
-                     [init : A  init])
-                (array-ref (array-fold arr (λ: ([arr : (Array A)] [k : Index])
-                                             (array-axis-op arr k init)))
-                           #()))]))))
-
-(define-fold array-all-sum array-axis-sum Float Real Float-Complex Number)
-(define-fold array-all-prod array-axis-prod Float Real Float-Complex Number)
-(define-fold array-all-min array-axis-min Float Real)
-(define-fold array-all-max array-axis-max Float Real)
+      (cond [(fx= k 0)  arr]
+            [else  (let ([k  (fx- k 1)])
+                     (loop k (f arr k)))])))
+  
+  (: array-all-fold (All (A) (case-> ((Array A) (A A -> A) -> A)
+                                     ((Array A) (A A -> A) A -> A))))
+  (define array-all-fold
+    (case-lambda
+      [(arr f)
+       (array-ref (array-fold arr (λ: ([arr : (Array A)] [k : Index])
+                                    (array-axis-fold arr k f)))
+                  #())]
+      [(arr f init)
+       (array-ref (array-fold arr (λ: ([arr : (Array A)] [k : Index])
+                                    (array-axis-fold arr k f init)))
+                  #())]))
+  
+  )  ; begin-encourage-inline
 
 ;; ===================================================================================================
 ;; Count
@@ -124,69 +116,38 @@
     (unsafe-array-axis-reduce
      arr k (λ: ([dk : Index] [proc : (Index -> A)])
              (let: loop : Index ([jk : Nonnegative-Fixnum  0] [acc : Nonnegative-Fixnum  0])
-               (if (jk . < . dk)
-                   (cond [(pred? (proc jk))  (loop (+ jk 1) (unsafe-fx+ acc 1))]
-                         [else  (loop (+ jk 1) acc)])
+               (if (jk . fx< . dk)
+                   (cond [(pred? (proc jk))  (loop (fx+ jk 1) (unsafe-fx+ acc 1))]
+                         [else  (loop (fx+ jk 1) acc)])
                    (assert acc index?)))))))
 
-(: array-all-count (All (A) ((Array A) (A -> Any) -> Index)))
-(define (array-all-count arr pred?)
-  (define: i : (Boxof Nonnegative-Fixnum) (box 0))
-  (define proc (unsafe-array-proc arr))
-  (define ds (array-shape arr))
-  (for-each-array-index ds (λ (js) (when (pred? (proc js)) (set-box! i (unsafe-fx+ (unbox i) 1)))))
-  (assert (unbox i) index?))
-
 ;; ===================================================================================================
-;; Short-cutting andmap
+;; Short-cutting axis folds
 
-(: array-axis-andmap (All (A) ((Array A) Integer (A -> Any) -> (Array Boolean))))
-(define (array-axis-andmap arr k pred?)
-  (let ([k  (check-array-axis 'array-axis-andmap arr k)])
+(: array-axis-and (All (A) ((Array A) Integer -> (Array (U A Boolean)))))
+(define (array-axis-and arr k)
+  (let ([k  (check-array-axis 'array-axis-and arr k)])
     (unsafe-array-axis-reduce
      arr k (λ: ([dk : Index] [proc : (Index -> A)])
-             (let: loop : Boolean ([jk : Nonnegative-Fixnum  0])
-               (cond [(jk . < . dk)  (if (pred? (proc jk)) (loop (+ jk 1)) #f)]
-                     [else  #t]))))))
+             (let: loop : (U A Boolean) ([jk : Nonnegative-Fixnum  0] [acc : (U A Boolean)  #t])
+               (cond [(jk . fx< . dk)  (define v (and acc (proc jk)))
+                                       (if v (loop (fx+ jk 1) v) v)]
+                     [else  acc]))))))
 
-(: array-all-andmap (All (A) ((Array A) (A -> Any) -> Boolean)))
-(define (array-all-andmap arr pred?)
-  (let/ec: return : Boolean
-    (define proc (unsafe-array-proc arr))
-    (define ds (array-shape arr))
-    (for-each-array-index ds (λ (js) (unless (pred? (proc js)) (return #f))))
-    #t))
-
-;; ===================================================================================================
-;; Short-cutting ormap
-
-(: array-axis-ormap (All (A) ((Array A) Integer (A -> Any) -> (Array Boolean))))
-(define (array-axis-ormap arr k pred?)
-  (let ([k  (check-array-axis 'array-axis-ormap arr k)])
+(: array-axis-or (All (A) ((Array A) Integer -> (Array (U A #f)))))
+(define (array-axis-or arr k)
+  (let ([k  (check-array-axis 'array-axis-or arr k)])
     (unsafe-array-axis-reduce
      arr k (λ: ([dk : Index] [proc : (Index -> A)])
-             (let: loop : Boolean ([jk : Nonnegative-Fixnum  0])
-               (cond [(jk . < . dk)  (if (pred? (proc jk)) #t (loop (+ jk 1)))]
-                     [else  #f]))))))
+             (let: loop : (U A #f) ([jk : Nonnegative-Fixnum  0] [acc : (U A #f)  #f])
+               (cond [(jk . fx< . dk)  (define v (or acc (proc jk)))
+                                       (if v v (loop (fx+ jk 1) v))]
+                     [else  acc]))))))
 
-(: array-all-ormap (All (A) ((Array A) (A -> Any) -> Boolean)))
-(define (array-all-ormap arr pred?)
-  (let/ec: return : Boolean
-    (define f (unsafe-array-proc arr))
-    (define ds (array-shape arr))
-    (for-each-array-index ds (λ (js) (when (pred? (f js)) (return #t))))
-    #f))
+(: array-all-and (All (A B) ((Array A) -> (U A Boolean))))
+(define (array-all-and arr)
+  (array-ref ((inst array-fold (U A Boolean)) arr array-axis-and) #()))
 
-;; ===================================================================================================
-
-(: array-all-equal? ((Array Any) (Array Any) -> Boolean))
-(define array-all-equal? equal?)
-
-(: array-all-eqv? ((Array Any) (Array Any) -> Boolean))
-(define array-all-eqv? (array-lift-comparison eqv?))
-
-(: array-all-eq? ((Array Any) (Array Any) -> Boolean))
-(define array-all-eq? (array-lift-comparison eq?))
-
-(: array-all= ((Array Number) (Array Number) -> Boolean))
-(define array-all= (array-lift-comparison =))
+(: array-all-or (All (A B) ((Array A) -> (U A #f))))
+(define (array-all-or arr)
+  (array-ref ((inst array-fold (U A #f)) arr array-axis-or) #()))
diff --git a/collects/math/private/matrix/matrix-types.rkt b/collects/math/private/matrix/matrix-types.rkt
index 6541fdad31..4bb1332680 100644
--- a/collects/math/private/matrix/matrix-types.rkt
+++ b/collects/math/private/matrix/matrix-types.rkt
@@ -37,7 +37,8 @@
   (vector-ref (array-shape a) 0))
 
 (: matrix-all= : (Matrix Number) (Matrix Number) -> Boolean)
-(define matrix-all= array-all=)
+(define (matrix-all= arr0 arr1)
+  (array-all-and (array= arr0 arr1)))
 
 (: matrix-dimensions : (Matrix Number) -> (Values Index Index))
 (define (matrix-dimensions a)
diff --git a/collects/math/scribblings/math-array.scrbl b/collects/math/scribblings/math-array.scrbl
index 4cb6819558..194c86cbba 100644
--- a/collects/math/scribblings/math-array.scrbl
+++ b/collects/math/scribblings/math-array.scrbl
@@ -5,9 +5,9 @@
           (for-label racket/base racket/vector racket/match racket/unsafe/ops racket/string
                      math plot
                      (only-in typed/racket/base
-                              ann inst : λ: define: make-predicate
+                              ann inst : λ: define: make-predicate ->
                               Flonum Real Boolean Any Integer Index Natural Exact-Positive-Integer
-                              Nonnegative-Real Sequenceof Fixnum Values
+                              Nonnegative-Real Sequenceof Fixnum Values Number
                               All U List Vector Listof Vectorof Struct))
           "utils.rkt")
 
@@ -16,7 +16,8 @@
                          (require racket/match
                                   racket/vector
                                   racket/string
-                                  racket/sequence)]
+                                  racket/sequence
+                                  racket/list)]
 
 @title[#:tag "arrays" #:style 'toc]{Arrays}
 @(author-neil)
@@ -52,6 +53,10 @@ applied to indexes to retrieve elements.
 
 @local-table-of-contents[]
 
+
+@;{==================================================================================================}
+
+
 @section{Preliminaries}
 
 @subsection{Definitions}
@@ -296,6 +301,10 @@ of stretching just singleton axes:
 Notice that @racket[(array #["+" "-"])] was repeated five times, and that
 @racket[arr3] was repeated three full times and once partially.
 
+
+@;{==================================================================================================}
+
+
 @section{Types, Predicates and Accessors}
 
 @defform[(Array A)]{
@@ -433,6 +442,10 @@ Returns the number of @racket[arr]'s dimensions. Equivalent to
 Returns the vector of data that @racket[arr] contains.
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Construction}
 
 @defform[(array #[#[...] ...])]{
@@ -614,6 +627,10 @@ have the value @racket[on-value]; the rest have the value @racket[off-value].
                  (diagonal-array 2 7 1 0)]
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Conversion}
 
 @deftogether[(@defproc[(list->array [lst (Listof A)]) (Array A)]
@@ -732,6 +749,10 @@ They call the current @racket[array-custom-printer]:
                      #t)]
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Comprehensions and Sequences}
 
 Sometimes sequential processing is unavoidable, so @racket[math/array] provides loops and sequences.
@@ -805,6 +826,10 @@ Returns a sequence of indexes for shape @racket[ds], in row-major order.
                  (indexes-array #(3 3))]
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Pointwise Operations}
 
 Most of the operations documented in this section are simple macros that apply @racket[array-map]
@@ -992,6 +1017,10 @@ This is used for both @racket[(array-broadcasting #t)] and @racket[(array-broadc
                  (array-broadcast (array #[0 1]) ((inst vector Index) 5))]
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Indexing and Slicing}
 
 @defproc[(array-ref [arr (Array A)] [js In-Indexes]) A]{
@@ -1236,6 +1265,10 @@ Given a slice @racket[s] and an axis length @racket[dk], returns the arguments t
 that would produce an equivalent slice specification.
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Transformations}
 
 @defproc[(array-transform [arr (Array A)] [ds In-Indexes] [proc (Indexes -> In-Indexes)])
@@ -1257,8 +1290,8 @@ Double an array in every dimension by duplicating elements:
                      (vector-map (λ: ([d : Index]) (* d 2)) (array-shape arr))
                      (λ: ([js : Indexes])
                        (vector-map (λ: ([j : Index]) (quotient j 2)) js)))]
-Recall that, because @racket[array-transform] returns @tech{non-strict} arrays, the above
-result takes little more space than the original array.
+Because @racket[array-transform] returns @tech{non-strict} arrays, the above result takes
+little more space than the original array.
 
 Almost all array transformations, including those effected by @secref{array:slice-specs}, are
 implemented using @racket[array-transform] or its unsafe counterpart.
@@ -1314,7 +1347,7 @@ Returns an array like @racket[arr], but with axes permuted according to @racket[
 
 The list @racket[perm] represents a mapping from source axis numbers to destination
 axis numbers: the source is the list position, the destination is the list element.
-For example, the permutation @racket['(0 1)] is the identity permutation for two-dimensional
+For example, the permutation @racket['(0 1 2)] is the identity permutation for three-dimensional
 arrays, @racket['(1 0)] swaps axes @racket[0] and @racket[1], and @racket['(3 1 2 0)] swaps
 axes @racket[0] and @racket[3].
 
@@ -1344,36 +1377,209 @@ Returns an array with shape @racket[(vector (array-size arr))], with the element
                  (array-flatten (array #[#[0 1] #[2 3]]))]
 }
 
-@section{Folds}
 
-@;{
-array-axis-fold
-array-fold
-array-all-sum
-array-all-prod
-array-all-min
-array-all-max
-array-axis-count
-array-all-count
-array-axis-andmap
-array-all-andmap
-array-axis-ormap
-array-all-ormap
-array-all-eq?
-array-all-equal?
-array-all-eqv?
-array-all=
-array-axis-sum
-array-axis-prod
-array-axis-min
-array-axis-max
-array-lift-comparison
-array-axis-fft
-array-fft
+@;{==================================================================================================}
+
+
+@section{Folds and Other Axis Reductions}
+
+@defproc*[([(array-axis-fold [arr (Array A)] [k Integer] [f (A A -> A)]) (Array A)]
+           [(array-axis-fold [arr (Array A)] [k Integer] [f (A B -> B)] [init B]) (Array B)])]{
+Folds a binary function @racket[f] over axis @racket[k] of @racket[arr]. The result has the
+shape of @racket[arr] but with axis @racket[k] removed.
+
+The three-argument variant uses the first value of each row in axis @racket[k] as @racket[init].
+It therefore requires axis @racket[k] to have positive length.
+
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 4)))
+                 arr
+                 (array-axis-fold arr 0 +)
+                 (array-axis-fold arr 1 (inst cons Index (Listof Index)) empty)]
+Notice that the second example returns an array of reversed lists.
+This is therefore a left fold; see @racket[foldl].
+}
+
+@deftogether[(@defform*[((array-axis-sum arr k)
+                         (array-axis-sum arr k init))]
+              @defform*[((array-axis-prod arr k)
+                         (array-axis-prod arr k init))
+                        #:contracts ([arr   (Array Number)]
+                                     [k     Integer]
+                                     [init  Number])])]{
+}
+@deftogether[(@defform*[((array-axis-min arr k)
+                         (array-axis-min arr k init))]
+              @defform*[((array-axis-max arr k)
+                         (array-axis-max arr k init))
+                        #:contracts ([arr   (Array Real)]
+                                     [k     Integer]
+                                     [init  Real])])]{
+Some standard per-axis folds, defined in terms of @racket[array-axis-fold]. The two-argument
+variants require axis @racket[k] to have positive length.
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 4)))
+                 arr
+                 (array-axis-fold arr 1 +)
+                 (array-axis-sum arr 1)
+                 (array-axis-sum arr 0 0.0)]
+}
+
+@defproc[(array-fold [arr (Array A)] [g ((Array A) Index -> (Array A))]) (Array A)]{
+Folds @racket[g] over @italic{each axis} of @racket[arr], in reverse order. The arguments
+to @racket[g] are an array (initially @racket[arr]) and the current axis.
+It should return an array with one fewer dimension than the array given, but does not have to.
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 4)))
+                 arr
+                 (array-fold arr (λ: ([arr : (Array Integer)] [k : Index])
+                                   (array-axis-sum arr k)))
+                 (+ 0 1 2 3 4 5 6 7 8 9 10 11)
+                 (array-fold
+                  arr
+                  (λ: ([arr : (Array (Listof* Index))] [k : Index])
+                    (array-map
+                     (inst reverse (Listof* Index))
+                     (array-axis-fold arr k
+                                      (inst cons (Listof* Index) (Listof (Listof* Index)))
+                                      empty))))]
+}
+
+@defproc*[([(array-all-fold [arr (Array A)] [f (A A -> A)]) A]
+           [(array-all-fold [arr (Array A)] [f (A A -> A)] [init A]) A])]{
+Folds @racket[f] over every element of @racket[arr] by folding @racket[f] over each axis in
+reverse order. The two-argument variant is equivalent to
+@racketblock[(array-ref (array-fold arr (λ: ([arr : (Array A)] [k : Index])
+                                          (array-axis-fold arr k f)))
+                        #())]
+and the three-argument variant is similar. The two-argument variant requires every axis to have
+positive length.
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 4)))
+                 arr
+                 (array-all-fold arr +)
+                 (array-all-fold (array #[]) + 0.0)]
+}
+
+@deftogether[(@defform*[((array-all-sum arr)
+                         (array-all-sum arr init))]
+              @defform*[((array-all-prod arr)
+                         (array-all-prod arr init))
+                        #:contracts ([arr   (Array Number)]
+                                     [init  Number])])]{
+}
+@deftogether[(@defform*[((array-all-min arr)
+                         (array-all-min arr init))]
+              @defform*[((array-all-max arr)
+                         (array-all-max arr init))
+                        #:contracts ([arr   (Array Real)]
+                                     [init  Real])])]{
+Some standard whole-array folds, defined in terms of @racket[array-all-fold]. The one-argument
+variants require each axis in @racket[arr] to have positive length.
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 4)))
+                 arr
+                 (array-all-fold arr +)
+                 (array-all-sum arr)
+                 (array-all-sum arr 0.0)]
+}
+
+@defproc[(array-axis-count [arr (Array A)] [k Integer] [pred? (A -> Any)]) (Array Index)]{
+Counts the elements @racket[x] in rows of axis @racket[k] for which @racket[(pred? x)] is true.
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 3)))
+                 arr
+                 (array-axis-count arr 1 odd?)]
+}
+
+@defproc*[([(array-count [pred? (A -> Any)] [arr0 (Array A)]) Index]
+           [(array-count [pred? (A B Ts ... -> Any)]
+                         [arr0 (Array A)]
+                         [arr1 (Array B)]
+                         [arrs (Array Ts)] ...)
+            Index])]{
+The two-argument variant returns the number of elements @racket[x] in @racket[arr] for which
+@racket[(pred? x)] is true. The other variant does the same with the corresponding elements from
+any number of arrays. If the arrays' shapes are not the same, they are @tech{broadcast} first.
+@examples[#:eval typed-eval
+                 (array-count zero? (array #[#[0 1 0 2] #[0 3 -1 4]]))
+                 (array-count equal?
+                              (array #[#[0 1] #[2 3] #[0 1] #[2 3]])
+                              (array #[0 1]))]
+}
+
+@deftogether[(@defproc[(array-axis-and [arr (Array A)] [k Integer]) (Array (U A Boolean))]
+              @defproc[(array-axis-or [arr (Array A)] [k Integer]) (Array (U A #f))])]{
+Apply @racket[and] or @racket[or] to each row in axis @racket[k] of array @racket[arr], using
+short-cut evaluation.
+
+Consider the following array, whose second element takes 10 seconds to compute:
+@interaction[#:eval typed-eval
+                    (define arr (build-array #(2) (λ: ([js : Indexes])
+                                                    (cond [(zero? (vector-ref js 0))  #f]
+                                                          [else  (sleep 10)
+                                                                 #t]))))]
+Printing it takes over 10 seconds, but this returns immediately:
+@interaction[#:eval typed-eval
+                    (array-axis-and arr 0)]
+}
+
+@deftogether[(@defproc[(array-all-and [arr (Array A)]) (U A Boolean)]
+              @defproc[(array-all-or [arr (Array A)]) (U A #f)])]{
+Apply @racket[and] or @racket[or] to each element in @racket[arr] using short-cut evaluation.
+
+@racket[(array-all-and arr)] is defined as
+@racketblock[(array-ref (array-fold arr array-axis-and) #())]
+and @racket[array-all-or] is defined similarly.
+
+@examples[#:eval typed-eval
+                 (define arr (index-array #(3 3)))
+                 (array-all-and (array= arr arr))
+                 (define brr (array+ arr (array 1)))
+                 (array-all-and (array= arr brr))
+                 (array-all-or (array= arr (array 0)))]
+}
+
+@defproc[(array-axis-reduce [arr (Array A)] [k Integer] [h (Index (Integer -> A) -> B)]) (Array B)]{
+Like @racket[array-axis-fold], but allows evaluation control (such as short-cutting @racket[and] and
+@racket[or]) by moving the loop into @racket[h]. The result has the shape of @racket[arr], but with
+axis @racket[k] removed.
+
+The arguments to @racket[h] are the length of axis @racket[k] and a procedure that retrieves
+elements from that axis's rows by their indexes in axis @racket[k]. It should return the elements
+of the resulting array.
+
+For example, summing the squares of the rows in axis @racket[1]:
+@interaction[#:eval typed-eval
+                    (define arr (index-array #(3 3)))
+                    arr
+                    (array-axis-reduce
+                     arr 1
+                     (λ: ([dk : Index] [proc : (Integer -> Real)])
+                       (for/fold: ([s : Real 0]) ([jk  (in-range dk)])
+                         (+ s (sqr (proc jk))))))
+                    (array-axis-sum (array-map sqr arr) 1)]
+
+Every fold, including @racket[array-axis-fold], is ultimately defined using
+@racket[array-axis-reduce] or its unsafe counterpart.
+}
+
+
+@;{==================================================================================================}
+
+
+@section{Other Array Operations}
+
+@;{array-fft
 array-inverse-fft
+array-axis-fft
 array-axis-inverse-fft
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Subtypes}
 
 @subsection{Flonum Arrays}
@@ -1382,7 +1588,7 @@ array-axis-inverse-fft
 }
 
 @defproc[(array->flarray [arr (Array Real)]) FlArray]{
-Returns an flarray that has approximately the same elements as @racket[arr].
+Returns an flonum array that has approximately the same elements as @racket[arr].
 The elements may lose precision during the conversion.
 }
 
@@ -1429,7 +1635,7 @@ flarray>=
 }
 
 @defproc[(array->fcarray [arr (Array Number)]) FCArray]{
-Returns an fcarray that has approximately the same elements as @racket[arr].
+Returns a float-complex array that has approximately the same elements as @racket[arr].
 The elements may lose precision during the conversion.
 }
 
@@ -1466,6 +1672,10 @@ fcarray-asin
 fcarray-atan
 }
 
+
+@;{==================================================================================================}
+
+
 @section{Unsafe Operations}
 
 @;{
@@ -1481,6 +1691,7 @@ unsafe-build-array
 unsafe-mutable-array
 unsafe-flarray
 unsafe-array-transform
+unsafe-array-axis-reduce
 }
 
 @;{