source: project/wiki/eggref/5/srfi-179 @ 39515

Last change on this file since 39515 was 39515, checked in by Diego, 5 months ago

Wrangle srfi-179 document into an overly verbose wiki page

File size: 123.2 KB
Line 
1<nowiki>
2<script>
3MathJax = {
4  tex: {
5    inlineMath: [['$', '$'], ['\\(', '\\)']]
6  }
7};
8</script>
9<script type="text/javascript" id="MathJax-script" async
10  src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js">
11</script>
12</nowiki>
13
14[[tags: egg]]
15
16[[toc:]]
17
18== srfi-179
19
20Nonempty Intervals and Generalized Arrays (Updated)
21
22=== Chicken implementation and documentation notes
23
24'''NOTE:''' These are just differences to the reference implementation, not
25deviations from the SRFI specification. If any behavior is different from the
26SRFI, please contact the egg maintainer.
27
28'''NOTE:''' This document is mostly the original SRFI html wrangled into wiki
29syntax using pandoc and manual edits. If anything looks wrong, please contact
30the egg maintainer.
31
32* Complex number vector operations ({{c64}}, {{c128}}) are provided by srfi-160
33* {{specialized-array-default-safe?}} and {{specialized-aray-default-mutable?}} are implemented in terms of chicken's {{make-parameter}}.
34* The original implementation's {{macro-absent-obj}} has been replaced by default optional parameter value {{#f}}.
35* Checkers for flonum storage classes ({{f32}}, {{f64}}, {{c64}}, {{c128}}) have been relaxed according to acceptable values to chicken's srfi-4 {{-set!}} procedures.
36
37=== Abstract
38
39This SRFI specifies an array mechanism for Scheme. Arrays as defined here are quite general; at their most basic, an array is simply a mapping, or function, from multi-indices of exact integers $i_0,\ldots,i_{d-1}$ to Scheme values. The set of multi-indices $i_0,\ldots,i_{d-1}$ that are valid for a given array form the ''domain'' of the array. In this SRFI, each array's domain consists of the cross product of nonempty intervals of exact integers $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ of $\mathbb Z^d$, $d$-tuples of integers. Thus, we introduce a data type called $d$-''intervals'', or more briefly ''intervals'', that encapsulates this notion. (We borrow this terminology from, e.g., Elias Zakon's [[http://www.trillia.com/zakon1.html|Basic Concepts of Mathematics]].) Specialized variants of arrays are specified to provide portable programs with efficient representations for common use cases.
40
41=== Rationale
42
43This SRFI was motivated by a number of somewhat independent notions, which we outline here and which are explained below.
44
45* Provide a '''general API''' (Application Program Interface) that specifies the minimal required properties of any given array, without requiring any specific implementation strategy from the programmer for that array.
46* Provide a '''single, efficient implementation for dense arrays''' (which we call ''specialized arrays'').
47* Provide '''useful array transformations''' by exploiting the algebraic structure of affine one-to-one mappings on multi-indices.
48* Separate '''the routines that specify the work to be done''' ({{array-map}}, {{array-outer-product}}, etc.) from '''the routines that actually do the work''' ({{array-copy}}, {{array-assign!}}, {{array-fold}}, etc.). This approach '''avoids temporary intermediate arrays''' in computations.
49* Encourage '''bulk processing of arrays''' rather than word-by-word operations.
50
51This SRFI differs from the finalized [[https://srfi.schemers.org/srfi-122/|SRFI 122]] in the following ways:
52
53* The procedures {{specialized-array-default-mutable?}}, {{interval-for-each}}, {{interval-cartesian-product}}, {{interval-rotate}} and {{array-elements-in-order?}}, {{array-outer-product}}, {{array-tile}}, {{array-rotate}}, {{array-reduce}}, {{array-assign!}}, {{array-ref}}, {{array-set!}}, and {{specialized-array-reshape}} have been added together with some examples.
54* Global variables {{f8-storage-class}} and {{f16-storage-class}} have been added.
55* Homogeneous storage classes must be implemented using homogeneous vectors, or be defined as {{#f}}.
56* The procedure {{make-interval}} now takes one or two arguments.
57* Specialized arrays can be mutable or immutable; the default, which can be changed, is mutable. Shared arrays inherit safety and mutability from source arrays.
58* The discussion of Haar transforms as examples of separable transforms has been corrected.
59* The documentation has a few more examples of image processing algorithms.
60* Some matrix examples have been added to this document.
61
62=== Overview
63
64==== Bawden-style arrays
65
66In a [[https://groups.google.com/forum/?hl=en#!msg/comp.lang.scheme/7nkx58Kv6RI/a5hdsduFL2wJ|1993 post]] to the news group comp.lang.scheme, Alan Bawden gave a simple implementation of multi-dimensional arrays in R4RS scheme. The only constructor of new arrays required specifying an initial value, and he provided the three low-level primitives {{array-ref}}, {{array-set!}}, and {{array?}}, as well as {{make-array}} and {{make-shared-array}}. His arrays were defined on rectangular intervals in $\mathbb Z^d$ of the form $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$. I'll note that his function {{array-set!}} put the value to be entered into the array at the front of the variable-length list of indices that indicate where to place the new value. He offered an intriguing way to "share" arrays in the form of a routine {{make-shared-array}} that took a mapping from a new interval of indices into the domain of the array to be shared. His implementation incorporated what he called an ''indexer'', which was a function from the interval $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$ to an interval $[0,N)$, where the ''body'' of the array consisted of a single Scheme vector of length $N$. Bawden called the mapping specified in {{make-shared-array}} ''linear'', but I prefer the term ''affine'', as I explain later.
67
68Mathematically, Bawden's arrays can be described as follows. We'll use the vector notation $\vec i$ for a multi-index $i_0,\ldots,i_{d-1}$. (Multi-indices correspond to Scheme {{values}}.) Arrays will be denoted by capital letters $A,B,\ldots$, the domain of the array $A$ will be denoted by $D_A$, and the indexer of $A$, mapping $D_A$ to the interval $[0,N)$, will be denoted by $I_A$. Initially, Bawden constructs $I_A$ such that $I_A(\vec i)$ steps consecutively through the values $0,1,\ldots,N-1$ as $\vec i$ steps through the multi-indices $(l_0,\ldots,l_{d-2},l_{d-1})$, $(l_0,\ldots,l_{d-2},l_{d-1}+1)$, $\ldots$, $(l_0,\ldots,l_{d-2}+1,l_{d-1})$, etc., in lexicographical order, which means that if $\vec i$ and $\vec j$ are two multi-indices, then $\vec i<\vec j$ iff the first coordinate $k$ where $\vec i$ and $\vec j$ differ satisfies $\vec i_k<\vec j_k$.
69
70In {{make-shared-array}}, Bawden allows you to specify a new $r$-dimensional interval $D_B$ as the domain of a new array $B$, and a mapping $T_{BA}:D_B\to D_A$ of the form $T_{BA}(\vec i)=M\vec i+\vec b$; here $M$ is a $d\times r$ matrix of integer values and $\vec b$ is a $d$-vector. So this mapping $T_{BA}$ is ''affine'', in that $T_{BA}(\vec i)-T_{BA}(\vec j)=M(\vec i-\vec j)$ is ''linear'' (in a linear algebra sense) in $\vec i-\vec j$. The new indexer of $B$ satisfies $I_B(\vec i)=I_A(T_{BA}(\vec i))$.
71
72A fact Bawden exploits in the code, but doesn't point out in the short post, is that $I_B$ is again an affine map, and indeed, the composition of ''any'' two affine maps is again affine.
73
74==== Our extensions of Bawden-style arrays
75
76We incorporate Bawden-style arrays into this SRFI, but extend them in one minor way that we find useful.
77
78We introduce the notion of a ''storage class'', an object that contains functions that manipulate, store, check, etc., different types of values. A {{generic-storage-class}} can manipulate any Scheme value, whereas, e.g., a {{u1-storage-class}} can store only the values 0 and 1 in each element of a body.
79
80We also require that our affine maps be one-to-one, so that if $\vec i\neq\vec j$ then $T(\vec i)\neq T(\vec j)$. Without this property, modifying the $\vec i$th component of $A$ would cause the $\vec j$th component to change.
81
82==== Common transformations on Bawden-style arrays
83
84Requiring the transformations $T_{BA}:D_B\to D_A$ to be affine may seem esoteric and restricting, but in fact many common and useful array transformations can be expressed in this way. We give several examples below:
85
86* '''Restricting the domain of an array:''' If the domain of $B$, $D_B$, is a subset of the domain of $A$, then $T_{BA}(\vec i)=\vec i$ is a one-to-one affine mapping. We define {{array-extract}} to define this common operation; it's like looking at a rectangular sub-part of a spreadsheet. We use it to extract the common part of overlapping domains of three arrays in an image processing example below.
87* '''Tiling an array:''' For various reasons (parallel processing, optimizing cache localization, GPU programming, etc.) one may wish to process a large array as a number of subarrays of the same dimensions, which we call ''tiling'' the array. The routine {{array-tile}} returns a new array, each entry of which is a subarray extracted (in the sense of {{array-extract}}) from the input array.
88* '''Translating the domain of an array:''' If $\vec d$ is a vector of integers, then $T_{BA}(\vec i)=\vec i-\vec d$ is a one-to-one affine map of $D_B=\{\vec i+\vec d\mid \vec i\in D_A\}$ onto $D_A$. We call $D_B$ the ''translate'' of $D_A$, and we define {{array-translate}} to provide this operation.
89* '''Permuting the coordinates of an array:''' If $\pi$ [[https://en.wikipedia.org/wiki/Permutation|permutes]] the coordinates of a multi-index $\vec i$, and $\pi^{-1}$ is the inverse of $\pi$, then $T_{BA}(\vec i)=\pi (\vec i)$ is a one-to-one affine map from $D_B=\{\pi^{-1}(\vec i)\mid \vec i\in D_A\}$ onto $D_A$. We provide {{array-permute}} for this operation. (The only nonidentity permutation of a two-dimensional spreadsheet turns rows into columns and vice versa.) We also provide {{array-rotate}} for the special permutations that rotate the axes. For example, in three dimensions we have the following three rotations: $i\ j\ k\to j\ k\ i$; $i\ j\ k\to k\ i\ j$; and the trivial (identity) rotation $i\ j\ k\to i\ j\ k$. The three-dimensional permutations that are not rotations are $i\ j\ k\to i\ k\ j$; $i\ j\ k\to j\ i\ k$; and $i\ j\ k\to k\ j\ i$.
90* '''Currying an array:''' Let's denote the cross product of two intervals $\text{Int}_1$ and $\text{Int}_2$ by $\text{Int}_1\times\text{Int}_2$; if $\vec j=(j_0,\ldots,j_{r-1})\in \text{Int}_1$ and $\vec i=(i_0,\ldots,i_{s-1})\in \text{Int}_2$, then $\vec j\times\vec i$, which we define to be $(j_0,\ldots,j_{r-1},i_0,\ldots,i_{s-1})$, is in $\text{Int}_1\times\text{Int}_2$. If $D_A=\text{Int}_1\times\text{Int}_2$ and $\vec j\in\text{Int}_1$, then $T_{BA}(\vec i)=\vec j\times\vec i$ is a one-to-one affine mapping from $D_B=\text{Int}_2$ into $D_A$. For each vector $\vec j$ we can compute a new array in this way; we provide {{array-curry}} for this operation, which returns an array whose domain is $\text{Int}_1$ and whose elements are themselves arrays, each of which is defined on $\text{Int}_2$. Currying a two-dimensional array would be like organizing a spreadsheet into a one-dimensional array of rows of the spreadsheet.
91* '''Traversing some indices in a multi-index in reverse order:''' Consider an array $A$ with domain $D_A=[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$. Fix $D_B=D_A$ and assume we're given a vector of booleans $F$ ($F$ for "flip?"). Then define $T_{BA}:D_B\to D_A$ by $i_j\to i_j$ if $F_j$ is {{#f}} and $i_j\to u_j+l_j-1-i_j$ if $F_j$ is {{#t}}. In other words, we reverse the ordering of the $j$th coordinate of $\vec i$ if and only if $F_j$ is true. $T_{BA}$ is an affine mapping from $D_B\to D_A$, which defines a new array $B$, and we can provide {{array-reverse}} for this operation. Applying {{array-reverse}} to a two-dimensional spreadsheet might reverse the order of the rows or columns (or both).
92* '''Uniformly sampling an array:''' Assume that $A$ is an array with domain $[0,u_1)\times\cdots\times[0,u_{d-1})$ (i.e., an interval all of whose lower bounds are zero). We'll also assume the existence of vector $S$ of scale factors, which are positive exact integers. Let $D_B$ be a new interval with $j$th lower bound equal to zero and $j$th upper bound equal to $\operatorname{ceiling}(u_j/S_j)$ and let $T_{BA}(\vec i)_j=i_j\times S_j$, i.e., the $j$th coordinate is scaled by $S_j$. ($D_B$ contains precisely those multi-indices that $T_{BA}$ maps into $D_A$.) Then $T_{BA}$ is an affine one-to-one mapping, and we provide {{interval-scale}} and {{array-sample}} for these operations.
93
94We make several remarks. First, all these operations could have been computed by specifying the particular mapping $T_{BA}$ explicitly, so that these routines are simply "convenience" procedures. Second, because the composition of any number of affine mappings is again affine, accessing or changing the elements of a restricted, translated, curried, permuted array is no slower than accessing or changing the elements of the original array itself. Finally, we note that by combining array currying and permuting, say, one can come up with simple expressions of powerful algorithms, such as extending one-dimensional transforms to multi-dimensional separable transforms, or quickly generating two-dimensional slices of three-dimensional image data. Examples are given below.
95
96==== Generalized arrays
97
98Bawden-style arrays are clearly useful as a programming construct, but they do not fulfill all our needs in this area. An array, as commonly understood, provides a mapping from multi-indices $(i_0,\ldots,i_{d-1})$ of exact integers in a nonempty, rectangular, $d$-dimensional interval $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ (the ''domain'' of the array) to Scheme objects. Thus, two things are necessary to specify an array: an interval and a mapping that has that interval as its domain.
99
100Since these two things are often sufficient for certain algorithms, we introduce in this SRFI a minimal set of interfaces for dealing with such arrays.
101
102Specifically, an array specifies a nonempty, multi-dimensional interval, called its ''domain'', and a mapping from this domain to Scheme objects. This mapping is called the ''getter'' of the array, accessed with the procedure {{array-getter}}; the domain of the array (more precisely, the domain of the array's getter) is accessed with the procedure {{array-domain}}.
103
104If this mapping can be changed, the array is said to be ''mutable'' and the mutation is effected by the array's ''setter'', accessed by the procedure {{array-setter}}. We call an object of this type a mutable array. Note: If an array does not have a setter, then we call it immutable even though the array's getter might not be a "pure" function, i.e., the value it returns may not depend solely on the arguments passed to the getter.
105
106In general, we leave the implementation of generalized arrays completely open. They may be defined simply by closures, or they may have hash tables or databases behind an implementation, one may read the values from a file, etc.
107
108In this SRFI, Bawden-style arrays are called ''specialized''. A specialized array is an example of a mutable array.
109
110==== Sharing generalized arrays
111
112Even if an array $A$ is not a specialized array, then it could be "shared" by specifying a new interval $D_B$ as the domain of a new array $B$ and an affine map $T_{BA}:D_B\to D_A$. Each call to $B$ would then be computed as $B(\vec i)=A(T_{BA}(\vec i))$.
113
114One could again "share" $B$, given a new interval $D_C$ as the domain of a new array $C$ and an affine transform $T_{CB}:D_C\to D_B$, and then each access $C(\vec i)=A(T_{BA}(T_{CB}(\vec i)))$. The composition $T_{BA}\circ T_{CB}:D_C\to D_A$, being itself affine, could be precomputed and stored as $T_{CA}:D_C\to D_A$, and $C(\vec i)=A(T_{CA}(\vec i))$ can be computed with the overhead of computing a single affine transformation.
115
116So, if we wanted, we could share generalized arrays with constant overhead by adding a single layer of (multi-valued) affine transformations on top of evaluating generalized arrays. Even though this could be done transparently to the user, we do not do that here; it would be a compatible extension of this SRFI to do so. We provide only the routine {{specialized-array-share}}, not a more general {{array-share}}.
117
118Certain ways of sharing generalized arrays, however, are relatively easy to code and not that expensive. If we denote {{(array-getter A)}} by {{A-getter}}, then if B is the result of {{array-extract}} applied to A, then {{(array-getter B)}} is simply {{A-getter}}. Similarly, if A is a two-dimensional array, and B is derived from A by applying the permutation $\pi((i,j))=(j,i)$, then {{(array-getter B)}} is {{(lambda (i j) (A-getter j i))}}. Translation and currying also lead to transformed arrays whose getters are relatively efficiently derived from {{A-getter}}, at least for arrays of small dimension.
119
120Thus, while we do not provide for sharing of generalized arrays for general one-to-one affine maps $T$, we do allow it for the specific functions {{array-extract}}, {{array-translate}}, {{array-permute}}, {{array-curry}}, {{array-reverse}}, {{array-tile}}, {{array-rotate}} and {{array-sample}}, and we provide relatively efficient implementations of these functions for arrays of dimension no greater than four.
121
122==== Array-map does not produce a specialized array
123
124Daniel Friedman and David Wise wrote a famous paper [[http://www.cs.indiana.edu/cgi-bin/techreports/TRNNN.cgi?trnum=TR44|CONS should not Evaluate its Arguments]]. In the spirit of that paper, our procedure {{array-map}} does not immediately produce a specialized array, but a simple immutable array, whose elements are recomputed from the arguments of {{array-map}} each time they are accessed. This immutable array can be passed on to further applications of {{array-map}} for further processing without generating the storage bodies for intermediate arrays.
125
126We provide the procedure {{array-copy}} to transform a generalized array (like that returned by {{array-map}}) to a specialized, Bawden-style array, for which accessing each element again takes $O(1)$ operations.
127
128==== Notational convention
129
130If {{A}} is an array, then we generally define {{A_}} to be {{(array-getter A)}} and {{A!}} to be {{(array-setter A)}}. The latter notation is motivated by the general Scheme convention that the names of functions that modify the contents of data structures end in {{!}}, while the notation for the getter of an array is motivated by the TeX notation for subscripts. See particularly the [[#Haar|Haar transform]] example.
131
132=== Notes
133
134* '''Relationship to [[https://docs.racket-lang.org/math/array_nonstrict.html#%28tech._nonstrict%29|nonstrict arrays]] in Racket.''' It appears that what we call simply arrays in this SRFI are called nonstrict arrays in the math/array library of Racket, which in turn was influenced by an [[https://research.microsoft.com/en-us/um/people/simonpj/papers/ndp/RArrays.pdf|array proposal for Haskell]]. Our "specialized" arrays are related to Racket's "strict" arrays.
135* '''Indexers.''' The argument {{new-domain->old-domain}} to {{specialized-array-share}} is, conceptually, a multi-valued array.
136* '''Source of function names.''' The function {{array-curry}} gets its name from the [[https://en.wikipedia.org/wiki/Currying|curry operator]] in programming—we are currying the getter of the array and keeping careful track of the domains. {{interval-projections}} can be thought of as currying the characteristic function of the interval, encapsulated here as {{interval-contains-multi-index?}}.
137* '''Choice of functions on intervals.''' The choice of functions for both arrays and intervals was motivated almost solely by what I needed for arrays.
138* '''No empty intervals.''' This SRFI considers arrays over only nonempty intervals of positive dimension. The author of this proposal acknowledges that other languages and array systems allow either zero-dimensional intervals or empty intervals of positive dimension, but prefers to leave such empty intervals as possibly compatible extensions to the current proposal.
139* '''Multi-valued arrays.''' While this SRFI restricts attention to single-valued arrays, wherein the getter of each array returns a single value, allowing multi-valued immutable arrays would be a compatible extension of this SRFI.
140* '''No low-level specialized array constructor.''' While the author of the SRFI uses mainly {{(make-array ...)}}, {{array-map}}, and {{array-copy}} to construct arrays, and while there are several other ways to construct arrays, there is no really low-level interface given for constructing specialized arrays (where one specifies a body, an indexer, etc.). It was felt that certain difficulties, some surmountable (such as checking that a given body is compatible with a given storage class) and some not (such as checking that an indexer is indeed affine), made a low-level interface less useful. At the same time, the simple {{(make-array ...)}} mechanism is so general, allowing one to specify getters and setters as general functions, as to cover nearly all needs.
141
142=== Specification
143
144; Miscellaneous Functions : [[#translation?|translation?]], [[#permutation?|permutation?]].
145; Intervals : [[#make-interval|make-interval]], [[#interval?|interval?]], [[#interval-dimension|interval-dimension]], [[#interval-lower-bound|interval-lower-bound]], [[#interval-upper-bound|interval-upper-bound]], [[#interval-lower-bounds-rarrow-list|interval-lower-bounds->list]], [[#interval-upper-bounds-rarrow-list|interval-upper-bounds->list]], [[#interval-lower-bounds-rarrow-vector|interval-lower-bounds->vector]], [[#interval-upper-bounds-rarrow-vector|interval-upper-bounds->vector]], [[#interval=|interval=]], [[#interval-volume|interval-volume]], [[#interval-subset?|interval-subset?]], [[#interval-contains-multi-index?|interval-contains-multi-index?]], [[#interval-projections|interval-projections]], [[#interval-for-each|interval-for-each]], [[#interval-dilate|interval-dilate]], [[#interval-intersect|interval-intersect]], [[#interval-translate|interval-translate]], [[#interval-permute|interval-permute]], [[#interval-rotate|interval-rotate]], [[#interval-scale|interval-scale]], [[#interval-cartesian-product|interval-cartesian-product]].
146; Storage Classes : [[#make-storage-class|make-storage-class]], [[#storage-class?|storage-class?]], [[#storage-class-getter|storage-class-getter]], [[#storage-class-setter|storage-class-setter]], [[#storage-class-checker|storage-class-checker]], [[#storage-class-maker|storage-class-maker]], [[#storage-class-copier|storage-class-copier]], [[#storage-class-length|storage-class-length]], [[#storage-class-default|storage-class-default]], [[#generic-storage-class|generic-storage-class]], [[#s8-storage-class|s8-storage-class]], [[#s16-storage-class|s16-storage-class]], [[#s32-storage-class|s32-storage-class]], [[#s64-storage-class|s64-storage-class]], [[#u1-storage-class|u1-storage-class]], [[#u8-storage-class|u8-storage-class]], [[#u16-storage-class|u16-storage-class]], [[#u32-storage-class|u32-storage-class]], [[#u64-storage-class|u64-storage-class]], [[#f8-storage-class|f8-storage-class]], [[#f16-storage-class|f16-storage-class]], [[#f32-storage-class|f32-storage-class]], [[#f64-storage-class|f64-storage-class]], [[#c64-storage-class|c64-storage-class]], [[#c128-storage-class|c128-storage-class]].
147; Arrays : [[#make-array|make-array]], [[#array?|array?]], [[#array-domain|array-domain]], [[#array-getter|array-getter]], [[#array-dimension|array-dimension]], [[#mutable-array?|mutable-array?]], [[#array-setter|array-setter]], [[#specialized-array-default-safe?|specialized-array-default-safe?]], [[#specialized-array-default-mutable?|specialized-array-default-mutable?]], [[#make-specialized-array|make-specialized-array]], [[#specialized-array?|specialized-array?]], [[#array-storage-class|array-storage-class]], [[#array-indexer|array-indexer]], [[#array-body|array-body]], [[#array-safe?|array-safe?]], [[#array-elements-in-order?|array-elements-in-order?]], [[#specialized-array-share|specialized-array-share]], [[#array-copy|array-copy]], [[#array-curry|array-curry]], [[#array-extract|array-extract]], [[#array-tile|array-tile]], [[#array-translate|array-translate]], [[#array-permute|array-permute]], [[#array-rotate|array-rotate]], [[#array-reverse|array-reverse]], [[#array-sample|array-sample]], [[#array-outer-product|array-outer-product]], [[#array-map|array-map]], [[#array-for-each|array-for-each]], [[#array-fold|array-fold]], [[#array-fold-right|array-fold-right]], [[#array-reduce|array-reduce]], [[#array-any|array-any]], [[#array-every|array-every]], [[#array-rarrow-list|array->list]], [[#list-rarrow-array|list->array]], [[#array-assign!|array-assign!]], [[#array-ref|array-ref]], [[#array-set!|array-set!]], [[#specialized-array-reshape|specialized-array-reshape]].
148
149=== Miscellaneous Functions
150
151This document refers to ''translations'' and ''permutations''. A translation is a vector of exact integers. A permutation of dimension $n$ is a vector whose entries are the exact integers $0,1,\ldots,n-1$, each occurring once, in any order.
152
153==== Procedures
154
155<procedure>(translation? object)</procedure>
156
157Returns {{#t}} if {{object}} is a translation, and {{#f}} otherwise.
158
159<procedure>(permutation? object)</procedure>
160
161Returns {{#t}} if {{object}} is a permutation, and {{#f}} otherwise.
162
163=== Intervals
164
165An interval represents the set of all multi-indices of exact integers $i_0,\ldots,i_{d-1}$ satisfying $l_0\leq i_0<u_0,\ldots,l_{d-1}\leq i_{d-1}<u_{d-1}$, where the ''lower bounds'' $l_0,\ldots,l_{d-1}$ and the ''upper bounds'' $u_0,\ldots,u_{d-1}$ are specified multi-indices of exact integers. The positive integer $d$ is the ''dimension'' of the interval. It is required that $l_0<u_0,\ldots,l_{d-1}<u_{d-1}$.
166
167Intervals are a data type distinct from other Scheme data types.
168
169==== Procedures
170
171<procedure>(make-interval arg1 #!optional arg2)</procedure>
172
173Create a new interval. {{arg1}} and {{arg2}} (if given) are nonempty vectors (of the same length) of exact integers.
174
175If {{arg2}} is not given, then the entries of {{arg1}} must be positive, and they are taken as the {{upper-bounds}} of the interval, and {{lower-bounds}} is set to a vector of the same length with exact zero entries.
176
177If {{arg2}} is given, then {{arg1}} is taken to be {{lower-bounds}} and {{arg2}} is taken to be {{upper-bounds}}, which must satisfy
178
179  (< (vector-ref lower-bounds i) (vector-ref upper-bounds i))
180for $0\leq i<{}${{(vector-length lower-bounds)}}. It is an error if {{lower-bounds}} and {{upper-bounds}} do not satisfy these conditions.
181
182<procedure>(interval? obj)</procedure>
183
184Returns {{#t}} if {{obj}} is an interval, and {{#f}} otherwise.
185
186<procedure>(interval-dimension interval)</procedure>
187
188If {{interval}} is an interval built with
189
190 (make-interval lower-bounds upper-bounds)
191then {{interval-dimension}} returns {{(vector-length lower-bounds)}}. It is an error to call {{interval-dimension}} if {{interval}} is not an interval.
192
193<procedure>(interval-lower-bound interval i)</procedure>
194
195<procedure>(interval-upper-bound interval i)</procedure>
196
197If {{interval}} is an interval built with
198
199 (make-interval lower-bounds upper-bounds)
200
201and {{i}} is an exact integer that satisfies
202
203$0 \leq i<$ {{(vector-length lower-bounds)}},
204
205then {{interval-lower-bound}} returns {{(vector-ref lower-bounds i)}} and {{interval-upper-bound}} returns {{(vector-ref upper-bounds i)}}. It is an error to call {{interval-lower-bound}} or {{interval-upper-bound}} if {{interval}} and {{i}} do not satisfy these conditions.
206
207<procedure>(interval-lower-bounds->list interval)</procedure>
208
209<procedure>(interval-upper-bounds->list interval)</procedure>
210
211If {{interval}} is an interval built with
212
213 (make-interval lower-bounds upper-bounds)
214
215then {{interval-lower-bounds->list}} returns {{(vector->list lower-bounds)}} and {{interval-upper-bounds->list}} returns {{(vector->list upper-bounds)}}. It is an error to call {{interval-lower-bounds->list}} or {{interval-upper-bounds->list}} if {{interval}} does not satisfy these conditions.
216
217<procedure>(interval-lower-bounds->vector interval)</procedure>
218
219<procedure>(interval-upper-bounds->vector interval)</procedure>
220
221If {{interval}} is an interval built with
222
223 (make-interval lower-bounds upper-bounds)
224then {{interval-lower-bounds->vector}} returns a copy of {{lower-bounds}} and {{interval-upper-bounds->vector}} returns a copy of {{upper-bounds}}. It is an error to call {{interval-lower-bounds->vector}} or {{interval-upper-bounds->vector}} if {{interval}} does not satisfy these conditions.
225
226<procedure>(interval-volume interval)</procedure>
227
228If {{interval}} is an interval built with
229
230 (make-interval lower-bounds upper-bounds)
231then, assuming the existence of {{vector-map}}, {{interval-volume}} returns
232
233 (apply * (vector->list (vector-map - upper-bounds lower-bounds)))
234It is an error to call {{interval-volume}} if {{interval}} does not satisfy this condition.
235
236<procedure>(interval= interval1 interval2)</procedure>
237
238If {{interval1}} and {{interval2}} are intervals built with
239
240 (make-interval lower-bounds1 upper-bounds1)
241and
242
243 (make-interval lower-bounds2 upper-bounds2)
244respectively, then {{interval=}} returns
245
246 (and (equal? lower-bounds1 lower-bounds2) (equal? upper-bounds1 upper-bounds2))
247It is an error to call {{interval=}} if {{interval1}} or {{interval2}} do not satisfy this condition.
248
249<procedure>(interval-subset? interval1 interval2)</procedure>
250
251If {{interval1}} and {{interval2}} are intervals of the same dimension $d$, then {{interval-subset?}} returns {{#t}} if
252
253 (>= (interval-lower-bound interval1 j) (interval-lower-bound interval2 j))
254and
255
256 (<= (interval-upper-bound interval1 j) (interval-upper-bound interval2 j))
257for all $0\leq j<d$, otherwise it returns {{#f}}. It is an error if the arguments do not satisfy these conditions.
258
259<procedure>(interval-contains-multi-index? interval index-0 index-1 ...)</procedure>
260
261If {{interval}} is an interval with dimension $d$ and {{index-0}}, {{index-1}}, ..., is a multi-index of length $d$, then {{interval-contains-multi-index?}} returns {{#t}} if
262
263{{(interval-lower-bound interval j)}} $\leq$ {{index-j}} $<$ {{(interval-upper-bound interval j)}}
264
265for $0\leq j < d$, and {{#f}} otherwise.
266
267It is an error to call {{interval-contains-multi-index?}} if {{interval}} and {{index-0}},..., do not satisfy this condition.
268
269<procedure>(interval-projections interval right-dimension)</procedure>
270
271Conceptually, {{interval-projections}} takes a $d$-dimensional interval $[l_0,u_0)\times [l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ and splits it into two parts
272
273$[l_0,u_0)\times\cdots\times[l_{d-\text{right-dimension}-1},u_{d-\text{right-dimension}-1})$
274
275and
276
277$[l_{d-\text{right-dimension}},u_{d-\text{right-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$
278
279This function, the inverse of Cartesian products or cross products of intervals, is used to keep track of the domains of curried arrays.
280
281More precisely, if {{interval}} is an interval and {{right-dimension}} is an exact integer that satisfies {{0 < right-dimension < d}} then {{interval-projections}} returns two intervals:
282
283<enscript highlight="scheme">(values
284 (make-interval
285  (vector (interval-lower-bound interval 0)
286          ...
287          (interval-lower-bound interval
288                                (- d right-dimension 1)))
289  (vector (interval-upper-bound interval 0)
290          ...
291          (interval-upper-bound interval
292                                (- d right-dimension 1))))
293 (make-interval
294  (vector (interval-lower-bound interval
295                                (- d right-dimension))
296          ...
297          (interval-lower-bound interval
298                                (- d 1)))
299  (vector (interval-upper-bound interval
300                                (- d right-dimension))
301          ...
302          (interval-upper-bound interval
303                                (- d 1)))))</enscript>
304It is an error to call {{interval-projections}} if its arguments do not satisfy these conditions.
305
306<procedure>(interval-for-each f interval)</procedure>
307
308This routine assumes that {{interval}} is an interval and {{f}} is a routine whose domain includes elements of {{interval}}. It is an error to call {{interval-for-each}} if {{interval}} and {{f}} do not satisfy these conditions.
309
310{{interval-for-each}} calls {{f}} with each multi-index of {{interval}} as arguments, all in lexicographical order.
311
312<procedure>(interval-dilate interval lower-diffs upper-diffs)</procedure>
313
314If {{interval}} is an interval with lower bounds $\ell_0,\dots,\ell_{d-1}$ and upper bounds $u_0,\dots,u_{d-1}$, and {{lower-diffs}} is a vector of exact integers $L_0,\dots,L_{d-1}$ and {{upper-diffs}} is a vector of exact integers $U_0,\dots,U_{d-1}$, then {{interval-dilate}} returns a new interval with lower bounds $\ell_0+L_0,\dots,\ell_{d-1}+L_{d-1}$ and upper bounds $u_0+U_0,\dots,u_{d-1}+U_{d-1}$, as long as this is a nonempty interval. It is an error if the arguments do not satisfy these conditions.
315
316Examples:
317
318<enscript highlight="scheme">(interval=
319 (interval-dilate (make-interval '#(100 100))
320                  '#(1 1) '#(1 1))
321 (make-interval '#(1 1) '#(101 101))) => #t
322(interval=
323 (interval-dilate (make-interval '#(100 100))
324                  '#(-1 -1) '#(1 1))
325 (make-interval '#(-1 -1) '#(101 101))) => #t
326(interval=
327 (interval-dilate (make-interval '#(100 100))
328                  '#(0 0) '#(-50 -50))
329 (make-interval '#(50 50))) => #t
330(interval-dilate
331 (make-interval '#(100 100))
332 '#(0 0) '#(-500 -50)) => error</enscript>
333<procedure>(interval-intersect interval-1 interval-2 ...)</procedure>
334
335If all the arguments are intervals of the same dimension and they have a nonempty intersection, then {{interval-intersect}} returns that intersection; otherwise it returns {{#f}}.
336
337It is an error if the arguments are not all intervals with the same dimension.
338
339<procedure>(interval-translate interval translation)</procedure>
340
341If {{interval}} is an interval with lower bounds $\ell_0,\dots,\ell_{d-1}$ and upper bounds $u_0,\dots,u_{d-1}$, and {{translation}} is a translation with entries $T_0,\dots,T_{d-1}$ , then {{interval-translate}} returns a new interval with lower bounds $\ell_0+T_0,\dots,\ell_{d-1}+T_{d-1}$ and upper bounds $u_0+T_0,\dots,u_{d-1}+T_{d-1}$. It is an error if the arguments do not satisfy these conditions.
342
343One could define {{(interval-translate interval translation)}} by {{(interval-dilate interval translation translation)}}.
344
345<procedure>(interval-permute interval permutation)</procedure>
346
347The argument {{interval}} must be an interval, and the argument {{permutation}} must be a valid permutation with the same dimension as {{interval}}. It is an error if the arguments do not satisfy these conditions.
348
349Heuristically, this function returns a new interval whose axes have been permuted in a way consistent with {{permutation}}. But we have to say how the entries of {{permutation}} are associated with the new interval.
350
351We have chosen the following convention: If the permutation is $(\pi_0,\ldots,\pi_{d-1})$, and the argument interval represents the cross product $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$, then the result represents the cross product $[l_{\pi_0},u_{\pi_0})\times[l_{\pi_1},u_{\pi_1})\times\cdots\times[l_{\pi_{d-1}},u_{\pi_{d-1}})$.
352
353For example, if the argument interval represents $[0,4)\times[0,8)\times[0,21)\times [0,16)$ and the permutation is {{#(3 0 1 2)}}, then the result of {{(interval-permute interval permutation)}} will be the representation of $[0,16)\times [0,4)\times[0,8)\times[0,21)$.
354
355<procedure>(interval-rotate interval dim)</procedure>
356
357Informally, {{(interval-rotate interval dim)}} rotates the axes of {{interval}} {{dim}} places to the left.
358
359More precisely, {{(interval-rotate interval dim)}} assumes that {{interval}} is an interval and {{dim}} is an exact integer between 0 (inclusive) and {{(interval-dimension interval)}} (exclusive). It computes the permutation {{(vector dim ... (- (interval-dimension interval) 1) 0 ... (- dim 1))}} (unless {{dim}} is zero, in which case it constructs the identity permutation) and returns {{(interval-permute interval permutation)}}. It is an error if the arguments do not satisfy these conditions.
360
361<procedure>(interval-scale interval scales)</procedure>
362
363If {{interval}} is a $d$-dimensional interval $[0,u_1)\times\cdots\times[0,u_{d-1})$ with all lower bounds zero, and {{scales}} is a length-$d$ vector of positive exact integers, which we'll denote by $\vec s$, then {{interval-scale}} returns the interval $[0,\operatorname{ceiling}(u_1/s_1))\times\cdots\times[0,\operatorname{ceiling}(u_{d-1}/s_{d-1}))$.
364
365It is an error if {{interval}} and {{scales}} do not satisfy this condition.
366
367<procedure>(interval-cartesian-product interval . intervals)</procedure>
368
369Implements the Cartesian product of the intervals in {{(cons interval intervals)}}. Returns
370
371<enscript highlight="scheme">(make-interval (list->vector (apply append (map array-lower-bounds->list (cons interval intervals))))
372               (list->vector (apply append (map array-upper-bounds->list (cons interval intervals)))))</enscript>
373It is an error if any argument is not an interval.
374
375=== Storage classes
376
377Conceptually, a storage-class is a set of functions to manage the backing store of a specialized array. The functions allow one to make a backing store, to get values from the store and to set new values, to return the length of the store, and to specify a default value for initial elements of the backing store. Typically, a backing store is a (heterogeneous or homogeneous) vector. A storage-class has a type distinct from other Scheme types.
378
379==== Procedures
380
381<procedure>(make-storage-class getter setter checker maker copier length default)</procedure>
382
383Here we assume the following relationships between the arguments of {{make-storage-class}}. Assume that the "elements" of the backing store are of some "type", either heterogeneous (all Scheme types) or homogeneous (of some restricted type).
384
385* {{(maker n value)}} returns a linearly addressed object containing {{n}} elements of value {{value}}.
386* {{copier}} may be #f or a procedure; if a procedure then if {{to}} and {{from}} were created by {{maker}}, then {{(copier to at from start end)}} copies elements from {{from}} beginning at {{start}} (inclusive) and ending at {{end}} (exclusive) to {{to}} beginning at {{at}}. It is assumed that all the indices involved are within the domain of {{from}} and {{to}}, as needed. The order in which the elements are copied is unspecified.
387* If {{v}} is an object created by {{(maker n value)}} and 0 <== {{i}} < {{n}}, then {{(getter v i)}} returns the current value of the {{i}}'th element of {{v}}, and {{(checker (getter v i))> #t}}.
388* If {{v}} is an object created by {{(maker n value)}}, 0 <== {{i}} < {{n}}, and {{(checker val)> #t}}, then {{(setter v i val)}} sets the value of the {{i}}'th element of {{v}} to {{val}}.
389* If {{v}} is an object created by {{(maker n value)}} then {{(length v)}} returns {{n}}.
390
391If the arguments do not satisfy these conditions, then it is an error to call {{make-storage-class}}.
392
393Note that we assume that {{getter}} and {{setter}} generally take ''O''(1) time to execute.
394
395<procedure>(storage-class? m)</procedure>
396
397Returns {{#t}} if {{m}} is a storage class, and {{#f}} otherwise.
398
399<procedure>(storage-class-getter m)</procedure>
400
401<procedure>(storage-class-setter m)</procedure>
402
403<procedure>(storage-class-checker m)</procedure>
404
405<procedure>(storage-class-maker m)</procedure>
406
407<procedure>(storage-class-copier m)</procedure>
408
409<procedure>(storage-class-length m)</procedure>
410
411<procedure>(storage-class-default m)</procedure>
412
413If {{m}} is an object created by
414
415 (make-storage-class getter setter checker maker copier length default)
416
417then {{storage-class-getter}} returns {{getter}}, {{storage-class-setter}} returns {{setter}}, {{storage-class-checker}} returns {{checker}}, {{storage-class-maker}} returns {{maker}}, {{storage-class-copier}} returns {{copier}}, {{storage-class-length}} returns {{length}}, and {{storage-class-default}} returns {{default}}. Otherwise, it is an error to call any of these routines.
418
419==== Global Variables
420
421<constant>generic-storage-class</constant>
422
423<constant>s8-storage-class</constant>
424
425<constant>s16-storage-class</constant>
426
427<constant>s32-storage-class</constant>
428
429<constant>s64-storage-class</constant>
430
431<constant>u1-storage-class</constant>
432
433<constant>u8-storage-class</constant>
434
435<constant>u16-storage-class</constant>
436
437<constant>u32-storage-class</constant>
438
439<constant>u64-storage-class</constant>
440
441<constant>f8-storage-class</constant>
442
443<constant>f16-storage-class</constant>
444
445<constant>f32-storage-class</constant>
446
447<constant>f64-storage-class</constant>
448
449<constant>c64-storage-class</constant>
450
451<constant>c128-storage-class</constant>
452
453{{generic-storage-class}} is defined as if by
454
455<enscript highlight="scheme">(define generic-storage-class
456  (make-storage-class vector-ref
457                      vector-set!
458                      (lambda (arg) #t)
459                      make-vector
460                      vector-copy!
461                      vector-length
462                      #f))</enscript>
463Implementations shall define {{sX-storage-class}} for {{X}}=8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between -2<sup>{{X}}-1</sup> and 2<sup>{{X}}-1</sup>-1 inclusive), {{uX-storage-class}} for {{X}}=1, 8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between 0 and 2<sup>{{X}}</sup>-1 inclusive), {{fX-storage-class}} for {{X}}= 8, 16, 32, and 64 (which have default value 0.0 and manipulate 8-, 16-, 32-, and 64-bit floating-point numbers), and {{cX-storage-class}} for {{X}}= 64 and 128 (which have default value 0.0+0.0i and manipulate complex numbers with, respectively, 32- and 64-bit floating-point numbers as real and imaginary parts).
464
465Implementations with an appropriate homogeneous vector type should define the associated global variable using {{make-storage-class}}, otherwise they shall define the associated global variable to {{#f}}.
466
467=== Arrays
468
469Arrays are a data type distinct from other Scheme data types.
470
471==== Procedures
472
473<procedure>(make-array interval getter [ setter ])</procedure>
474
475Assume first that the optional argument {{setter}} is not given.
476
477If {{interval}} is an interval and {{getter}} is a function from {{interval}} to Scheme objects, then {{make-array}} returns an array with domain {{interval}} and getter {{getter}}.
478
479It is an error to call {{make-array}} if {{interval}} and {{getter}} do not satisfy these conditions.
480
481If now {{setter}} is specified, assume that it is a procedure such that getter and setter satisfy: If
482
483{{(i1,...,in)}} $\neq$ {{(j1,...,jn)}}
484
485are elements of {{interval}} and
486
487 (getter j1 ... jn) => x
488
489then "after"
490
491  (setter v i1 ... in)
492
493we have
494
495 (getter j1 ... jn) => x
496
497and
498
499 (getter i1,...,in) => v
500
501Then {{make-array}} builds a mutable array with domain {{interval}}, getter {{getter}}, and setter {{setter}}. It is an error to call {{make-array}} if its arguments do not satisfy these conditions.
502
503Example:
504
505<enscript highlight="scheme">  (define a (make-array (make-interval '#(1 1) '#(11 11))
506                        (lambda (i j)
507                          (if (= i j)
508                              1
509                              0))))</enscript>
510defines an array for which {{(array-getter a)}} returns 1 when i=j and 0 otherwise.
511
512Example:
513
514<enscript highlight="scheme">(define a   ;; a sparse array
515  (let ((domain
516         (make-interval '#(1000000 1000000)))
517        (sparse-rows
518         (make-vector 1000000 '())))
519    (make-array
520     domain
521     (lambda (i j)
522       (cond ((assv j (vector-ref sparse-rows i))
523              => cdr)
524             (else
525              0.0)))
526     (lambda (v i j)
527       (cond ((assv j (vector-ref sparse-rows i))
528              => (lambda (pair)
529                   (set-cdr! pair v)))
530             (else
531              (vector-set!
532               sparse-rows
533               i
534               (cons (cons j v)
535                     (vector-ref sparse-rows i)))))))))
536
537(define a_ (array-getter a))
538(define a! (array-setter a))
539
540(a_ 12345 6789)  => 0.
541(a_ 0 0) => 0.
542(a! 1.0 0 0) => undefined
543(a_ 12345 6789)  => 0.
544(a_ 0 0) => 1.</enscript>
545<procedure>(array? obj)</procedure>
546
547Returns {{#t}} if {{obj}} is an array and {{#f}} otherwise.
548
549<procedure>(array-domain array)</procedure>
550
551<procedure>(array-getter array)</procedure>
552
553If {{array}} is an array built by
554
555 (make-array interval getter [setter])
556(with or without the optional {{setter}} argument) then {{array-domain}} returns {{interval}} and {{array-getter}} returns {{getter}}. It is an error to call {{array-domain}} or {{array-getter}} if {{array}} is not an array.
557
558Example:
559
560<enscript highlight="scheme">(define a (make-array (make-interval '#(1 1) '#(11 11))
561                      (lambda (i j)
562                        (if (= i j)
563                            1
564                            0))))
565(define a_ (array-getter a))
566
567(a_ 3 3) => 1
568(a_ 2 3) => 0
569(a_ 11 0) => is an error</enscript>
570<procedure>(array-dimension array)</procedure>
571
572Shorthand for {{(interval-dimension (array-domain array))}}. It is an error to call this function if {{array}} is not an array.
573
574<procedure>(mutable-array? obj)</procedure>
575
576Returns {{#t}} if {{obj}} is a mutable array and {{#f}} otherwise.
577
578<procedure>(array-setter array)</procedure>
579
580If {{array}} is an array built by
581
582 (make-array interval getter setter)
583then {{array-setter}} returns {{setter}}. It is an error to call {{array-setter}} if {{array}} is not a mutable array.
584
585<parameter>(specialized-array-default-safe?)</parameter>
586
587With no argument, returns {{#t}} if newly constructed specialized arrays check the arguments of setters and getters by default, and {{#f}} otherwise.
588
589If {{bool}} is {{#t}} then the next call to {{specialized-array-default-safe?}} will return {{#t}}; if {{bool}} is {{#f}} then the next call to {{specialized-array-default-safe?}} will return {{#f}}; otherwise it is an error.
590
591<parameter>(specialized-array-default-mutable? [ bool ])</parameter>
592
593With no argument, returns {{#t}} if newly constructed specialized arrays are mutable by default, and {{#f}} otherwise.
594
595If {{bool}} is {{#t}} then the next call to {{specialized-array-default-mutable?}} will return {{#t}}; if {{bool}} is {{#f}} then the next call to {{specialized-array-default-mutable?}} will return {{#f}}; otherwise it is an error.
596
597<procedure>(make-specialized-array interval [ storage-class generic-storage-class ] [ safe? (specialized-array-default-safe?) ])</procedure>
598
599Constructs a mutable specialized array from its arguments.
600
601{{interval}} must be given as a nonempty interval. If given, {{storage-class}} must be a storage class; if it is not given it defaults to {{generic-storage-class}}. If given, {{safe?}} must be a boolean; if it is not given it defaults to the current value of {{(specialized-array-default-safe?)}}.
602
603The body of the result is constructed as
604
605<enscript highlight="scheme">  ((storage-class-maker storage-class)
606   (interval-volume interval)
607   (storage-class-default storage-class))
608  </enscript>
609The indexer of the resulting array is constructed as the lexicographical mapping of {{interval}} onto the interval {{[0,(interval-volume interval))}}.
610
611If {{safe}} is {{#t}}, then the arguments of the getter and setter (including the value to be stored) of the resulting array are always checked for correctness.
612
613After correctness checking (if needed), {{(array-getter array)}} is defined simply as
614
615<enscript highlight="scheme">  (lambda multi-index
616    ((storage-class-getter storage-class)
617     (array-body array)
618     (apply (array-indexer array) multi-index)))
619  </enscript>
620and {{(array-setter array)}} is defined as
621
622<enscript highlight="scheme">  (lambda (val . multi-index)
623    ((storage-class-setter storage-class)
624     (array-body array)
625     (apply (array-indexer array) multi-index)
626     val))
627  </enscript>
628It is an error if the arguments of {{make-specialized-array}} do not satisfy these conditions.
629
630'''Examples.''' A simple array that can hold any type of element can be defined with {{(make-specialized-array (make-interval '#(3 3)))}}. If you find that you're using a lot of unsafe arrays of unsigned 16-bit integers, one could define
631
632<enscript highlight="scheme">  (define (make-u16-array interval)
633    (make-specialized-array interval u16-storage-class #f))</enscript>
634and then simply call, e.g., {{(make-u16-array (make-interval '#(3 3)))}}.
635
636<procedure>(specialized-array? obj)</procedure>
637
638Returns {{#t}} if {{obj}} is a specialized-array, and {{#f}} otherwise. A specialized-array is an array.
639
640<procedure>(array-storage-class array)</procedure>
641
642<procedure>(array-indexer array)</procedure>
643
644<procedure>(array-body array)</procedure>
645
646<procedure>(array-safe? array)</procedure>
647
648{{array-storage-class}} returns the storage-class of {{array}}. {{array-safe?}} is true if and only if the arguments of {{(array-getter array)}} and {{(array-setter array)}} (including the value to be stored in the array) are checked for correctness.
649
650{{(array-body array)}} is a linearly indexed, vector-like object (e.g., a vector, string, u8vector, etc.) indexed from 0.
651
652{{(array-indexer array)}} is assumed to be a one-to-one, but not necessarily onto, affine mapping from {{(array-domain array)}} into the indexing domain of {{(array-body array)}}.
653
654Please see [[#make-specialized-array|{{make-specialized-array}}]] for how {{(array-body array)}}, etc., are used.
655
656It is an error to call any of these routines if {{array}} is not a specialized array.
657
658<procedure>(array-elements-in-order? A)</procedure>
659
660Assumes that {{A}} is a specialized array, in which case it returns {{#t}} if the elements of {{A}} are in order and stored adjacently in {{(array-body A)}} and {{#f}} otherwise.
661
662It is an error if {{A}} is not a specialized array.
663
664<procedure>(specialized-array-share array new-domain new-domain->old-domain)</procedure>
665
666Constructs a new specialized array that shares the body of the specialized array {{array}}. Returns an object that is behaviorally equivalent to a specialized array with the following fields:
667
668<enscript highlight="scheme">domain:        new-domain
669storage-class: (array-storage-class array)
670body:          (array-body array)
671indexer:       (lambda multi-index
672                 (call-with-values
673                     (lambda ()
674                       (apply new-domain->old-domain
675                              multi-index))
676                   (array-indexer array)))</enscript>
677The resulting array inherits its safety and mutability from {{array}}.
678
679Note: It is assumed that the affine structure of the composition of {{new-domain->old-domain}} and {{(array-indexer array)}} will be used to simplify:
680
681<enscript highlight="scheme">  (lambda multi-index
682    (call-with-values
683        (lambda ()
684          (apply new-domain->old-domain multi-index))
685      (array-indexer array)))</enscript>
686It is an error if {{array}} is not a specialized array, or if {{new-domain}} is not an interval, or if {{new-domain->old-domain}} is not a one-to-one affine mapping from {{new-domain}} to {{(array-domain array)}}.
687
688'''Example:''' One can apply a "shearing" operation to an array as follows:
689
690<enscript highlight="scheme">(define a
691  (array-copy
692   (make-array (make-interval '#(5 10))
693               list)))
694(define b
695  (specialized-array-share
696   a
697   (make-interval '#(5 5))
698   (lambda (i j)
699     (values i (+ i j)))))
700;; Print the "rows" of b
701(array-for-each (lambda (row)
702                  (pretty-print (array->list row)))
703                (array-curry b 1))
704
705;; which prints
706;; ((0 0) (0 1) (0 2) (0 3) (0 4))
707;; ((1 1) (1 2) (1 3) (1 4) (1 5))
708;; ((2 2) (2 3) (2 4) (2 5) (2 6))
709;; ((3 3) (3 4) (3 5) (3 6) (3 7))
710;; ((4 4) (4 5) (4 6) (4 7) (4 8))</enscript>
711This "shearing" operation cannot be achieved by combining the procedures {{array-extract}}, {{array-translate}}, {{array-permute}}, {{array-translate}}, {{array-curry}}, {{array-reverse}}, and {{array-sample}}.
712
713<procedure>(array-copy array [ result-storage-class generic-storage-class ] [ new-domain #f ] [ mutable? (specialized-array-default-mutable?) ] [ safe? (specialized-array-default-safe?) ])</procedure>
714
715Assumes that {{array}} is an array, {{result-storage-class}} is a storage class that can manipulate all the elements of {{array}}, {{new-domain}} is either {{#f}} or an interval with the same volume as {{(array-domain array)}}, and {{mutable?}} and {{safe?}} are booleans.
716
717If {{new-domain}} is {{#f}}, then it is set to {{(array-domain array)}}.
718
719The specialized array returned by {{array-copy}} can be defined conceptually by:
720
721<enscript highlight="scheme">(list->array (array->list array)
722             new-domain
723             result-storage-class
724             mutable?
725             safe?)</enscript>
726It is an error if the arguments do not satisfy these conditions.
727
728'''Note:''' If {{new-domain}} is not the same as {{(array-domain array)}}, one can think of the resulting array as a ''reshaped'' version of {{array}}.
729
730<procedure>(array-curry array inner-dimension)</procedure>
731
732If {{array}} is an array whose domain is an interval $[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$, and {{inner-dimension}} is an exact integer strictly between $0$ and $d$, then {{array-curry}} returns an immutable array with domain $[l_0,u_0)\times\cdots\times[l_{d-\text{inner-dimension}-1},u_{d-\text{inner-dimension}-1})$, each of whose entries is in itself an array with domain $[l_{d-\text{inner-dimension}},u_{d-\text{inner-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$.
733
734For example, if {{A}} and {{B}} are defined by
735
736<enscript highlight="scheme">(define interval (make-interval '#(10 10 10 10)))
737(define A (make-array interval list))
738(define B (array-curry A 1))
739
740(define A_ (array-getter A))
741(define B_ (array-getter B))
742  </enscript>
743so
744
745 (A_ i j k l) => (list i j k l)
746then {{B}} is an immutable array with domain {{(make-interval '#(10 10 10))}}, each of whose elements is itself an (immutable) array and
747
748<enscript highlight="scheme">(equal?
749 (A_ i j k l)
750 ((array-getter (B_ i j k)) l)) => #t</enscript>
751for all multi-indices {{i j k l}} in {{interval}}.
752
753The subarrays are immutable, mutable, or specialized according to whether the array argument is immutable, mutable, or specialized.
754
755More precisely, if
756
757 0 < inner-dimension < (interval-dimension (array-domain array))
758then {{array-curry}} returns a result as follows.
759
760If the input array is specialized, then array-curry returns
761
762<enscript highlight="scheme">(call-with-values
763    (lambda () (interval-projections (array-domain array)
764                                     inner-dimension))
765  (lambda (outer-interval inner-interval)
766    (make-array
767     outer-interval
768     (lambda outer-multi-index
769       (specialized-array-share
770        array
771        inner-interval
772        (lambda inner-multi-index
773          (apply values
774                 (append outer-multi-index
775                         inner-multi-index))))))))</enscript>
776Otherwise, if the input array is mutable, then array-curry returns
777
778<enscript highlight="scheme">(call-with-values
779    (lambda () (interval-projections (array-domain array)
780                                     inner-dimension))
781  (lambda (outer-interval inner-interval)
782    (make-array
783     outer-interval
784     (lambda outer-multi-index
785       (make-array
786        inner-interval
787        (lambda inner-multi-index
788          (apply (array-getter array)
789                 (append outer-multi-index
790                         inner-multi-index)))
791        (lambda (v . inner-multi-index)
792          (apply (array-setter array)
793                 v
794                 (append outer-multi-index
795                         inner-multi-index))))))))</enscript>
796Otherwise, array-curry returns
797
798<enscript highlight="scheme">(call-with-values
799    (lambda () (interval-projections (array-domain array)
800                                     inner-dimension))
801  (lambda (outer-interval inner-interval)
802    (make-array
803     outer-interval
804     (lambda outer-multi-index
805       (make-array
806        inner-interval
807        (lambda inner-multi-index
808          (apply (array-getter array)
809                 (append outer-multi-index
810                         inner-multi-index))))))))</enscript>
811It is an error to call {{array-curry}} if its arguments do not satisfy these conditions.
812
813If {{array}} is a specialized array, the subarrays of the result inherit their safety and mutability from {{array}}.
814
815'''Note:''' Let's denote by {{B}} the result of {{(array-curry A k)}}. While the result of calling {{(array-getter B)}} is an immutable, mutable, or specialized array according to whether {{A}} itself is immutable, mutable, or specialized, {{B}} is always an immutable array, where {{(array-getter B)}}, which returns an array, is computed anew for each call. If {{(array-getter B)}} will be called multiple times with the same arguments, it may be useful to store these results in a specialized array for fast repeated access.
816
817Please see the note in the discussion of [[#array-tile|array-tile]].
818
819Example:
820
821<enscript highlight="scheme">(define a (make-array (make-interval '#(10 10))
822                      list))
823(define a_ (array-getter a))
824(a_ 3 4)  => (3 4)
825(define curried-a (array-curry a 1))
826(define curried-a_ (array-getter curried-a))
827((array-getter (curried-a_ 3)) 4)
828                    => (3 4)</enscript>
829<procedure>(array-extract array new-domain)</procedure>
830
831Returns a new array with the same getter (and setter, if appropriate) of the first argument, defined on the second argument.
832
833Assumes that {{array}} is an array and {{new-domain}} is an interval that is a sub-interval of {{(array-domain array)}}. If {{array}} is a specialized array, then returns
834
835<enscript highlight="scheme">  (specialized-array-share array
836                           new-domain
837                           values)
838  </enscript>
839Otherwise, if {{array}} is a mutable array, then {{array-extract}} returns
840
841<enscript highlight="scheme">  (make-array new-domain
842              (array-getter array)
843              (array-setter array))
844</enscript>
845Finally, if {{array}} is an immutable array, then {{array-extract}} returns
846
847<enscript highlight="scheme">  (make-array new-domain
848              (array-getter array))</enscript>
849It is an error if the arguments of {{array-extract}} do not satisfy these conditions.
850
851If {{array}} is a specialized array, the resulting array inherits its mutability and safety from {{array}}.
852
853<procedure>(array-tile A S)</procedure>
854
855Assume that {{A}} is an array and {{S}} is a vector of positive, exact integers. The routine {{array-tile}} returns a new immutable array $T$, each entry of which is a subarray of {{A}} whose domain has sidelengths given (mostly) by the entries of {{S}}. These subarrays completely "tile" {{A}}, in the sense that every entry in {{A}} is an entry of precisely one entry of the result $T$.
856
857More formally, if {{S}} is the vector $(s_0,\ldots,s_{d-1})$, and the domain of {{A}} is the interval $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$, then $T$ is an immutable array with all lower bounds zero and upper bounds given by $$ \operatorname{ceiling}((u_0-l_0)/s_0),\ldots,\operatorname{ceiling}((u_{d-1}-l_{d-1})/s_{d-1}). $$ The $i_0,\ldots,i_{d-1}$ entry of $T$ is {{(array-extract A D_i)}} with the interval {{D_i}} given by $$ [l_0+i_0*s_0,\min(l_0+(i_0+1)s_0,u_0))\times\cdots\times[l_{d-1}+i_{d-1}*s_{d-1},\min(l_{d-1}+(i_{d-1}+1)s_{d-1},u_{d-1})). $$ (The "minimum" operators are necessary if $u_j-l_j$ is not divisible by $s_j$.) Thus, each entry of $T$ will be a specialized, mutable, or immutable array, depending on the type of the input array {{A}}.
858
859It is an error if the arguments of {{array-tile}} do not satisfy these conditions.
860
861If {{A}} is a specialized array, the subarrays of the result inherit safety and mutability from {{A}}.
862
863'''Note:''' The routines {{array-tile}} and {{array-curry}} both decompose an array into subarrays, but in different ways. For example, if {{A}} is defined as {{(make-array (make-interval '#(10 10)) list)}}, then {{(array-tile A '#(1 10))}} returns an array with domain {{(make-interval '#(10 1))}}, each element of which is an array with domain {{(make-interval '#(1 10))}} (i.e., a two-dimensional array whose elements are two-dimensional arrays), while {{(array-curry A 1)}} returns an array with domain {{(make-interval '#(10))}}, each element of which has domain {{(make-interval '#(10))}} (i.e., a one-dimensional array whose elements are one-dimensional arrays).
864
865<procedure>(array-translate array translation)</procedure>
866
867Assumes that {{array}} is a valid array, {{translation}} is a valid translation, and that the dimensions of the array and the translation are the same. The resulting array will have domain {{(interval-translate (array-domain array) translation)}}.
868
869If {{array}} is a specialized array, returns a new specialized array
870
871<enscript highlight="scheme">(specialized-array-share
872 array
873 (interval-translate (array-domain array)
874                     translation)
875 (lambda multi-index
876   (apply values
877          (map -
878               multi-index
879               (vector->list translation)))))</enscript>
880that shares the body of {{array}}, as well as inheriting its safety and mutability.
881
882If {{array}} is not a specialized array but is a mutable array, returns a new mutable array
883
884<enscript highlight="scheme">(make-array
885 (interval-translate (array-domain array)
886                     translation)
887 (lambda multi-index
888   (apply (array-getter array)
889          (map -
890               multi-index
891               (vector->list translation))))
892 (lambda (val . multi-index)
893   (apply (array-setter array)
894          val
895          (map -
896               multi-index
897               (vector->list translation)))))
898 </enscript>
899that employs the same getter and setter as the original array argument.
900
901If {{array}} is not a mutable array, returns a new array
902
903<enscript highlight="scheme">(make-array
904 (interval-translate (array-domain array)
905                     translation)
906 (lambda multi-index
907   (apply (array-getter array)
908          (map - multi-index (vector->list translation)))))</enscript>
909that employs the same getter as the original array.
910
911It is an error if the arguments do not satisfy these conditions.
912
913<procedure>(array-permute array permutation)</procedure>
914
915Assumes that {{array}} is a valid array, {{permutation}} is a valid permutation, and that the dimensions of the array and the permutation are the same. The resulting array will have domain {{(interval-permute (array-domain array) permutation)}}.
916
917We begin with an example. Assume that the domain of {{array}} is represented by the interval $[0,4)\times[0,8)\times[0,21)\times [0,16)$, as in the example for {{interval-permute}}, and the permutation is {{#(3 0 1 2)}}. Then the domain of the new array is the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$.
918
919So the multi-index argument of the {{getter}} of the result of {{array-permute}} must lie in the new domain of the array, the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$. So if we define {{old-getter}} as {{(array-getter array)}}, the definition of the new array must be in fact
920
921<enscript highlight="scheme">(make-array (interval-permute (array-domain array)
922                              '#(3 0 1 2))
923            (lambda (l i j k)
924              (old-getter i j k l)))</enscript>
925So you see that if the first argument if the new getter is in $[0,16)$, then indeed the fourth argument of {{old-getter}} is also in $[0,16)$, as it should be. This is a subtlety that I don't see how to overcome. It is the listing of the arguments of the new getter, the {{lambda}}, that must be permuted.
926
927Mathematically, we can define $\pi^{-1}$, the inverse of a permutation $\pi$, such that $\pi^{-1}$ composed with $\pi$ gives the identity permutation. Then the getter of the new array is, in pseudo-code, {{(lambda multi-index (apply old-getter (}}$\pi^{-1}${{ multi-index)))}}. We have assumed that $\pi^{-1}$ takes a list as an argument and returns a list as a result.
928
929Employing this same pseudo-code, if {{array}} is a specialized array and we denote the permutation by $\pi$, then {{array-permute}} returns the new specialized array
930
931<enscript highlight="scheme">(specialized-array-share array
932                         (interval-permute (array-domain array) π)
933                         (lambda multi-index
934                           (apply values (π-1 multi-index))))</enscript>
935The resulting array shares the body of {{array}}, as well as its safety and mutability.
936
937Again employing this same pseudo-code, if {{array}} is not a specialized array, but is a mutable-array, then {{array-permute}} returns the new mutable
938
939<enscript highlight="scheme">(make-array (interval-permute (array-domain array) π)
940            (lambda multi-index
941              (apply (array-getter array)
942                     (π-1 multi-index)))
943            (lambda (val . multi-index)
944              (apply (array-setter array)
945                     val
946                     (π-1 multi-index))))</enscript>
947which employs the setter and the getter of the argument to {{array-permute}}.
948
949Finally, if {{array}} is not a mutable array, then {{array-permute}} returns
950
951<enscript highlight="scheme">(make-array (interval-permute (array-domain array) π)
952            (lambda multi-index
953              (apply (array-getter array)
954                     (π-1 multi-index))))</enscript>
955It is an error to call {{array-permute}} if its arguments do not satisfy these conditions.
956
957<procedure>(array-rotate array dim)</procedure>
958
959Informally, {{(array-rotate array dim)}} rotates the axes of {{array}} {{dim}} places to the left.
960
961More precisely, {{(array-rotate array dim)}} assumes that {{array}} is an array and {{dim}} is an exact integer between 0 (inclusive) and {{(array-dimension array)}} (exclusive). It computes the permutation {{(vector dim ... (- (array-dimension array) 1) 0 ... (- dim 1))}} (unless {{dim}} is zero, in which case it constructs the identity permutation) and returns {{(array-permute array permutation)}}. It is an error if the arguments do not satisfy these conditions.
962
963<procedure>(array-reverse array #!optional flip?)</procedure>
964
965We assume that {{array}} is an array and {{flip?}}, if given, is a vector of booleans whose length is the same as the dimension of {{array}}. If {{flip?}} is not given, it is set to a vector with length the same as the dimension of {{array}}, all of whose elements are {{#t}}.
966
967{{array-reverse}} returns a new array that is specialized, mutable, or immutable according to whether {{array}} is specialized, mutable, or immutable, respectively. Informally, if {{(vector-ref flip? k)}} is true, then the ordering of multi-indices in the k'th coordinate direction is reversed, and is left undisturbed otherwise.
968
969More formally, we introduce the function
970
971<enscript highlight="scheme">(define flip-multi-index
972  (let* ((domain (array-domain array))
973         (lowers (interval-lower-bounds->list domain))
974         (uppers (interval-upper-bounds->list domain)))
975    (lambda (multi-index)
976      (map (lambda (i_k flip?_k l_k u_k)
977             (if flip?_k
978                 (- (+ l_k u_k -1) i_k)
979                 i_k))
980           multi-index
981           (vector->list flip?)
982           lowers
983           uppers))))</enscript>
984Then if {{array}} is specialized, {{array-reverse}} returns
985
986<enscript highlight="scheme">(specialized-array-share
987 array
988 domain
989 (lambda multi-index
990   (apply values
991          (flip-multi-index multi-index))))</enscript>
992and the result inherits the safety and mutability of {{array}}.
993
994Otherwise, if {{array}} is mutable, then {{array-reverse}} returns
995
996<enscript highlight="scheme">(make-array
997 domain
998 (lambda multi-index
999   (apply (array-getter array)
1000          (flip-multi-index multi-index)))
1001   (lambda (v . multi-index)
1002     (apply (array-setter array)
1003            v
1004            (flip-multi-index multi-index)))))</enscript>
1005Finally, if {{array}} is immutable, then {{array-reverse}} returns
1006
1007<enscript highlight="scheme">(make-array
1008 domain
1009 (lambda multi-index
1010   (apply (array-getter array)
1011          (flip-multi-index multi-index))))) </enscript>
1012It is an error if {{array}} and {{flip?}} don't satisfy these requirements.
1013
1014The following example using {{array-reverse}} was motivated by [[https://funcall.blogspot.com/2020/01/palindromes-redux-and-sufficiently.html|a blog post by Joe Marshall]].
1015
1016<enscript highlight="scheme">(define (palindrome? s)
1017  (let ((n (string-length s)))
1018    (or (< n 2)
1019        (let* ((a
1020                ;; an array accessing the characters of s
1021                (make-array (make-interval (vector n))
1022                            (lambda (i)
1023                              (string-ref s i))))
1024               (ra
1025                ;; the array accessed in reverse order
1026                (array-reverse a))
1027               (half-domain
1028                (make-interval (vector (quotient n 2)))))
1029          (array-every
1030           char=?
1031           ;; the first half of s
1032           (array-extract a half-domain)
1033           ;; the reversed second half of s
1034           (array-extract ra half-domain))))))
1035
1036(palindrome? "") => #t
1037(palindrome? "a") => #t
1038(palindrome? "aa") => #t
1039(palindrome? "ab") => #f
1040(palindrome? "aba") => #t
1041(palindrome? "abc") => #f
1042(palindrome? "abba") => #t
1043(palindrome? "abca") => #f
1044(palindrome? "abbc") => #f</enscript>
1045<procedure>(array-sample array scales)</procedure>
1046
1047We assume that {{array}} is an array all of whose lower bounds are zero, and {{scales}} is a vector of positive exact integers whose length is the same as the dimension of {{array}}.
1048
1049Informally, if we construct a new matrix $S$ with the entries of {{scales}} on the main diagonal, then the $\vec i$th element of {{(array-sample array scales)}} is the $S\vec i$th element of {{array}}.
1050
1051More formally, if {{array}} is specialized, then {{array-sample}} returns
1052
1053<enscript highlight="scheme">(specialized-array-share
1054 array
1055 (interval-scale (array-domain array)
1056                 scales)
1057 (lambda multi-index
1058   (apply values
1059          (map * multi-index (vector->list scales)))))</enscript>
1060with the result inheriting the safety and mutability of {{array}}.
1061
1062Otherwise, if {{array}} is mutable, then {{array-sample}} returns
1063
1064<enscript highlight="scheme">(make-array
1065 (interval-scale (array-domain array)
1066                 scales)
1067 (lambda multi-index
1068   (apply (array-getter array)
1069          (map * multi-index (vector->list scales))))
1070 (lambda (v . multi-index)
1071   (apply (array-setter array)
1072          v
1073          (map * multi-index (vector->list scales)))))</enscript>
1074Finally, if {{array}} is immutable, then {{array-sample}} returns
1075
1076<enscript highlight="scheme">(make-array
1077 (interval-scale (array-domain array)
1078                 scales)
1079 (lambda multi-index
1080   (apply (array-getter array)
1081          (map * multi-index (vector->list scales)))))</enscript>
1082It is an error if {{array}} and {{scales}} don't satisfy these requirements.
1083
1084<procedure>(array-outer-product op array1 array2)</procedure>
1085
1086Implements the outer product of {{array1}} and {{array2}} with the operator {{op}}, similar to the APL function with the same name.
1087
1088Assume that {{array1}} and {{array2}} are arrays and that {{op}} is a function of two arguments. Assume that {{(list-tail l n)}} returns the list remaining after the first {{n}} items of the list {{l}} have been removed, and {{(list-take l n)}} returns a new list consisting of the first {{n}} items of the list {{l}}. Then {{array-outer-product}} returns the immutable array
1089
1090<enscript highlight="scheme">(make-array (interval-cartesian-product (array-domain array1)
1091                                        (array-domain array2))
1092            (lambda args
1093              (op (apply (array-getter array1) (list-take args (array-dimension array1)))
1094                  (apply (array-getter array2) (list-tail args (array-dimension array1))))))</enscript>
1095This operation can be considered a partial inverse to {{array-curry}}. It is an error if the arguments do not satisfy these assumptions.
1096
1097'''Note:''' You can see from the above definition that if {{C}} is {{(array-outer-product op A B)}}, then each call to {{(array-getter C)}} will call {{op}} as well as {{(array-getter A)}} and {{(array-getter B)}}. This implies that if all elements of {{C}} are eventually accessed, then {{(array-getter A)}} will be called {{(array-volume B)}} times; similarly {{(array-getter B)}} will be called {{(array-volume A)}} times.
1098
1099This implies that if {{(array-getter A)}} is expensive to compute (for example, if it's returning an array, as does {{array-curry}}) then the elements of {{A}} should be pre-computed if necessary and stored in a specialized array, typically using {{array-copy}}, before that specialized array is passed as an argument to {{array-outer-product}}. In the examples below, the code for Gaussian elimination applies {{array-outer-product}} to shared specialized arrays, which are of course themselves specialized arrays; the code for matrix multiplication and {{inner-product}} applies {{array-outer-product}} to curried arrays, so we apply {{array-copy}} to the arguments before passage to {{array-outer-product}}.
1100
1101<procedure>(array-map f array . arrays)</procedure>
1102
1103If {{array}}, {{(car arrays)}}, ... all have the same domain and {{f}} is a procedure, then {{array-map}} returns a new immutable array with the same domain and getter
1104
1105<enscript highlight="scheme">(lambda multi-index
1106  (apply f
1107         (map (lambda (g)
1108                (apply g multi-index))
1109              (map array-getter
1110                   (cons array arrays)))))</enscript>
1111It is assumed that {{f}} is appropriately defined to be evaluated in this context.
1112
1113It is expected that {{array-map}} and {{array-for-each}} will specialize the construction of
1114
1115<enscript highlight="scheme">(lambda multi-index
1116  (apply f
1117         (map (lambda (g)
1118                (apply g multi-index))
1119              (map array-getter
1120                   (cons array
1121                         arrays)))))</enscript>
1122It is an error to call {{array-map}} if its arguments do not satisfy these conditions.
1123
1124'''Note:''' The ease of constructing temporary arrays without allocating storage makes it trivial to imitate, e.g., Javascript's map with index. For example, we can add the index to each element of an array {{a}} by
1125
1126<enscript highlight="scheme">(array-map +
1127           a
1128           (make-array (array-domain a)
1129                       (lambda (i) i)))</enscript>
1130or even
1131
1132<enscript highlight="scheme">(make-array (array-domain a)
1133            (let ((a_ (array-getter a)))
1134              (lambda (i)
1135                (+ (a_ i) i))))</enscript>
1136<procedure>(array-for-each f array . arrays)</procedure>
1137
1138If {{array}}, {{(car arrays)}}, ... all have the same domain and {{f}} is an appropriate procedure, then {{array-for-each}} calls
1139
1140<enscript highlight="scheme">(interval-for-each
1141 (lambda multi-index
1142   (apply f
1143          (map (lambda (g)
1144                 (apply g multi-index))
1145               (map array-getter
1146                    (cons array
1147                          arrays)))))
1148 (array-domain array))</enscript>
1149In particular, {{array-for-each}} always walks the indices of the arrays in lexicographical order.
1150
1151It is expected that {{array-map}} and {{array-for-each}} will specialize the construction of
1152
1153<enscript highlight="scheme">(lambda multi-index
1154  (apply f
1155         (map (lambda (g)
1156                (apply g multi-index))
1157              (map array-getter
1158                   (cons array
1159                         arrays)))))</enscript>
1160It is an error to call {{array-for-each}} if its arguments do not satisfy these conditions.
1161
1162<procedure>(array-fold kons knil array)</procedure>
1163
1164If we use the defining relations for fold over lists from [[https://srfi.schemers.org/srfi-1/|SRFI 1]]:
1165
1166<enscript highlight="scheme">(fold kons knil lis)
1167    = (fold kons (kons (car lis) knil) (cdr lis))
1168(fold kons knil '())
1169    = knil
1170 </enscript>
1171then {{(array-fold kons knil array)}} returns
1172
1173 (fold kons knil (array->list array))
1174It is an error if {{array}} is not an array, or if {{kons}} is not a procedure.
1175
1176<procedure>(array-fold-right kons knil array)</procedure>
1177
1178If we use the defining relations for fold-right over lists from [[https://srfi.schemers.org/srfi-1/|SRFI 1]]:
1179
1180<enscript highlight="scheme">(fold-right kons knil lis)
1181    = (kons (car lis) (fold-right kons knil (cdr lis)))
1182(fold-right kons knil '())
1183    = knil</enscript>
1184then {{(array-fold-right kons knil array)}} returns
1185
1186 (fold-right kons knil (array->list array))
1187It is an error if {{array}} is not an array, or if {{kons}} is not a procedure.
1188
1189<procedure>(array-reduce op A)</procedure>
1190
1191We assume that {{A}} is an array and {{op}} is a procedure of two arguments that is associative, i.e., {{(op (op x y) z)}} is the same as {{(op x (op  y z))}}.
1192
1193Then {{(array-reduce op A)}} returns
1194
1195<enscript highlight="scheme">(let ((box '())
1196      (A_ (array-getter A)))
1197  (interval-for-each
1198   (lambda args
1199     (if (null? box)
1200         (set! box (list (apply A_ args)))
1201         (set-car! box (op (car box)
1202                           (apply A_ args)))))
1203   (array-domain A))
1204  (car box))</enscript>
1205The implementation is allowed to use the associativity of {{op}} to reorder the computations in {{array-reduce}}. It is an error if the arguments do not satisfy these conditions.
1206
1207As an example, we consider the finite sum: $$ S_m=\sum_{k=1}^m \frac 1{k^2}. $$ One can show that $$ \frac 1 {m+1}<\frac{\pi^2}6-S_m<\frac 1m. $$ We attempt to compute this in floating-point arithmetic in two ways. In the first, we apply {{array-reduce}} to an array containing the terms of the series, basically a serial computation. In the second, we divide the series into blocks of no more than 1,000 consecutive terms, apply {{array-reduce}} to get a new sequence of terms, and repeat the process. The second way is approximately what might happen with GPU computing.
1208
1209We find with $m=1{,}000{,}000{,}000$:
1210
1211<enscript highlight="scheme">(define A (make-array (make-interval '#(1) '#(1000000001))
1212                      (lambda (k)
1213                        (fl/ (flsquare (inexact k))))))
1214(define (block-sum A)
1215  (let ((N (interval-volume (array-domain A))))
1216    (cond ((<= N 1000)
1217           (array-reduce fl+ A))
1218          ((<= N (square 1000))
1219           (block-sum (array-map block-sum
1220                                 (array-tile A (vector (integer-sqrt N))))))
1221          (else
1222           (block-sum (array-map block-sum
1223                                 (array-tile A (vector (quotient N 1000)))))))))
1224(array-reduce fl+ A) => 1.644934057834575
1225(block-sum A)        => 1.6449340658482325</enscript>
1226Since $\pi^2/6\approx{}${{1.6449340668482264}}, we see using the first method that the difference $\pi^2/6-{}${{1.644934057834575}}${}\approx{}${{9.013651380840315e-9}} and with the second we have $\pi^2/6-{}${{1.6449340658482325}}${}\approx{}${{9.99993865491433e-10}}. The true difference should be between $\frac 1{1{,}000{,}000{,}001}\approx{}${{9.99999999e-10}} and $\frac 1{1{,}000{,}000{,}000}={}${{1e-9}}. The difference for the first method is about 10 times too big, and, in fact, will not change further because any further terms, when added to the partial sum, are too small to increase the sum after rounding-to-nearest in double-precision IEEE-754 floating-point arithmetic.
1227
1228<procedure>(array-any pred array1 array2 ...)</procedure>
1229
1230Assumes that {{array1}}, {{array2}}, etc., are arrays, all with the same domain, which we'll call {{interval}}. Also assumes that {{pred}} is a procedure that takes as many arguments as there are arrays and returns a single value.
1231
1232{{array-any}} first applies {{(array-getter array1)}}, etc., to the first element of {{interval}} in lexicographical order, to which values it then applies {{pred}}.
1233
1234If the result of {{pred}} is not {{#f}}, then that result is returned by {{array-any}}. If the result of {{pred}} is {{#f}}, then {{array-any}} continues with the second element of {{interval}}, etc., returning the first nonfalse value of {{pred}}.
1235
1236If {{pred}} always returns {{#f}}, then {{array-any}} returns {{#f}}.
1237
1238If it happens that {{pred}} is applied to the results of applying {{(array-getter array1)}}, etc., to the last element of {{interval}}, then this last call to {{pred}} is in tail position.
1239
1240The functions {{(array-getter array1)}}, etc., are applied only to those values of {{interval}} necessary to determine the result of {{array-any}}.
1241
1242It is an error if the arguments do not satisfy these assumptions.
1243
1244<procedure>(array-every pred array1 array2 ...)</procedure>
1245
1246Assumes that {{array1}}, {{array2}}, etc., are arrays, all with the same domain, which we'll call {{interval}}. Also assumes that {{pred}} is a procedure that takes as many arguments as there are arrays and returns a single value.
1247
1248{{array-every}} first applies {{(array-getter array1)}}, etc., to the first element of {{interval}} in lexicographical order, to which values it then applies {{pred}}.
1249
1250If the result of {{pred}} is {{#f}}, then that result is returned by {{array-every}}. If the result of {{pred}} is nonfalse, then {{array-every}} continues with the second element of {{interval}}, etc., returning the first value of {{pred}} that is {{#f}}.
1251
1252If {{pred}} always returns a nonfalse value, then the last nonfalse value returned by {{pred}} is also returned by {{array-every}}.
1253
1254If it happens that {{pred}} is applied to the results of applying {{(array-getter array1)}}, etc., to the last element of {{interval}}, then this last call to {{pred}} is in tail position.
1255
1256The functions {{(array-getter array1)}}, etc., are applied only to those values of {{interval}} necessary to determine the result of {{array-every}}.
1257
1258It is an error if the arguments do not satisfy these assumptions.
1259
1260<procedure>(array->list array)</procedure>
1261
1262Stores the elements of {{array}} into a newly allocated list in lexicographical order. It is an error if {{array}} is not an array.
1263
1264It is guaranteed that {{(array-getter array)}} is called precisely once for each multi-index in {{(array-domain array)}} in lexicographical order.
1265
1266<procedure>(list->array l domain [ result-storage-class generic-storage-class ] [ mutable? (specialized-array-default-mutable?) ] [ safe? (specialized-array-default-safe?) ])</procedure>
1267
1268Assumes that {{l}} is an list, {{domain}} is an interval with volume the same as the length of {{l}}, {{result-storage-class}} is a storage class that can manipulate all the elements of {{l}}, and {{mutable?}} and {{safe?}} are booleans.
1269
1270Returns a specialized array with domain {{domain}} whose elements are the elements of the list {{l}} stored in lexicographical order. The result is mutable or safe depending on the values of {{mutable?}} and {{safe?}}.
1271
1272It is an error if the arguments do not satisfy these assumptions, or if any element of {{l}} cannot be stored in the body of {{result-storage-class}}, and this last error shall be detected and raised.
1273
1274<procedure>(array-assign! destination source)</procedure>
1275
1276Assumes that {{destination}} is a mutable array and {{source}} is an array, and that the elements of {{source}} can be stored into {{destination}}.
1277
1278The array {{destination}} must be compatible with {{source}}, in the sense that either {{destination}} and {{source}} have the same domain, or {{destination}} is a specialized array whose elements are stored adjacently and in order in its body and whose domain has the same volume as the domain of {{source}}.
1279
1280Evaluates {{(array-getter source)}} on the multi-indices in {{(array-domain source)}} in lexicographical order, and assigns each value to the multi-index in {{destination}} in the same lexicographical order.
1281
1282It is an error if the arguments don't satisfy these assumptions.
1283
1284If assigning any element of {{destination}} affects the value of any element of {{source}}, then the result is undefined.
1285
1286'''Note:''' If the domains of {{destination}} and {{source}} are not the same, one can think of {{destination}} as a ''reshaped'' copy of {{source}}.
1287
1288<procedure>(array-ref A i0 . i-tail)</procedure>
1289
1290Assumes that {{A}} is an array, and every element of {{(cons i0 i-tail)}} is an exact integer.
1291
1292Returns {{(apply (array-getter A) i0 i-tail)}}.
1293
1294It is an error if {{A}} is not an array, or if the number of arguments specified is not the correct number for {{(array-getter A)}}.
1295
1296<procedure>(array-set! A v i0 . i-tail)</procedure>
1297
1298Assumes that {{A}} is a mutable array, that {{v}} is a value that can be stored within that array, and that every element of {{(cons i0 i-tail)}} is an exact integer.
1299
1300Returns {{(apply (array-setter A) v i0 i-tail)}}.
1301
1302It is an error if {{A}} is not a mutable array, if {{v}} is not an appropriate value to be stored in that array, or if the number of arguments specified is not the correct number for {{(array-setter A)}}.
1303
1304'''Note:''' In the sample implementation, because {{array-ref}} and {{array-set!}} take a variable number of arguments and they must check that {{A}} is an array of the appropriate type, programs written in a style using these functions, rather than the style in which {{1D-Haar-loop}} is coded below, can take up to three times as long runtime.
1305
1306'''Note:''' In the sample implementation, checking whether the multi-indices are exact integers and within the domain of the array, and checking whether the value is appropriate for storage into the array, is delegated to the underlying definition of the array argument. If the first argument is a safe specialized array, then these items are checked; if it is an unsafe specialized array, they are not. If it is a generalized array, it is up to the programmer whether to define the getter and setter of the array to check the correctness of the arguments.
1307
1308<procedure>(specialized-array-reshape array new-domain [ copy-on-failure? #f ])</procedure>
1309
1310Assumes that {{array}} is a specialized array, {{new-domain}} is an interval with the same volume as {{(array-domain array)}}, and {{copy-on-failure?}}, if given, is a boolean.
1311
1312If there is an affine map that takes the multi-indices in {{new-domain}} to the cells in {{(array-body array)}} storing the elements of {{array}} in lexicographical order, returns a new specialized array, with the same body and elements as {{array}} and domain {{new-domain}}. The result inherits its mutability and safety from {{array}}.
1313
1314If there is not an affine map that takes the multi-indices in {{new-domain}} to the cells storing the elements of {{array}} in lexicographical order and {{copy-on-failure?}} is {{#t}}, then returns a specialized array copy of {{array}} with domain {{new-domain}}, storage class {{(array-storage-class array)}}, mutability {{(mutable-array? array)}}, and safety {{(array-safe? array)}}.
1315
1316It is an error if these conditions on the arguments are not met.
1317
1318'''Note:''' The code in the sample implementation to determine whether there exists an affine map from {{new-domain}} to the multi-indices of the elements of {{array}} in lexicographical order is modeled on the corresponding code in the Python library NumPy.
1319
1320'''Note:''' In the sample implementation, if an array cannot be reshaped and {{copy-on-failure?}} is {{#f}}, an error is raised in tail position. An implementation might want to replace this error call with a continuable exception to give the programmer more flexibility.
1321
1322'''Examples:''' Reshaping an array is not a Bawden-type array transform. For example, we use {{array-display}} defined below to see:
1323
1324<enscript highlight="scheme">;;; The entries of A are the multi-indices of the locations
1325
1326(define A (array-copy (make-array (make-interval '#(3 4)) list)))
1327
1328(array-display A)
1329
1330;;; Displays
1331
1332;;; (0 0)   (0 1)   (0 2)   (0 3)
1333;;; (1 0)   (1 1)   (1 2)   (1 3)
1334;;; (2 0)   (2 1)   (2 2)   (2 3)
1335
1336(array-display (array-rotate A 1))
1337
1338;;; Displays
1339
1340;;; (0 0)   (1 0)   (2 0)
1341;;; (0 1)   (1 1)   (2 1)
1342;;; (0 2)   (1 2)   (2 2)
1343;;; (0 3)   (1 3)   (2 3)
1344
1345(array-display (specialized-array-reshape A (make-interval '#(4 3))))
1346
1347;;; Displays
1348
1349;;; (0 0)   (0 1)   (0 2)
1350;;; (0 3)   (1 0)   (1 1)
1351;;; (1 2)   (1 3)   (2 0)
1352;;; (2 1)   (2 2)   (2 3)
1353
1354(define B (array-sample A '#(2 1)))
1355
1356(array-display B)
1357
1358;;; Displays
1359
1360;;; (0 0)   (0 1)   (0 2)   (0 3)
1361;;; (2 0)   (2 1)   (2 2)   (2 3)
1362
1363(array-display (specialized-array-reshape B (make-interval '#(8)))) => fails
1364
1365(array-display (specialized-array-reshape B (make-interval '#(8)) #t))
1366
1367;;; Displays
1368
1369;;; (0 0)   (0 1)   (0 2)   (0 3)   (2 0)   (2 1)   (2 2)   (2 3)</enscript>
1370The following examples succeed:
1371
1372<enscript highlight="scheme">(specialized-array-reshape
1373 (array-copy (make-array (make-interval '#(2 1 3 1)) list))
1374 (make-interval '#(6)))
1375(specialized-array-reshape
1376 (array-copy (make-array (make-interval '#(2 1 3 1)) list))
1377 (make-interval '#(3 2)))
1378(specialized-array-reshape
1379 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
1380 (make-interval '#(6)))
1381(specialized-array-reshape
1382 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
1383 (make-interval '#(3 2)))
1384(specialized-array-reshape
1385 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
1386 (make-interval '#(3 2)))
1387(specialized-array-reshape
1388 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
1389 (make-interval '#(3 1 2 1)))
1390(specialized-array-reshape
1391 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#f #f #f #t)) '#(1 1 2 1))
1392 (make-interval '#(4)))
1393(specialized-array-reshape
1394 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#t #f #t #t)) '#(1 1 2 1))
1395 (make-interval '#(4)))</enscript>
1396The following examples raise an exception:
1397
1398<enscript highlight="scheme">(specialized-array-reshape
1399 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
1400 (make-interval '#(6)))
1401(specialized-array-reshape
1402 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
1403 (make-interval '#(3 2)))
1404(specialized-array-reshape
1405 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #f))
1406 (make-interval '#(6)))
1407(specialized-array-reshape
1408 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #t))
1409 (make-interval '#(3 2)))
1410(specialized-array-reshape
1411 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t)) '#(1 1 2 1))
1412 (make-interval '#(4)) )
1413(specialized-array-reshape
1414 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#f #f #t #t)) '#(1 1 2 1))
1415 (make-interval '#(4)))</enscript>
1416In the next examples, we start with vector fields, $100\times 100$ arrays of 4-vectors. In one example, we reshape each large array to $100\times 100\times2\times2$ vector fields (so we consider each 4-vector as a $2\times 2$ matrix), and multiply the $2\times 2$ matrices together. In the second example, we reshape each 4-vector to a $2\times 2$ matrix individually, and compare the times.
1417
1418<enscript highlight="scheme">(define interval-flat (make-interval '#(100 100 4)))
1419
1420(define interval-2x2  (make-interval '#(100 100 2 2)))
1421
1422(define A (array-copy (make-array interval-flat (lambda args (random-integer 5)))))
1423
1424(define B (array-copy (make-array interval-flat (lambda args (random-integer 5)))))
1425
1426(define C (array-copy (make-array interval-flat (lambda args 0))))
1427
1428(define (2x2-matrix-multiply-into! A B C)
1429  (let ((C! (array-setter C))
1430        (A_ (array-getter A))
1431        (B_ (array-getter B)))
1432    (C! (+ (* (A_ 0 0) (B_ 0 0))
1433           (* (A_ 0 1) (B_ 1 0)))
1434        0 0)
1435    (C! (+ (* (A_ 0 0) (B_ 0 1))
1436           (* (A_ 0 1) (B_ 1 1)))
1437        0 1)
1438    (C! (+ (* (A_ 1 0) (B_ 0 0))
1439           (* (A_ 1 1) (B_ 1 0)))
1440        1 0)
1441    (C! (+ (* (A_ 1 0) (B_ 0 1))
1442           (* (A_ 1 1) (B_ 1 1)))
1443        1 1)))
1444
1445;;; Reshape A, B, and C to change all the 4-vectors to 2x2 matrices
1446
1447(time
1448 (array-for-each 2x2-matrix-multiply-into!
1449                 (array-curry (specialized-array-reshape A interval-2x2) 2)
1450                 (array-curry (specialized-array-reshape B interval-2x2) 2)
1451                 (array-curry (specialized-array-reshape C interval-2x2) 2)))
1452;;; Displays
1453
1454;;;    0.015186 secs real time
1455;;;    0.015186 secs cpu time (0.015186 user, 0.000000 system)
1456;;;    4 collections accounting for 0.004735 secs real time (0.004732 user, 0.000000 system)
1457;;;    46089024 bytes allocated
1458;;;    no minor faults
1459;;;    no major faults
1460
1461;;; Reshape each 4-vector to a 2x2 matrix individually
1462
1463(time
1464 (array-for-each (lambda (A B C)
1465                   (2x2-matrix-multiply-into!
1466                    (specialized-array-reshape A (make-interval '#(2 2)))
1467                    (specialized-array-reshape B (make-interval '#(2 2)))
1468                    (specialized-array-reshape C (make-interval '#(2 2)))))
1469                 (array-curry A 1)
1470                 (array-curry B 1)
1471                 (array-curry C 1)))
1472
1473;;; Displays
1474
1475;;;    0.039193 secs real time
1476;;;    0.039193 secs cpu time (0.039191 user, 0.000002 system)
1477;;;    6 collections accounting for 0.006855 secs real time (0.006851 user, 0.000001 system)
1478;;;    71043024 bytes allocated
1479;;;    no minor faults
1480;;;    no major faults</enscript>
1481=== Implementation
1482
1483We provide an implementation in [[https://github.com/gambit/gambit|Gambit Scheme]]; the nonstandard techniques used in the implementation are: DSSSL-style optional and keyword arguments; a unique object to indicate absent arguments; {{define-structure}}; and {{define-macro}}.
1484
1485There is a [[https://github.com/scheme-requests-for-implementation/srfi-179|git repository]] of this document, a sample implementation, a test file, and other materials.
1486
1487=== Relationship to other SRFIs
1488
1489Final SRFIs [[#SRFI-25|25]], [[#SRFI-47|47]], [[#SRFI-58|58]], and [[#SRFI-63|63]] deal with "Multi-dimensional Array Primitives", "Array", "Array Notation", and "Homogeneous and Heterogeneous Arrays", respectively. Each of these previous SRFIs deal with what we call in this SRFI specialized arrays. Many of the functions in these previous SRFIs have corresponding forms in this SRFI. For example, from [[https://srfi.schemers.org/srfi-63/|SRFI 63]], we can translate:
1490
1491; {{(array? obj)}}
1492: {{(array? obj)}}
1493; {{(array-rank A)}}
1494: {{(array-dimension A)}}
1495; {{(make-array prototype k1 ...)}}
1496: {{(make-specialized-array (make-interval (vector k1 ...)) storage-class)}}.
1497; {{(make-shared-array A mapper k1 ...)}}
1498: {{(specialized-array-share A (make-interval (vector k1 ...)) mapper)}}
1499; {{(array-in-bounds? A index1 ...)}}
1500: {{(interval-contains-multi-index? (array-domain A) index1 ...)}}
1501; {{(array-ref A k1 ...)}}
1502: {{(let ((A_ (array-getter A))) ... (A_ k1 ...) ... )}} or {{(array-ref A k1 ...)}}
1503; {{(array-set! A obj k1 ...)}}
1504: {{(let ((A! (array-setter A))) ... (A! obj k1 ...) ...)}} or {{(array-set! A obj k1 ...)}}
1505
1506At the same time, this SRFI has some special features:
1507
1508* Intervals, used as the domains of arrays in this SRFI, are useful objects in their own rights, with their own procedures. We make a sharp distinction between the domains of arrays and the arrays themselves.
1509* Intervals can have nonzero lower bounds in each dimension.
1510* Intervals cannot be empty.
1511* Arrays must have a getter, but may have no setter.
1512* There are many predefined array transformations: {{array-extract}}, {{array-tile}}, {{array-translate}}, {{array-permute}}, {{array-rotate}}, {{array-sample}}, {{array-reverse}}.
1513
1514=== Other examples
1515
1516Image processing applications provided significant motivation for this SRFI.
1517
1518'''Manipulating images in PGM format.''' On a system with eight-bit chars, one can write routines to read and write greyscale images in the PGM format of the netpbm package as follows. The lexicographical order in {{array-copy}} guarantees the the correct order of execution of the input procedures:
1519
1520<enscript highlight="scheme">(define make-pgm   cons)
1521(define pgm-greys  car)
1522(define pgm-pixels cdr)
1523
1524(define (read-pgm file)
1525
1526  (define (read-pgm-object port)
1527    (skip-white-space port)
1528    (let ((o (read port)))
1529      ;; to skip the newline or next whitespace
1530      (read-char port)
1531      (if (eof-object? o)
1532          (error "reached end of pgm file")
1533          o)))
1534
1535  (define (skip-to-end-of-line port)
1536    (let loop ((ch (read-char port)))
1537      (if (not (eq? ch #\newline))
1538          (loop (read-char port)))))
1539
1540  (define (white-space? ch)
1541    (case ch
1542      ((#\newline #\space #\tab) #t)
1543      (else #f)))
1544
1545  (define (skip-white-space port)
1546    (let ((ch (peek-char port)))
1547      (cond ((white-space? ch)
1548             (read-char port)
1549             (skip-white-space port))
1550            ((eq? ch #\#)
1551             (skip-to-end-of-line port)
1552             (skip-white-space port))
1553            (else #f))))
1554
1555  ;; The image file formats defined in netpbm
1556  ;; are problematical because they read the data
1557  ;; in the header as variable-length ISO-8859-1 text,
1558  ;; including arbitrary whitespace and comments,
1559  ;; and then they may read the rest of the file
1560  ;; as binary data.
1561  ;; So we give here a solution of how to deal
1562  ;; with these subtleties in Gambit Scheme.
1563
1564  (call-with-input-file
1565      (list path:          file
1566            char-encoding: 'ISO-8859-1
1567            eol-encoding:  'lf)
1568    (lambda (port)
1569
1570      ;; We're going to read text for a while,
1571      ;; then switch to binary.
1572      ;; So we need to turn off buffering until
1573      ;; we switch to binary.
1574
1575      (port-settings-set! port '(buffering: #f))
1576
1577      (let* ((header (read-pgm-object port))
1578             (columns (read-pgm-object port))
1579             (rows (read-pgm-object port))
1580             (greys (read-pgm-object port)))
1581
1582        ;; Now we switch back to buffering
1583        ;; to speed things up.
1584
1585        (port-settings-set! port '(buffering: #t))
1586
1587        (make-pgm
1588         greys
1589         (array-copy
1590          (make-array
1591           (make-interval (vector rows columns))
1592           (cond ((or (eq? header 'p5)
1593                      (eq? header 'P5))
1594                  ;; pgm binary
1595                  (if (< greys 256)
1596                      ;; one byte/pixel
1597                      (lambda (i j)
1598                        (char->integer
1599                         (read-char port)))
1600                      ;; two bytes/pixel,
1601                      ;;little-endian
1602                      (lambda (i j)
1603                        (let* ((first-byte
1604                                (char->integer
1605                                 (read-char port)))
1606                               (second-byte
1607                                (char->integer
1608                                 (read-char port))))
1609                          (+ (* second-byte 256)
1610                             first-byte)))))
1611                 ;; pgm ascii
1612                 ((or (eq? header 'p2)
1613                      (eq? header 'P2))
1614                  (lambda (i j)
1615                      (read port)))
1616                   (else
1617                    (error "not a pgm file"))))
1618          (if (< greys 256)
1619              u8-storage-class
1620              u16-storage-class)))))))
1621
1622(define (write-pgm pgm-data file #!optional force-ascii)
1623  (call-with-output-file
1624      (list path:          file
1625            char-encoding: 'ISO-8859-1
1626            eol-encoding:  'lf)
1627    (lambda (port)
1628      (let* ((greys
1629              (pgm-greys pgm-data))
1630             (pgm-array
1631              (pgm-pixels pgm-data))
1632             (domain
1633              (array-domain pgm-array))
1634             (rows
1635              (fx- (interval-upper-bound domain 0)
1636                   (interval-lower-bound domain 0)))
1637             (columns
1638              (fx- (interval-upper-bound domain 1)
1639                   (interval-lower-bound domain 1))))
1640        (if force-ascii
1641            (display "P2" port)
1642            (display "P5" port))
1643        (newline port)
1644        (display columns port) (display  port)
1645        (display rows port) (newline port)
1646        (display greys port) (newline port)
1647        (array-for-each
1648         (if force-ascii
1649             (let ((next-pixel-in-line 1))
1650               (lambda (p)
1651                 (write p port)
1652                 (if (fxzero? (fxand next-pixel-in-line 15))
1653                     (begin
1654                       (newline port)
1655                       (set! next-pixel-in-line 1))
1656                     (begin
1657                       (display  port)
1658                       (set! next-pixel-in-line
1659                             (fx+ 1 next-pixel-in-line))))))
1660             (if (fx< greys 256)
1661                 (lambda (p)
1662                   (write-u8 p port))
1663                 (lambda (p)
1664                   (write-u8 (fxand p 255) port)
1665                   (write-u8 (fxarithmetic-shift-right p 8)
1666                             port))))
1667         pgm-array)))))</enscript>
1668One can write a a routine to convolve an image with a filter as follows:
1669
1670<enscript highlight="scheme">(define (array-convolve source filter)
1671  (let* ((source-domain
1672          (array-domain source))
1673         (S_
1674          (array-getter source))
1675         (filter-domain
1676          (array-domain filter))
1677         (F_
1678          (array-getter filter))
1679         (result-domain
1680          (interval-dilate
1681           source-domain
1682           ;; the left bound of an interval is an equality,
1683           ;; the right bound is an inequality, hence the
1684           ;; the difference in the following two expressions
1685           (vector-map -
1686                       (interval-lower-bounds->vector filter-domain))
1687           (vector-map (lambda (x)
1688                         (- 1 x))
1689                       (interval-upper-bounds->vector filter-domain)))))
1690    (make-array result-domain
1691                (lambda (i j)
1692                  (array-fold
1693                   (lambda (p q)
1694                     (+ p q))
1695                   0
1696                   (make-array
1697                    filter-domain
1698                    (lambda (k l)
1699                      (* (S_ (+ i k)
1700                             (+ j l))
1701                         (F_ k l))))))
1702                )))</enscript>
1703together with some filters
1704
1705<enscript highlight="scheme">(define sharpen-filter
1706  (list->array
1707   '(0 -1  0
1708    -1  5 -1
1709     0 -1  0)
1710   (make-interval '#(-1 -1) '#(2 2))))
1711
1712(define edge-filter
1713  (list->array
1714   '(0 -1  0
1715    -1  4 -1
1716     0 -1  0)
1717   (make-interval '#(-1 -1) '#(2 2))))</enscript>
1718Our computations might results in pixel values outside the valid range, so we define
1719
1720<enscript highlight="scheme">(define (round-and-clip pixel max-grey)
1721  (max 0 (min (exact (round pixel)) max-grey)))</enscript>
1722We can then compute edges and sharpen a test image as follows:
1723
1724<enscript highlight="scheme">(define test-pgm (read-pgm "girl.pgm"))
1725
1726(let ((greys (pgm-greys test-pgm)))
1727  (write-pgm
1728   (make-pgm
1729    greys
1730    (array-map (lambda (p)
1731                 (round-and-clip p greys))
1732               (array-convolve
1733                (pgm-pixels test-pgm)
1734                sharpen-filter)))
1735   "sharper-test.pgm"))
1736
1737(let* ((greys (pgm-greys test-pgm))
1738       (edge-array
1739        (array-copy
1740         (array-map
1741          abs
1742          (array-convolve
1743           (pgm-pixels test-pgm)
1744           edge-filter))))
1745       (max-pixel
1746        (array-fold max 0 edge-array))
1747       (normalizer
1748        (inexact (/ greys max-pixel))))
1749  (write-pgm
1750   (make-pgm
1751    greys
1752    (array-map (lambda (p)
1753                 (- greys
1754                    (round-and-clip (* p normalizer) greys)))
1755               edge-array))
1756   "edge-test.pgm"))</enscript>
1757'''Viewing two-dimensional slices of three-dimensional data.''' One example might be viewing two-dimensional slices of three-dimensional data in different ways. If one has a $1024 \times 512\times 512$ 3D image of the body stored as a variable {{body}}, then one could get 1024 axial views, each $512\times512$, of this 3D body by {{(array-curry body 2)}}; or 512 median views, each $1024\times512$, by {{(array-curry (array-permute body '#(1 0 2)) 2)}}; or finally 512 frontal views, each again $1024\times512$ pixels, by {{(array-curry (array-permute body '#(2 0 1)) 2)}}; see [[https://en.wikipedia.org/wiki/Anatomical_plane|Anatomical plane]]. Note that the first permutation is not a rotation—you want to have the head up in both the median and frontal views.
1758
1759'''Calculating second differences of images.''' For another example, if a real-valued function is defined on a two-dimensional interval $I$, its second difference in the direction $d$ at the point $x$ is defined as $\Delta^2_df(x)=f(x+2d)-2f(x+d)+f(x)$, and this function is defined only for those $x$ for which $x$, $x+d$, and $x+2d$ are all in $I$. See the beginning of the section on "Moduli of smoothness" in [[https://www.math.purdue.edu/~lucier/692/related_papers_summaries.html#Wavelets-and-approximation-theory|these notes on wavelets and approximation theory]] for more details.
1760
1761Using this definition, the following code computes all second-order forward differences of an image in the directions $d,2 d,3 d,\ldots$, defined only on the domains where this makes sense:
1762
1763<enscript highlight="scheme">(define (all-second-differences image direction)
1764  (let ((image-domain (array-domain image)))
1765    (let loop ((i 1)
1766               (result '()))
1767      (let ((negative-scaled-direction
1768             (vector-map (lambda (j)
1769                           (* -1 j i))
1770                         direction))
1771            (twice-negative-scaled-direction
1772             (vector-map (lambda (j)
1773                           (* -2 j i))
1774                         direction)))
1775        (cond ((interval-intersect
1776                image-domain
1777                (interval-translate
1778                 image-domain
1779                 negative-scaled-direction)
1780                (interval-translate
1781                 image-domain
1782                 twice-negative-scaled-direction))
1783               =>
1784               (lambda (subdomain)
1785                 (loop
1786                  (+ i 1)
1787                  (cons
1788                   (array-copy
1789                    (array-map
1790                     (lambda (f_i f_i+d f_i+2d)
1791                       (+ f_i+2d
1792                          (* -2. f_i+d)
1793                          f_i))
1794                     (array-extract
1795                      image
1796                      subdomain)
1797                     (array-extract
1798                      (array-translate
1799                       image
1800                       negative-scaled-direction)
1801                      subdomain)
1802                     (array-extract
1803                      (array-translate
1804                       image
1805                       twice-negative-scaled-direction)
1806                      subdomain)))
1807                   result))))
1808              (else
1809               (reverse result)))))))</enscript>
1810We can define a small synthetic image of size 8x8 pixels and compute its second differences in various directions:
1811
1812<enscript highlight="scheme">(define image
1813 (array-copy
1814  (make-array (make-interval '#(8 8))
1815              (lambda (i j)
1816                (exact->inexact (+ (* i i) (* j j)))))))
1817
1818(define (expose difference-images)
1819  (pretty-print (map (lambda (difference-image)
1820                       (list (array-domain difference-image)
1821                             (array->list difference-image)))
1822                     difference-images)))
1823
1824(begin
1825  (display
1826   "\nSecond-differences in the direction $k\times (1,0)$:\n")
1827  (expose (all-second-differences image '#(1 0)))
1828  (display
1829   "\nSecond-differences in the direction $k\times (1,1)$:\n")
1830  (expose (all-second-differences image '#(1 1)))
1831  (display
1832   "\nSecond-differences in the direction $k\times (1,-1)$:\n")
1833  (expose (all-second-differences image '#(1 -1))))</enscript>
1834On Gambit 4.8.5, this yields (after some hand editing):
1835
1836<enscript highlight="scheme">Second-differences in the direction $k\times (1,0)$:
1837((#<##interval #2 lower-bounds: #(0 0) upper-bounds: #(6 8)>
1838 (2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
1839  2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
1840  2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.))
1841 (#<##interval #3 lower-bounds: #(0 0) upper-bounds: #(4 8)>
1842  (8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8.
1843   8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8.))
1844 (#<##interval #4 lower-bounds: #(0 0) upper-bounds: #(2 8)>
1845  (18. 18. 18. 18. 18. 18. 18. 18. 18.
1846   18. 18. 18. 18. 18. 18. 18.)))
1847
1848Second-differences in the direction $k\times (1,1)$:
1849((#<##interval #5 lower-bounds: #(0 0) upper-bounds: #(6 6)>
1850  (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
1851   4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.))
1852 (#<##interval #6 lower-bounds: #(0 0) upper-bounds: #(4 4)>
1853  (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.
1854   16. 16.))
1855 (#<##interval #7 lower-bounds: #(0 0) upper-bounds: #(2 2)>
1856  (36. 36. 36. 36.)))
1857
1858Second-differences in the direction $k\times (1,-1)$:
1859((#<##interval #8 lower-bounds: #(0 2) upper-bounds: #(6 8)>
1860  (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
1861   4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.))
1862 (#<##interval #9 lower-bounds: #(0 4) upper-bounds: #(4 8)>
1863  (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.
1864   16. 16.))
1865 (#<##interval #10 lower-bounds: #(0 6) upper-bounds: #(2 8)>
1866  (36. 36. 36. 36.)))
1867</enscript>
1868You can see that with differences in the direction of only the first coordinate, the domains of the difference arrays get smaller in the first coordinate while staying the same in the second coordinate, and with differences in the diagonal directions, the domains of the difference arrays get smaller in both coordinates.
1869
1870'''Separable operators.''' Many multi-dimensional transforms in signal processing are ''separable'', in that the multi-dimensional transform can be computed by applying one-dimensional transforms in each of the coordinate directions. Examples of such transforms include the Fast Fourier Transform and the [[https://arxiv.org/abs/1210.1944|Fast Hyperbolic Wavelet Transform]]. Each one-dimensional subdomain of the complete domain is called a ''pencil'', and the same one-dimensional transform is applied to all pencils in a given direction. Given the one-dimensional array transform, one can define the multidimensional transform as follows:
1871
1872<enscript highlight="scheme">(define (make-separable-transform 1D-transform)
1873  (lambda (a)
1874    (let ((n (array-dimension a)))
1875      (do ((d 0 (fx+ d 1)))
1876          ((fx= d n))
1877        (array-for-each
1878         1D-transform
1879         (array-curry (array-rotate a d) 1))))))</enscript>
1880Here we have cycled through all rotations, putting each axis in turn at the end, and then applied {{1D-transform}} to each of the pencils along that axis.
1881
1882Wavelet transforms in particular are calculated by recursively applying a transform to an array and then downsampling the array; the inverse transform recursively downsamples and then applies a transform. So we define the following primitives:
1883
1884<enscript highlight="scheme">(define (recursively-apply-transform-and-downsample transform)
1885  (lambda (a)
1886    (let ((sample-vector (make-vector (array-dimension a) 2)))
1887      (define (helper a)
1888        (if (fx< 1 (interval-upper-bound (array-domain a) 0))
1889            (begin
1890              (transform a)
1891              (helper (array-sample a sample-vector)))))
1892      (helper a))))
1893
1894(define (recursively-downsample-and-apply-transform transform)
1895  (lambda (a)
1896    (let ((sample-vector (make-vector (array-dimension a) 2)))
1897      (define (helper a)
1898        (if (fx< 1 (interval-upper-bound (array-domain a) 0))
1899            (begin
1900              (helper (array-sample a sample-vector))
1901              (transform a))))
1902      (helper a))))</enscript>
1903By adding a single loop that calculates scaled sums and differences of adjacent elements in a one-dimensional array, we can define various <span id="Haar">Haar wavelet transforms</span> as follows:
1904
1905<enscript highlight="scheme">(define (1D-Haar-loop a)
1906  (let ((a_ (array-getter a))
1907        (a! (array-setter a))
1908        (n (interval-upper-bound (array-domain a) 0)))
1909    (do ((i 0 (fx+ i 2)))
1910        ((fx= i n))
1911      (let* ((a_i               (a_ i))
1912             (a_i+1             (a_ (fx+ i 1)))
1913             (scaled-sum        (fl/ (fl+ a_i a_i+1) (flsqrt 2.0)))
1914             (scaled-difference (fl/ (fl- a_i a_i+1) (flsqrt 2.0))))
1915        (a! scaled-sum i)
1916        (a! scaled-difference (fx+ i 1))))))
1917
1918(define 1D-Haar-transform
1919  (recursively-apply-transform-and-downsample 1D-Haar-loop))
1920
1921(define 1D-Haar-inverse-transform
1922  (recursively-downsample-and-apply-transform 1D-Haar-loop))
1923
1924(define hyperbolic-Haar-transform
1925  (make-separable-transform 1D-Haar-transform))
1926
1927(define hyperbolic-Haar-inverse-transform
1928  (make-separable-transform 1D-Haar-inverse-transform))
1929
1930(define Haar-transform
1931  (recursively-apply-transform-and-downsample
1932   (make-separable-transform 1D-Haar-loop)))
1933
1934(define Haar-inverse-transform
1935  (recursively-downsample-and-apply-transform
1936   (make-separable-transform 1D-Haar-loop)))</enscript>
1937We then define an image that is a multiple of a single, two-dimensional hyperbolic Haar wavelet, compute its hyperbolic Haar transform (which should have only one nonzero coefficient), and then the inverse transform:
1938
1939<enscript highlight="scheme">(let ((image
1940       (array-copy
1941        (make-array (make-interval '#(4 4))
1942                    (lambda (i j)
1943                      (case i
1944                        ((0) 1.)
1945                        ((1) -1.)
1946                        (else 0.)))))))
1947  (display "
1948Initial image:
1949")
1950  (pretty-print (list (array-domain image)
1951                      (array->list image)))
1952  (hyperbolic-Haar-transform image)
1953  (display "\nArray of hyperbolic Haar wavelet coefficients: \n")
1954  (pretty-print (list (array-domain image)
1955                      (array->list image)))
1956  (hyperbolic-Haar-inverse-transform image)
1957  (display "\nReconstructed image: \n")
1958  (pretty-print (list (array-domain image)
1959                      (array->list image))))</enscript>
1960This yields:
1961
1962<enscript highlight="scheme">Initial image:
1963(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
1964 (1. 1. 1. 1. -1. -1. -1. -1. 0. 0. 0. 0. 0. 0. 0. 0.))
1965
1966Array of hyperbolic Haar wavelet coefficients:
1967(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
1968 (0. 0. 0. 0. 2.8284271247461894 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.))
1969
1970Reconstructed image:
1971(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
1972 (.9999999999999996
1973  .9999999999999996
1974  .9999999999999996
1975  .9999999999999996
1976  -.9999999999999996
1977  -.9999999999999996
1978  -.9999999999999996
1979  -.9999999999999996
1980  0.
1981  0.
1982  0.
1983  0.
1984  0.
1985  0.
1986  0.
1987  0.))</enscript>
1988In perfect arithmetic, this hyperbolic Haar transform is ''orthonormal'', in that the sum of the squares of the elements of the image is the same as the sum of the squares of the hyperbolic Haar coefficients of the image. We can see that this is approximately true here.
1989
1990We can apply the (nonhyperbolic) Haar transform to the same image as follows:
1991
1992<enscript highlight="scheme"> (let ((image
1993       (array-copy
1994        (make-array (make-interval '#(4 4))
1995                    (lambda (i j)
1996                      (case i
1997                        ((0) 1.)
1998                        ((1) -1.)
1999                        (else 0.)))))))
2000  (display "\nInitial image:\n")
2001  (pretty-print (list (array-domain image)
2002                      (array->list image)))
2003  (Haar-transform image)
2004  (display "\nArray of Haar wavelet coefficients: \n")
2005  (pretty-print (list (array-domain image)
2006                      (array->list image)))
2007  (Haar-inverse-transform image)
2008  (display "\nReconstructed image: \n")
2009  (pretty-print (list (array-domain image)
2010                      (array->list image))))</enscript>
2011This yields:
2012
2013<enscript highlight="scheme">Initial image:
2014(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
2015 (1. 1. 1. 1. -1. -1. -1. -1. 0. 0. 0. 0. 0. 0. 0. 0.))
2016
2017Array of Haar wavelet coefficients:
2018(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
2019 (0. 0. 0. 0. 1.9999999999999998 0. 1.9999999999999998 0. 0. 0. 0. 0. 0. 0. 0. 0.))
2020
2021Reconstructed image:
2022(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
2023 (.9999999999999997
2024  .9999999999999997
2025  .9999999999999997
2026  .9999999999999997
2027  -.9999999999999997
2028  -.9999999999999997
2029  -.9999999999999997
2030  -.9999999999999997
2031  0.
2032  0.
2033  0.
2034  0.
2035  0.
2036  0.
2037  0.
2038  0.))</enscript>
2039You see in this example that this particular image has two, not one, nonzero coefficients in the two-dimensional Haar transform, which is again approximately orthonormal.
2040
2041'''Matrix multiplication and Gaussian elimination.''' While we have avoided conflating matrices and arrays, we give here some examples of matrix operations defined using operations from this SRFI.
2042
2043Given a nonsingular square matrix $A$ we can overwrite $A$ with lower-triangular matrix $L$ with ones on the diagonal and upper-triangular matrix $U$ so that $A=LU$ as follows. (We assume "pivoting" isn't needed.) For example, if $$A=\begin{pmatrix} a_{11}&a_{12}&a_{13}\\ a_{21}&a_{22}&a_{23}\\ a_{31}&a_{32}&a_{33}\end{pmatrix}=\begin{pmatrix} 1&0&0\\ \ell_{21}&1&0\\ \ell_{31}&\ell_{32}&1\end{pmatrix}\begin{pmatrix} u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\ 0&0&u_{33}\end{pmatrix}$$ then $A$ is overwritten with $$ \begin{pmatrix} u_{11}&u_{12}&u_{13}\\ \ell_{21}&u_{22}&u_{23}\\ \ell_{31}&\ell_{32}&u_{33}\end{pmatrix}. $$ The code uses {{array-assign!}}, {{specialized-array-share}}, {{array-extract}}, and {{array-outer-product}}.
2044
2045<enscript highlight="scheme">(define (LU-decomposition A)
2046  ;; Assumes the domain of A is [0,n)\times [0,n)
2047  ;; and that Gaussian elimination can be applied
2048  ;; without pivoting.
2049  (let ((n
2050         (interval-upper-bound (array-domain A) 0))
2051        (A_
2052         (array-getter A)))
2053    (do ((i 0 (fx+ i 1)))
2054        ((= i (fx- n 1)) A)
2055      (let* ((pivot
2056              (A_ i i))
2057             (column/row-domain
2058              ;; both will be one-dimensional
2059              (make-interval (vector (+ i 1))
2060                             (vector n)))
2061             (column
2062              ;; the column below the (i,i) entry
2063              (specialized-array-share A
2064                                       column/row-domain
2065                                       (lambda (k)
2066                                         (values k i))))
2067             (row
2068              ;; the row to the right of the (i,i) entry
2069              (specialized-array-share A
2070                                       column/row-domain
2071                                       (lambda (k)
2072                                         (values i k))))
2073
2074             ;; the subarray to the right and
2075             ;; below the (i,i) entry
2076             (subarray
2077              (array-extract
2078               A (make-interval
2079                  (vector (fx+ i 1) (fx+ i 1))
2080                  (vector n         n)))))
2081        ;; Compute multipliers.
2082        (array-assign!
2083         column
2084         (array-map (lambda (x)
2085                      (/ x pivot))
2086                    column))
2087        ;; Subtract the outer product of i'th
2088        ;; row and column from the subarray.
2089        (array-assign!
2090         subarray
2091         (array-map -
2092                    subarray
2093                    (array-outer-product * column row)))))))</enscript>
2094We then define a $4\times 4$ [[https://en.wikipedia.org/wiki/Hilbert_matrix|Hilbert matrix]]:
2095
2096<enscript highlight="scheme">(define A
2097  (array-copy
2098   (make-array (make-interval '#(4 4))
2099               (lambda (i j)
2100                 (/ (+ 1 i j))))))
2101
2102(define (array-display A)
2103
2104  (define (display-item x)
2105    (display x) (display "\t"))
2106
2107  (newline)
2108  (case (array-dimension A)
2109    ((1) (array-for-each display-item A) (newline))
2110    ((2) (array-for-each (lambda (row)
2111                           (array-for-each display-item row)
2112                           (newline))
2113                         (array-curry A 1)))
2114    (else
2115     (error "array-display can't handle > 2 dimensions: " A))))
2116
2117(display "\nHilbert matrix:\n\n")
2118
2119(array-display A)
2120
2121;;; which displays:
2122;;; 1       1/2     1/3     1/4
2123;;; 1/2     1/3     1/4     1/5
2124;;; 1/3     1/4     1/5     1/6
2125;;; 1/4     1/5     1/6     1/7
2126
2127(LU-decomposition A)
2128
2129(display "\nLU decomposition of Hilbert matrix:\n\n")
2130
2131(array-display A)
2132
2133;;; which displays:
2134;;; 1       1/2     1/3     1/4
2135;;; 1/2     1/12    1/12    3/40
2136;;; 1/3     1       1/180   1/120
2137;;; 1/4     9/10    3/2     1/2800
2138</enscript>
2139We can now define matrix multiplication as follows to check our result:
2140
2141<enscript highlight="scheme">;;; Functions to extract the lower- and upper-triangular
2142;;; matrices of the LU decomposition of A.
2143
2144(define (L a)
2145  (let ((a_ (array-getter a))
2146        (d  (array-domain a)))
2147    (make-array
2148     d
2149     (lambda (i j)
2150       (cond ((= i j) 1)        ;; diagonal
2151             ((> i j) (a_ i j)) ;; below diagonal
2152             (else 0))))))      ;; above diagonal
2153
2154(define (U a)
2155  (let ((a_ (array-getter a))
2156        (d  (array-domain a)))
2157    (make-array
2158     d
2159     (lambda (i j)
2160       (cond ((<= i j) (a_ i j)) ;; diagonal and above
2161             (else 0))))))       ;; below diagonal
2162
2163(display "\nLower triangular matrix of decomposition of Hilbert matrix:\n\n")
2164(array-display (L A))
2165
2166;;; which displays:
2167;;; 1       0       0       0
2168;;; 1/2     1       0       0
2169;;; 1/3     1       1       0
2170;;; 1/4     9/10    3/2     1
2171
2172
2173(display "\nUpper triangular matrix of decomposition of Hilbert matrix:\n\n")
2174(array-display (U A))
2175
2176;;; which displays:
2177;;; 1       1/2     1/3     1/4
2178;;; 0       1/12    1/12    3/40
2179;;; 0       0       1/180   1/120
2180;;; 0       0       0       1/2800
2181
2182;;; We'll define a brief, not-very-efficient matrix multiply routine.
2183
2184(define (array-dot-product a b)
2185  (array-fold + 0 (array-map * a b)))
2186
2187(define (matrix-multiply a b)
2188  (let ((a-rows
2189         (array-copy (array-curry a 1)))
2190        (b-columns
2191         (array-copy (array-curry (array-rotate b 1) 1))))
2192    (array-outer-product array-dot-product a-rows b-columns)))
2193
2194;;; We'll check that the product of the result of LU
2195;;; decomposition of A is again A.
2196
2197(define product (matrix-multiply (L A) (U A)))
2198
2199(display "\nProduct of lower and upper triangular matrices \n")
2200(display "of LU decomposition of Hilbert matrix:\n\n")
2201(array-display product)
2202
2203;;; which displays:
2204;;; 1       1/2     1/3     1/4
2205;;; 1/2     1/3     1/4     1/5
2206;;; 1/3     1/4     1/5     1/6
2207;;; 1/4     1/5     1/6     1/7</enscript>
2208'''Inner products.''' One can define an APL-style inner product as
2209
2210<enscript highlight="scheme">(define (inner-product A f g B)
2211  (array-outer-product
2212   (lambda (a b)
2213     (array-reduce f (array-map g a b)))
2214   (array-copy (array-curry A 1))
2215   (array-copy (array-curry (array-rotate B 1) 1))))</enscript>
2216This routine differs from that found in APL in several ways: The arguments {{A}} and {{B}} must each have two or more dimensions, and the result is always an array, never a scalar.
2217
2218We take some examples from the [[https://www.dyalog.com/uploads/aplx/APLXLangRef.pdf|APLX Language Reference]]:
2219
2220<enscript highlight="scheme">;; Examples from
2221;; http://microapl.com/apl_help/ch_020_020_880.htm
2222
2223(define TABLE1
2224  (list->array
2225   '(1 2
2226     5 4
2227     3 0)
2228   (make-interval '#(3 2))))
2229
2230(define TABLE2
2231  (list->array
2232   '(6 2 3 4
2233     7 0 1 8)
2234   (make-interval '#(2 4))))
2235
2236(array-display (inner-product TABLE1 + * TABLE2))
2237
2238;;; Displays
2239;;; 20      2       5       20
2240;;; 58      10      19      52
2241;;; 18      6       9       12
2242
2243(define X   ;; a "row vector"
2244  (list->array '(1 3 5 7) (make-interval '#(1 4))))
2245
2246(define Y   ;; a "column vector"
2247  (list->array '(2 3 6 7) (make-interval '#(4 1))))
2248
2249(array-display (inner-product X + (lambda (x y) (if (= x y) 1 0)) Y))
2250
2251;;; Displays
2252;;; 2</enscript>
2253=== Acknowledgments
2254
2255The SRFI author thanks Edinah K Gnang, John Cowan, Sudarshan S Chawathe, Jamison Hope, and Per Bothner for their comments and suggestions, and Arthur A. Gleckler, SRFI Editor, for his guidance and patience.
2256
2257=== References
2258
2259# [[https://groups.google.com/forum/?hl=en#!msg/comp.lang.scheme/7nkx58Kv6RI/a5hdsduFL2wJ|"multi-dimensional arrays in R5RS?"]], by Alan Bawden.
2260# [[https://srfi.schemers.org/srfi-4/|SRFI 4: Homogeneous Numeric Vector Datatypes]], by Marc Feeley.
2261# [[https://srfi.schemers.org/srfi-25/|SRFI 25: Multi-dimensional Array Primitives]], by Jussi Piitulainen.
2262# [[https://srfi.schemers.org/srfi-47/|SRFI 47: Array]], by Aubrey Jaffer.
2263# [[https://srfi.schemers.org/srfi-58/|SRFI 58: Array Notation]], by Aubrey Jaffer.
2264# [[https://srfi.schemers.org/srfi-63/|SRFI 63: Homogeneous and Heterogeneous Arrays]], by Aubrey Jaffer.
2265# [[https://srfi.schemers.org/srfi-164/|SRFI 164: Enhanced multi-dimensional Arrays]], by Per Bothner.
2266
2267=== Egg Maintainer
2268
2269Diego A. Mundo
2270
2271=== Copyright
2272
2273 Â© 2016, 2018, 2020 Bradley J Lucier. All Rights Reserved.
2274
2275 Permission is hereby granted, free of charge, to any person obtaining a copy of
2276 this software and associated documentation files (the "Software"), to deal in
2277 the Software without restriction, including without limitation the rights to
2278 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
2279 of the Software, and to permit persons to whom the Software is furnished to do
2280 so, subject to the following conditions:
2281
2282 The above copyright notice and this permission notice (including the next
2283 paragraph) shall be included in all copies or substantial portions of the
2284 Software.
2285
2286 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
2287 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
2288 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
2289 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
2290 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
2291 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
2292 SOFTWARE.
Note: See TracBrowser for help on using the repository browser.