by Bradley J. Lucier

This SRFI is currently in *final* status. Here is an explanation of each status that a SRFI can hold. To provide input on this SRFI, please send email to `srfi-231@nospamsrfi.schemers.org`

. To subscribe to the list, follow these instructions. You can access previous messages via the mailing list archive.

- Received: 2022-01-05
- Draft #1 published: 2022-01-07
- Draft #2 published: 2022-01-20
- Draft #3 published: 2022-01-26
- Draft #4 published: 2022-02-24
- Draft #5 published: 2022-03-16
- Draft #6 published: 2022-03-17
- Draft #7 published: 2022-04-25
- Draft #8 published: 2022-05-23
- Draft #9 published: 2022-05-26
- Draft #10 published: 2022-06-01
- Draft #11 published: 2022-06-09
- Draft #12 published: 2022-06-20
- Draft #13 published: 2022-06-24
- Draft #14 published: 2022-07-19
- Draft #15 published: 2022-08-12
- Draft #16 published: 2022-08-16
- Draft #17 published: 2022-09-17
- Draft #18 published: 2022-09-22
- Draft #19 published: 2022-09-23
- Finalized: 2022-09-25
- Revised to fix errata:
- 2023-01-16 (Add note about invoking the continuations of getters and setters in
`array-assign!`

more than once.)

- 2023-01-16 (Add note about invoking the continuations of getters and setters in

This 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 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 Basic Concepts of Mathematics.) Specialized variants of arrays provide portable programs with efficient representations for common use cases.

This is a revised and improved version of SRFI 179.

- Rationale
- Overview
- Notes
- Specification
- Preliminary notes
- Miscellaneous procedures
- Intervals
- Storage classes
- Arrays
- Implementation
- Relationship to other array libraries
- Other examples
- Acknowledgments
- References
- Copyright

This SRFI was motivated by a number of somewhat independent notions, which we outline here and which are explained below.

- 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. - Provide a
**single, efficient implementation for dense arrays**(which we call*specialized arrays*). - Provide
**useful array transformations**. - Separate
**the procedures that specify the work to be done**(`array-map`

,`array-outer-product`

, etc.) from**the procedures that actually do the work**(`array-copy`

,`array-assign!`

,`array-fold-left`

, etc.). This approach**avoids temporary intermediate arrays**in computations. - Encourage
**bulk processing of arrays**rather than word-by-word operations.

This SRFI differs from the finalized SRFI 179 in the following ways:

- Empty and zero-dimensional arrays are incorporated into this SRFI. It is an error to call the setter or getter of an empty array. Zero-dimensional arrays specify their sole element with an empty multi-index.
`specialized-array-default-safe?`

and`specialized-array-default-mutable?`

are now SRFI 39 parameters.`list->array`

is now called as`(list->array interval list ...)`

; i.e., the order of the first two arguments has been reversed.`array-copy`

no longer allows changing the domain of the result. Use`(specialized-array-reshape (array-copy ...)`

instead.`new-domain`)`make-specialized-array`

now accepts an optional initial value with which to fill the new array.- The SRFI 179 procedure
`array-fold`

has been replaced by`array-fold-left`

. Now`array-fold-left`

and`array-fold-right`

follow the definition of the left and right folds in R6RS (as well as Ocaml and Haskell).`array-fold-left`

from this SRFI has different semantics to`array-fold`

from SRFI 179. `array-assign!`

now requires that the source and destination have the same domain. Use`specialized-array-reshape`

on the destination array to mimic the SRFI 179 version.- If the first argument to
`array-copy`

is a specialized array, then omitted arguments are taken from the argument array and do not default to`generic-storage-class`

,`(specialized-array-default-mutable?)`

, and`(specialized-array-default-safe?)`

. Thus, by default,`array-copy`

makes a true copy of a specialized array. - Procedures that generate useful permutations have been added:
`index-rotate`

,`index-first`

,`index-last`

, and`index-swap`

. `interval-rotate`

and`array-rotate`

have been removed; use`(array-permute A (index-rotate (array-dimension A) k))`

instead of`(array-rotate A k)`

.`array-tile`

is now more flexible in how you can decompose an array.`array-elements-in-order?`

has been renamed`array-packed?`

.`interval-cartesian-product`

can now take zero arguments, in which case it returns`(make-interval '#())`

.`char-storage-class`

is provided.- Introduced new procedures
`interval-width`

,`interval-widths`

,`interval-empty?`

,`interval-fold-left`

,`interval-fold-right`

,`storage-class-data?`

,`storage-class-data->body`

,`array-empty?`

,`make-specialized-array-from-data`

,`vector->array`

,`array->vector`

,`list*->array`

,`array->list*`

,`vector*->array`

,`array->vector*`

,`array-inner-product`

,`array-stack`

,`array-append`

,`array-block`

,`array-freeze!`

, and`array-decurry`

. - The sample implementation now provides call/cc-safe implementations of all procedures whose names do not end with
`!`

. We define new procedures`array-copy!`

,`array-stack!`

,`array-decurry!`

,`array-append!`

, and`array-block!`

, which are not guaranteed to be call/cc-safe but which may be faster or use less memory than the corresponding call/cc-safe versions. See further discussion, with our definition of "call/cc-safe", in the notes below. - A new set of "Introductory remarks" surveys some of the more important procedures in this SRFI.

The next few sections talk perhaps too much about the mathematical ideas that underpin many of the procedures in this SRFI, so I discuss here some of the procedures and compare them to operations on spreadsheets, matrices, and imaging.

There are two procedures that simply create new arrays:

`make-array`

: Takes as arguments a specification of the valid indices $i\ j\ k$ etc. of the array, together with a Scheme procedure, which, when presented with indices in the valid range, computes the array element. The elements of the array are not precomputed and stored somewhere; the specified procedure is recalculated each time that element is needed. A procedure that modifies which element is returned at a given set of indices is allowed as a third argument. See the sparse matrix example below to see how this is useful. We call the result a*generalized array*.`make-specialized-array`

: Takes as an argument a specification of a valid range of indices and reserves a block of memory in which to store elements of the matrix; optionally, one can restrict which objects can be stored as elements in the array or generate code to precheck that all the indices are in range on each access, and to precheck that values stored as array elements actually comply with any given restrictions. Elements are stored in row-major order, as in C. We call the result a*specialized array*.

In the next group of procedures, the new and old arrays share elements, so modifications to one affect the others. Also, none of these procedures move any data: for specialized arrays they just change how the data are indexed, while for generalized arrays they manipulate the arguments of the getter and setter. For specialized arrays, these procedures can be combined in any way without increasing unreasonably the number of operations required to access an array element.

`array-extract`

: Constructs a rectangular "window" or "view" into an existing array, like a rectangular region of a spreadsheet, or a submatrix of a matrix.`array-translate`

: Slides an array around, like changing the zero-based indexing of C arrays to the 1-based indexing of Fortran arrays. If you wanted to compare two subimages of the same number of rows and columns of pixels, for example, you could use array-extract to select each of the subimages, and then use array-translate to overlay one on the other, i.e., to use the same indexing for both.`array-permute`

: Swaps rows, columns, sheets, etc., of the original array, like swapping rows and columns in a spreadsheet or transposing a matrix. The auxiliary procedures`index-rotate`

,`index-first`

,`index-last`

, and`index-swap`

create commonly used permutations.`array-reverse`

: Reverses the order of rows or columns (or both) of a spreadsheet. Like flipping an image vertically or horizontally.`array-sample`

: Accesses every second (or third, etc.) row or column, or both, of an array.

The following two procedures decompose arrays in different ways. They return *generalized arrays* whose elements are themselves arrays. Like the procedures described immediately above, the resulting subarrays share elements with their argument.

`array-curry`

: Represents a $d$-dimensional array as a $d'$-dimensional array whose entries are themselves arrays of dimension $d-d'$. Like thinking of a three-dimensional CT scan as a one-dimensional array of two-dimensional slices, or thinking of a matrix as a one-dimensional array of one-dimensional rows. You could combine this operation with`array-permute`

to think of a matrix as an array of columns, or look at slices in different orientations of a three-dimensional CT scan. Considering a video as a one-dimensional sequence (in time) of two-dimensional stills (in space) is another example of currying. The subarrays share elements with the original array. The procedures`array-decurry`

and`array-stack`

, described below, reverse this process.`array-tile`

: Decomposes a $d$-dimensional array into $d$-dimensional sub-blocks with cuts parallel to the coordinate axes, and returns the subarrays in an array. Like breaking a large matrix into smaller matrices for block matrix operations. The subarrays share elements with the original array. The procedures`array-block`

and`array-append`

, described below, reverse this process.

The next few procedures set up operations to be executed in the future. They build *generalized* arrays.

`array-map`

: Specifies an operation to be applied componentwise on arrays, so if

and`A`

are matrices,`B``(array-map +`

sets up a new generalized array that adds elements of the arrays componentwise. You can chain these operations, so have`A``B`)`(array-map + (array-map (lambda (`

without immediately computing and storing all the values of those arrays.`x`) (*`alpha x`))`A`)`B`)`array-outer-product`

: Applies an operation to all possible pairs of elements of two original arrays. Like considering an $m$-vector as a column vector and an $n$-vector as a row vector, and multiplying them together to compute an $m\times n$ matrix.`array-inner-product`

: Like APL's inner product; multiplying two matrices is an example of an operation implemented using`array-inner-product`

.

Then, there are procedures that *do* generate all elements of an array and either store them somewhere, or combine them in some way:

`array-copy`

: Evaluates the argument array at all valid indices and stores those values into a new specialized array.`array-assign!`

: Evaluates the argument array at all valid indices and assigns their values to the elements of an existing array. In the Gaussian Elimination example below, we combine`array-map`

,`array-outer-product`

,`array-extract`

, and`array-assign!`

to do one step of the elimination.`array-stack`

: Like taking the individually rendered frames of an animated movie and combining them in time to make a complete video. Can be considered a partial inverse to`array-curry`

. Returns a specialized array.`array-decurry`

: Takes a "curried" array of arrays, and returns a single array with the same elements. An inverse to`array-curry`

; a multi-dimensional version of`array-stack`

.`array-append`

: Like concatenating a number of images left to right, or top to bottom. Returns a specialized array. A partial inverse to`array-tile`

.`array-block`

: Assumes that an array has been decomposed into blocks by cuts perpendicular to each coordinate axis; takes an array of those blocks as an argument, and returns a reconstructed array. An inverse to`array-tile`

; a multi-dimensional version of array-append.`array-fold-left`

,`array-fold-right`

,`array-reduce`

,`array-for-each`

,`array-any`

, and`array-every`

: Evaluates all elements of an array (for`array-every`

and`array-any`

, as many as needed to know the result), and combines them in certain ways.

Finally, we have procedures that convert between other data and arrays:

`make-specialized-array-from-data`

: Construct a specialized array whose body shares elements with an existing data structure.`array->list`

,`list->array`

,`array->vector`

, and`vector->array`

: Either transfer the elements of an array to a list or vector, or construct a specialized array from the elements of a list or vector.`array->list*`

,`list*->array`

,`array->vector*`

, and`vector*->array`

: Either transfer the elements of an array to a nested list or vector, or construct a specialized array from the elements of a nested list or vector.

I hope this brief discussion gives a flavor for the design of this SRFI.

In a 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 procedure `array-set!`

puts 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 procedure `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 procedure 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.

Mathematically, 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$ if and only if the least coordinate $k$ where $\vec i$ and $\vec j$ differ satisfies $\vec i_k<\vec j_k$. This ordering of multi-indices is also known as row-major order, which is used in the programming language C to order the elements of multi-dimensional arrays. In contrast, the programming language Fortran uses column-major order to order the elements of multi-dimensional arrays.

In `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))$.

A 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.

We incorporate Bawden-style arrays into this SRFI, but extend them in one minor way that we find useful.

We introduce the notion of a *storage class*, an object that contains procedures that manipulate, store, check, etc., different types of values. The `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.

We 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.

Requiring 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:

**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.**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 dimension. We call this*tiling*the array. The procedure`array-tile`

returns a new array, each entry of which is a subarray extracted (in the sense of`array-extract`

) from the input array.**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.**Permuting the coordinates of an array:**If $\pi$ 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. Several procedures build commonly used permutations:`(index-rotate`

rotates`n``k`)

indices`n`

places to the left;`k``(index-first`

moves the $k$th of $n$ indices to be the first index, leaving the others in order;`n``k`)`(index-last`

moves the $k$th of $n$ indices to be the last index, again leaving the others in order;`n``k`)`(index-swap`

swaps the $i$th and $j$th of $n$ indices, again leaving the others in order.`n i j`)**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.**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).**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.

We make several remarks. First, all these operations could have been computed by specifying the particular mapping $T_{BA}$ explicitly, so that these procedures are in some sense "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.

Bawden-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})$, $l_k<u_k$ — 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.

Since these two things are often sufficient for certain algorithms, we introduce in this SRFI a minimal set of interfaces for dealing with such arrays.

We also consider as Scheme objects the case when $d>0$ and some $l_k=u_k$; in this case the mathematical cross product is empty, and arrays with such a domain have no elements but still "dimension" $d$. Applying the function associated with this array is an error.

Finally, we allow $d=0$; such an array would have one element, and the function that accesses it is a function with no arguments (i.e., a "thunk").

So an array specifies a 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`

.

If 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" procedure, i.e., the value it returns may not depend solely on the arguments passed to the getter.

In 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, or may read the values from a file, etc.

In this SRFI, Bawden-style arrays are called *specialized*. A specialized array may be either mutable or immutable.

Even if an array $A$ is not a specialized array, 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))$.

One 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))$ could be computed with the overhead of computing a single affine transformation.

So, 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 procedure `specialized-array-share`

, not a more general `array-share`

.

Certain ways of sharing generalized arrays, however, are relatively easy to code and are not expensive. If we denote `(array-getter `

by `A`)

, then if `A_`

is the result of `B``array-extract`

applied to

, then `A``(array-getter `

is simply `B`)

. Similarly, if `A_`

is a two-dimensional array, and `A`

is derived from `B`

by applying the permutation $\pi((i,j))=(j,i)$, then `A``(array-getter `

is `B`)`(lambda (i j) (`

. Translation and currying also lead to transformed arrays whose getters are relatively efficiently derived from `A_` j i))

, at least for arrays of small dimension.`A_`

Thus, 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 procedures `array-extract`

, `array-translate`

, `array-permute`

, `array-curry`

, `array-reverse`

, `array-tile`

, and `array-sample`

, and we provide relatively efficient implementations of these procedures for arrays of dimension no greater than four.

If

is an array, then we generally define `A`

to be `A_``(array-getter `

and `A`)

to be `A!``(array-setter `

. The latter notation is motivated by the general Scheme convention that the names of procedures that modify the contents of data structures end in `A`)

, while the notation for the getter of an array is motivated by the TeX notation for subscripts. See particularly the Haar transform example.`!`

**Empty and zero-dimensional arrays:**The vectors of upper and lower bounds of an interval can have zero elements, in which case the zero-dimensional interval itself has no elements, but zero-dimensional arrays with this domain have getters and setters that take zero indices as arguments, and which return or set a single element, much like a Scheme`box`

. If an interval has at least one upper and lower bound, and at least one of these upper bounds equals the associated lower bound, then that interval is empty, and arrays with empty intervals as domains have getters and setters that should raise an exception when called.**This SRFI and**The Scheme procedure`call-with-current-continuation`

:`call-with-current-continuation`

captures and encapsulates as a procedure the continuation of the current computation, which, perforce, includes a certain amount of state that consists of the values of captured variables at the point the continuation is captured. This captured procedure can be invoked multiple times, as any procedure can.

No procedure in the sample implementation itself calls`call-with-current-continuation`

, but the procedural arguments to, e.g.,`make-array`

,`specialized-array-share`

,`array-map`

, etc., may themselves call`call-with-current-continuation`

.

All procedures in the sample implementation whose names do not end with an exclamation point (!) are written in a way that does not modify the state of any data captured by a continuation. We call such procedures*call/cc safe*.

It is intended that all procedures in this library whose names do not end with an exclamation point (!) be implemented in a call/cc-safe way.

Besides the procedures`array-set!`

and`array-assign!`

, which explicitly mutate state and which, therefore, are*not*call/cc safe, we provide the procedures`array-copy!`

,`array-stack!`

,`array-decurry!`

,`array-append!`

, and`array-block!`

, which are not necessarily call/cc safe, but which may be faster or use less memory than the corresponding call/cc-safe versions.

For the procedures`array-copy!`

,`array-stack!`

,`array-decurry!`

,`array-append!`

, and`array-block!`

, it is an error if the continuation of each call to an array's getter is invoked more than once. For the procedure`array-assign!`

, it is an error if the continuation of each call to the destination's setter or the source's getter is invoked more than once.**Relationship to 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 array proposal for Haskell. Our "specialized" arrays are related to Racket's "strict" arrays.**Indexers.**The argument

to`new-domain->old-domain``specialized-array-share`

is, conceptually, the getter of a multi-valued array.**Source of procedure names.**The procedure`array-curry`

gets its name from the 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 procedure of the interval, encapsulated here as`interval-contains-multi-index?`

.**Choice of procedures on intervals.**The choice of procedures for both arrays and intervals was motivated almost solely by what I needed for arrays.**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.

- Miscellaneous Procedures
- translation?, permutation?, index-rotate, index-first, index-last, index-swap.
- Intervals
- make-interval, interval?, interval-dimension, interval-lower-bound, interval-upper-bound, interval-width, interval-lower-bounds->list, interval-upper-bounds->list, interval-lower-bounds->vector, interval-upper-bounds->vector, interval=, interval-widths, interval-volume, interval-empty?, interval-subset?, interval-contains-multi-index?, interval-projections, interval-for-each, interval-fold-left, interval-fold-right, interval-dilate, interval-intersect, interval-translate, interval-permute, interval-scale, interval-cartesian-product.
- Storage Classes
- make-storage-class, storage-class?, storage-class-getter, storage-class-setter, storage-class-checker, storage-class-maker, storage-class-copier, storage-class-length, storage-class-default, storage-class-data?, storage-class-data->body, generic-storage-class, char-storage-class, s8-storage-class, s16-storage-class, s32-storage-class, s64-storage-class, u1-storage-class, u8-storage-class, u16-storage-class, u32-storage-class, u64-storage-class, f8-storage-class, f16-storage-class, f32-storage-class, f64-storage-class, c64-storage-class, c128-storage-class.
- Arrays
- specialized-array-default-safe?, specialized-array-default-mutable?, make-array, array?, array-domain, array-getter, array-dimension, mutable-array?, array-setter, array-freeze!, array-empty?, make-specialized-array, make-specialized-array-from-data, specialized-array?, array-storage-class, array-indexer, array-body, array-safe?, array-packed?, specialized-array-share, array-copy, array-copy!, array-curry, array-extract, array-tile, array-translate, array-permute, array-reverse, array-sample, array-outer-product, array-inner-product, array-map, array-for-each, array-fold-left, array-fold-right, array-reduce, array-any, array-every, array->list, list->array, array->list*, list*->array, array->vector, vector->array, vector*->array, array->vector*, array-assign!, array-stack, array-stack!, array-decurry, array-decurry!, array-append, array-append!, array-block, array-block!, array-ref, array-set!, specialized-array-reshape.

We use `take`

and `drop`

from SRFI 1 to define various procedures.

This 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.
We provide four procedures that return useful permutations.

**Procedure: **`translation? `

`object`

Returns `#t`

if

is a translation, and `object``#f`

otherwise.

**Procedure: **`permutation? `

`object`

Returns `#t`

if

is a permutation, and `object``#f`

otherwise.

**Procedure: **`index-rotate `

`n` `k`

Assumes that `n` is a nonnegative exact integer and that `k` is an exact integer between 0 and `n` (inclusive). Returns a permutation that rotates `n` indices `k` places to the left:

```
(define (index-rotate n k)
(let ((identity-permutation (iota n)))
(list->vector (append (drop identity-permutation k)
(take identity-permutation k)))))
```

For example, `(index-rotate 5 3)`

returns `'#(3 4 0 1 2)`

. It is an error of the arguments do not satisfy these conditions.

**Procedure: **`index-first `

`n` `k`

Assumes that `n` is a positive exact integer and that `k` is an exact integer between 0 (inclusive) and `n` (exclusive). Returns a permutation of length `n` that moves index `k` (with count beginning at 0) to be first and leaves the other indices in order:

```
(define (index-first n k)
(let ((identity-permutation (iota n)))
(list->vector (cons k
(append (take identity-permutation k)
(drop identity-permutation (fx+ k 1)))))))
```

For example, `(index-first 5 3)`

returns `'#(3 0 1 2 4)`

. It is an error if the arguments do not satisfy these conditions.

**Procedure: **`index-last `

`n` `k`

Assumes that `n` is a positive exact integer and that `k` is an exact integer between 0 (inclusive) and `n` (exclusive). Returns a permutation of length `n` that moves index `k` (with count beginning at 0) to be last and leaves the other indices in order:

```
(define (index-last n k)
(let ((identity-permutation (iota n)))
(list->vector (append (take identity-permutation k)
(drop identity-permutation (fx+ k 1))
(list k)))))
```

For example, `(index-last 5 3)`

returns `'#(0 1 2 4 3)`

. It is an error if the arguments do not satisfy these conditions.

**Procedure: **`index-swap `

`n` `i` `j`

Assumes that

is a positive exact integer and that `n`

and `i`

are exact integers between 0 (inclusive) and `j`

(exclusive). Returns a permutation of length `n`

that swaps index `n`

and index `i`

and leaves the other indices in order.`j`

For example, `(index-swap 5 3 0)`

returns `#(3 1 2 0 4)`

. It is an error if the arguments do not satisfy these assumptions.

An 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 exact integers. The nonnegative integer $d$ is the *dimension*
of the interval.

If $l_k=u_k$ for some $k$ then the interval is *empty*; if $d=0$ then the interval is *zero-dimensional*. So rather than mathematical objects, it is perhaps better to think of intervals as pairs of vectors $L$ and $U$ for which $L_k\leq U_k$; $L$ or $U$ could be empty (hence $d=0$).

Intervals are a data type distinct from other Scheme data types.

**Procedure: **`make-interval `

`arg1` #!optional `arg2`

Create a new interval. Assumes that

and `arg1`

(if given) are vectors (of the same length) of exact integers.`arg2`

If

is not given, then the entries of `arg2`

, if any, must be nonnegative, and they are taken as the `arg1`

of the interval, and `upper-bounds`

is set to a vector of the same length with exact zero entries.`lower-bounds`

If

is given, then `arg2`

is taken to be `arg1`

and `lower-bounds`

is taken to be `arg2`

, which must satisfy`upper-bounds`

`(<= (vector-ref ``lower-bounds` i) (vector-ref `upper-bounds` i))

for
$0\leq i<{}$`(vector-length `

. It is an error if
`lower-bounds`)

and `lower-bounds`

do not satisfy these conditions.`upper-bounds`

**Example: **

```
(let ((A (make-interval '#(3 4)))
(B (make-interval '#(0 0) '#(3 4))))
(interval= A B)) ;; => #t
```

**Procedure: **`interval? `

`object`

Returns `#t`

if

is an interval, and `object``#f`

otherwise.

**Example: **

```
(let ((A (make-interval '#(3 4)))
(B 1))
(interval? A) ;; => #t
(interval? B)) ;; => #f
```

**Procedure: **`interval-dimension `

`interval`

Assumes

is an interval; if `interval`

is built with `interval`

`(make-interval ``lower-bounds` `upper-bounds`)

then `interval-dimension`

returns `(vector-length `

.`lower-bounds`)

It is an error to call `interval-dimension`

if

is not an interval.`interval`

**Example: **

```
(let ((A (make-interval '#(3 4)))
(B (make-interval '#())))
(interval-dimension A) ;; => 2
(interval-dimension B)) ;; => 0
```

**Procedure: **`interval-lower-bound `

`interval` `i`

**Procedure: **`interval-upper-bound `

`interval` `i`

**Procedure: **`interval-width `

`interval` `i`

Assumes that

is an interval made, e.g., with `interval``(make-interval `

, and that `lower-bounds upper-bounds`)

is an exact integer that satisfies`i`

`$0 \leq i<$ ``(interval-dimension ``interval`)

.

Then `interval-lower-bound`

returns
`(vector-ref `

, `lower-bounds` `i`)`interval-upper-bound`

returns
`(vector-ref `

, and `upper-bounds` `i`)`interval-width`

returns
`(- (vector-ref `

.`upper-bounds` `i`) (vector-ref `lower-bounds` `i`))

It is an error to call `interval-lower-bound`

, `interval-upper-bound`

, or `interval-width`

if

and `interval`

do not satisfy these conditions.`i`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4))))
(interval-lower-bound A 0) ;; => 1
(interval-upper-bound A 0) ;; => 3
(interval-width A 0)) ;; => 2
```

**Procedure: **`interval-lower-bounds->list `

`interval`

**Procedure: **`interval-upper-bounds->list `

`interval`

Assumes

is an interval built, e.g., by `interval`

`(make-interval ``lower-bounds` `upper-bounds`)

Then `interval-lower-bounds->list`

returns `(vector->list `

and `lower-bounds`)`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

does not satisfy these conditions.`interval`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4))))
(interval-lower-bounds->list A) ;; => (1 0)
(interval-upper-bounds->list A)) ;; => (3 4)
```

**Procedure: **`interval-lower-bounds->vector `

`interval`

**Procedure: **`interval-upper-bounds->vector `

`interval`

Assumes

is an interval built, e.g., with `interval`

`(make-interval ``lower-bounds` `upper-bounds`)

Then `interval-lower-bounds->vector`

returns a copy of

and `lower-bounds``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

does not satisfy these conditions.`interval`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4))))
(interval-lower-bounds->vector A) ;; => '#(1 0)
(interval-upper-bounds->vector A)) ;; => '#(3 4)
```

**Procedure: **`interval-widths `

`interval`

Assumes

is an interval built, e.g., with `interval`

`(make-interval ``lower-bounds` `upper-bounds`)

Then, assuming the existence of `vector-map`

, `interval-widths`

returns

`(vector-map - ``upper-bounds` `lower-bounds`)

It is an error to call `interval-widths`

if

does not satisfy this condition.`interval`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4))))
(interval-widths A)) ;; => '#(2 4)
```

**Procedure: **`interval-volume `

`interval`

Assumes

is an interval built, e.g., with `interval`

`(make-interval ``lower-bounds` `upper-bounds`)

Then, assuming the existence of `vector-map`

, `interval-volume`

returns

`(apply * (vector->list (vector-map - ``upper-bounds` `lower-bounds`)))

It is an error to call `interval-volume`

if

does not satisfy this condition.`interval`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4)))
(B (make-interval '#())))
(interval-volume A) ;; => 8
(interval-volume B)) ;; => 1
```

**Procedure: **`interval-empty? `

`interval`

Assumes

is an interval; returns `interval``(eqv? 0 (interval-volume `

.`interval`))

It is an error to call `interval-empty?`

if

does not satisfy this condition.`interval`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(3 4)))
(B (make-interval '#()))
(C (make-interval '#(1 0) '#(1 4))))
(interval-empty? A) ;; => #f
(interval-empty? B) ;; => #f
(interval-empty? C)) ;; => #t
```

**Procedure: **`interval= `

`interval1` `interval2`

Assumes that

and `interval1`

are intervals built, e.g., with `interval2`

`(make-interval ``lower-bounds1` `upper-bounds1`)

and

`(make-interval ``lower-bounds2` `upper-bounds2`)

respectively. Then `interval=`

returns

`(and (equal? ``lower-bounds1` `lower-bounds2`) (equal? `upper-bounds1` `upper-bounds2`))

It is an error to call `interval=`

if

or `interval1`

do not satisfy this condition.`interval2`

**Example: **

```
(let ((A (make-interval '#(1)))
(B (make-interval '#(1 1)))
(C (make-interval '#(0) '#(1)))
(D (make-interval '#(0 0)))
(E (make-interval '#(0))))
(interval= A B) ;; => #f
(interval= A C) ;; => #t
(interval= D E)) ;; => #f
```

**Procedure: **`interval-subset? `

`interval1` `interval2`

Assumes that

and `interval1`

are intervals of the same dimension $d$. Then `interval2``interval-subset?`

returns `#t`

if

`(>= (interval-lower-bound ``interval1` j) (interval-lower-bound `interval2` j))

and

`(<= (interval-upper-bound ``interval1` j) (interval-upper-bound `interval2` j))

for all $0\leq j<d$. Otherwise, it returns `#f`

. It is an error if the arguments do not satisfy these conditions.

**Example: **

```
(let ((A (make-interval '#(2 3)))
(B (make-interval '#(1 1)))
(C (make-interval '#(3 1) '#(3 3))))
(interval-subset? A B) ;; => #f
(interval-subset? B A) ;; => #t
(interval-subset? C A)) ;; => #f
```

**Procedure: **`interval-contains-multi-index? `

`interval` . `multi-index`

Assumes that

is an interval with dimension $d$ and `interval`

is a multi-index (a sequence of exact integers) of length $d$. Then `multi-index``interval-contains-multi-index?`

returns `(every <= (interval-lower-bounds->list `

.`interval`) `multi-index` (interval-upper-bounds->list `interval`))

It is an error to call `interval-contains-multi-index?`

if

and `interval`

do not satisfy this condition.`multi-index`

**Example: **

```
(let ((A (make-interval '#(1 0) '#(4 5))))
(interval-contains-multi-index? A 2 1) ;; => #t
(interval-contains-multi-index? A 0 3)) ;; => #f
```

**Procedure: **`interval-projections `

`interval` `right-dimension`

Conceptually, `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

$[l_0,u_0)\times\cdots\times[l_{d-\text{right-dimension}-1},u_{d-\text{right-dimension}-1})$

and

$[l_{d-\text{right-dimension}},u_{d-\text{right-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$

This procedure, the inverse of Cartesian products or cross products of intervals, is used to keep track of the domains of curried arrays.

More precisely, this procedure assumes that

is an interval and `interval`

is an exact integer that satisfies `right-dimension``0 <= `

, in which case `right-dimension` <= `d``interval-projections`

returns two intervals:

```
(let ((left-dimension
(- (interval-dimension interval right-dimension)))
(lowers
(interval-lower-bounds->list interval))
(uppers
(interval-upper-bounds->list interval)))
(values
(make-interval
(list->vector (take lowers left-dimension))
(list->vector (take uppers left-dimension)))
(make-interval
(list->vector (drop lowers left-dimension))
(list->vector (drop uppers left-dimension)))))
```

It is an error to call `interval-projections`

if its arguments do not satisfy these conditions.

**Example: **

```
(let ((A (make-interval '#(2 3 1 5 4))))
(call-with-values
(lambda ()
(interval-projections A 2))
(lambda (left right)
(interval= (make-interval '#(2 3 1)) left) ;; => #t
(interval= (make-interval '#(5 4)) right)))) ;; => #t
```

**Procedure: **`interval-for-each `

`f` `interval`

This procedure assumes that

is an interval and `interval`

is a procedure whose domain includes elements of `f`

. It is an error to call `interval``interval-for-each`

if

and `interval`

do not satisfy these conditions.`f`

`interval-for-each`

calls

with each multi-index of `f`

as arguments, all in lexicographical order.`interval`

In particular, if

is zero-dimensional, `interval`

is called as a thunk; if the interval is empty, then `f`

is never called.`f`

**Example: **

```
(let ((A (make-interval '#(3 2))))
(interval-for-each (lambda (i j)
(display i) (display #\space)
(display j) (display #\space)
(display "=> ")
(display (and (even? i)
(even? j)))
(newline))
A))
```

Displays:

```
0 0 => #t
0 1 => #f
1 0 => #f
1 1 => #f
2 0 => #t
2 1 => #f
```

**Procedure: **`interval-fold-left `

`f` `operator` `identity` `interval`

**Procedure: **`interval-fold-right `

`f` `operator` `identity` `interval`

These procedures assume that

is a procedure whose domain includes elements of `f`

, that `interval`

is a procedure of two arguments, and that `operator`

is an interval.`interval`

If

is empty, these procedures return `interval`

. If `identity`

is zero-dimensional, then `interval``interval-fold-left`

returns `(`

and `operator` `identity` (`f`))`interval-fold-right`

returns `(`

.`operator` (`f`) `identity`)

Otherwise, assume that there is a procedure `(next-multi-index multi-index interval)`

which, given an interval and a list representing a multi-index in that interval, returns either a list representing the next valid multi-index in the interval or `#f`

if no such multi-index exists.

Then these procedures can be defined as

```
(define (interval-fold-left f operator identity interval)
(cond ((interval-empty? interval)
identity)
((zero? (interval-dimension interval))
(operator identity (f)))
(else
(let loop ((result identity)
(multi-index (interval-lower-bounds->list interval)))
(let* ((item (apply f multi-index))
(new-result (operator result item))
(next (next-multi-index multi-index interval)))
(if next
(loop new-result next)
new-result))))))
(define (interval-fold-right f operator identity interval)
(cond ((interval-empty? interval)
identity)
((zero? (interval-dimension interval))
(operator (f) identity))
(else
(let loop ((multi-index (interval-lower-bounds->list interval)))
(if multi-index
(let* ((item (apply f multi-index))
(tail-result (loop (next-multi-index multi-index interval))))
(operator item tail-result))
identity)))))
```

Note that `interval-fold-left`

alternates evaluations of

and `f`

, while `operator``interval-fold-right`

evaluates

at all multi-indices before applying `f`

to any arguments.`operator`

It is an error if the arguments do not satisfy these assumptions.

**Example: **One can compute the Sieve of Eratosthenes with

```
(define (eratosthenes n)
;; Compute all primes <= n
(let* ((sqrt-n (exact (floor (sqrt n))))
(A (make-specialized-array (make-interval '#(2)
(vector (+ n 1)))
u1-storage-class
1))
(A_ (array-getter A))
(A! (array-setter A)))
(do ((i 2 (+ i 1)))
((> i sqrt-n)
(interval-fold-right identity
(lambda (i result)
(if (eqv? (A_ i) 1)
(cons i result)
result))
'()
(array-domain A)))
(if (eqv? (A_ i) 1)
(do ((j (square i) (+ j i)))
((> j n))
(A! 0 j))))))
(length (eratosthenes 1000000)) => 78498
```

**Procedure: **`interval-dilate `

`interval` `lower-diffs` `upper-diffs`

Assumes that

is an interval with
lower bounds $\ell_0,\dots,\ell_{d-1}$ and
upper bounds $u_0,\dots,u_{d-1}$, and `interval`

is a vector of exact integers $L_0,\dots,L_{d-1}$ and `lower-diffs`

is a vector of exact integers $U_0,\dots,U_{d-1}$. Then `upper-diffs``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
valid interval. It is an error if the arguments do not satisfy these conditions.

**Examples:**

```
(interval=
(interval-dilate (make-interval '#(100 100))
'#(1 1) '#(1 1))
(make-interval '#(1 1) '#(101 101))) => #t
(interval=
(interval-dilate (make-interval '#(100 100))
'#(-1 -1) '#(1 1))
(make-interval '#(-1 -1) '#(101 101))) => #t
(interval=
(interval-dilate (make-interval '#(100 100))
'#(0 0) '#(-50 -50))
(make-interval '#(50 50))) => #t
(interval-dilate
(make-interval '#(100 100))
'#(0 0) '#(-500 -50)) => error
```

**Procedure: **`interval-intersect `

`interval` . `intervals`

Assumes that all the arguments are intervals of the same dimension. If they have a valid intersection,
then `interval-intersect`

returns that intersection; otherwise it returns `#f`

.

More precisely, `interval-intersect`

calculates

```
(let* ((intervals (cons interval intervals))
(lower-bounds (apply vector-map max (map interval-lower-bounds intervals)))
(upper-bounds (apply vector-map min (map interval-upper-bounds intervals))))
(and (vector-every (lambda (x y) (<= x y)) lower-bounds upper-bounds)
(make-interval lower-bounds upper-bounds)))
```

It is an error if the arguments are not all intervals with the same dimension.

**Example: **

```
(let ((A (make-interval '#(2 5) '#(10 7)))
(B (make-interval '#(0 6) '#(8 11)))
(C (make-interval '#(2 6) '#(8 7)))
(D (make-interval '#(1 1))))
(interval= (interval-intersect A B) C) ;; => #t
(interval-intersect A D)) ;; => #f
```

**Procedure: **`interval-translate `

`interval` `translation`

Assumes that

is an interval, with, e.g.,
lower bounds $\ell_0,\dots,\ell_{d-1}$ and
upper bounds $u_0,\dots,u_{d-1}$, and `interval`

is a translation with entries $T_0,\dots,T_{d-1}$.
Then `translation``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.

One could define `(interval-translate interval translation)`

by `(interval-dilate interval translation translation)`

.

**Example: **

```
(let ((A (make-interval '#(2 5) '#(10 7)))
(B (make-interval '#(1 6) '#(9 8))))
(interval= (interval-translate A '#(-1 1)) B)) ;; => #t
```

**Procedure: **`interval-permute `

`interval` `permutation`

Assumes that

is an interval and `interval`

is a permutation with the same dimension as `permutation`

. It is an error if the arguments do not satisfy these conditions.`interval`

Heuristically, this procedure returns a new interval whose axes have been permuted in a way consistent with

.
But we have to say how the entries of `permutation`

are associated with the new interval.`permutation`

We 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}})$.

For 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 `

will be
the representation of $[0,16)\times [0,4)\times[0,8)\times[0,21)$.`interval` `permutation`)

**Example: **

```
(let ((A (make-interval '#(4 8 21 16)))
(B (make-interval '#(16 4 8 21))))
(interval= (interval-permute A '#(3 0 1 2)) B)) ;; => #t
```

**Procedure: **`interval-scale `

`interval` `scales`

Assumes that

is a $d$-dimensional interval $[0,u_1)\times\cdots\times[0,u_{d-1})$ with all lower bounds zero, and `interval`

is a length-$d$ vector of positive exact integers, which we'll denote by $\vec s$. Then `scales``interval-scale`

returns the interval $[0,\operatorname{ceiling}(u_1/s_1))\times\cdots\times[0,\operatorname{ceiling}(u_{d-1}/s_{d-1}))$.

It is an error if

and `interval`

do not satisfy these conditions.`scales`

**Example: **

```
(let ((A (make-interval '#(4 7)))
(B (make-interval '#(2 4))))
(interval= (interval-scale A '#(3 2)) B)) ;; => #t
```

**Procedure: **`interval-cartesian-product . `

`intervals`

Assumes that all the arguments are intervals. Implements the Cartesian product of the intervals in

. Returns:`intervals`

```
(make-interval
(list->vector
(apply append (map interval-lower-bounds->list intervals)))
(list->vector
(apply append (map interval-upper-bounds->list intervals))))
```

It is an error if any argument is not an interval.

**Example: **

```
(let ((A (make-interval '#(3 4)))
(B (make-interval '#(1 2 3) '#(7 8 9)))
(C (make-interval '#(0 0 1 2 3) '#(3 4 7 8 9))))
(interval= (interval-cartesian-product A B) C)) ;; => #t
```

Conceptually, a storage-class is a set of procedures to manage the backing store of a specialized array. The procedures allow one to make a backing store, to get values from the store, to set new values, to return the length of the store, to specify a default value for initial elements of the backing store, to recognize which data can be converted to a backing store of this storage class, and to convert data to a backing store of this storage class. Typically, a backing store is a (heterogeneous or homogeneous) vector. A storage-class has a type distinct from other Scheme types.

**Procedure: **`make-storage-class `

`getter` `setter` `checker` `maker` `copier` `length` `default` `data?` `data->body`

Here 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).

`(`

returns a linearly addressed object containing`maker n``value`)

elements of value`n`

.`value``copier`may be #f or a procedure; if a procedure then if

and`to`

were created by`from`

, then`maker``(`

copies elements from`copier to at from start end`)

beginning at`from`

(inclusive) and ending at`start`

(exclusive) to`end`

beginning at`to`

. It is assumed that all the indices involved are within the domain of`at`

and`from`

, as needed. The order in which the elements are copied is unspecified.`to`- If

is an object created by`v``(`

and 0 <=`maker n value`)

<`i`

, then`n``(`

returns the current value of the`getter v i`)

'th element of`i`

, and`v``(`

.`checker`(`getter v i`)) => #t - If

is an object created by`v``(`

, 0 <=`maker n value`)

<`i`

, and`n``(`

, then`checker``val`) => #t`(`

sets the value of the`setter v i val`)

'th element of`i`

to`v`

.`val` - If

is an object created by`v``(`

then`maker n value`)`(`

returns`length v`)

.`n` - The

and`data?`

entries are low-level procedures.`data->body``((storage-class-data?`

returns`storage-class`)`data`)`#t`

if and only if`((storage-class-data->body`

returns a body sharing data with`storage-class`)`data`)

, without copying. See the discussion of`data``make-specialized-array-from-data`

.

If the arguments do not satisfy these conditions, then it is an error to call `make-storage-class`

.

Note that we assume that

and `getter`

generally take `setter`*O*(1) time to execute.

**Procedure: **`storage-class? `

`m`

Returns `#t`

if

is a storage class, and `m``#f`

otherwise.

**Procedure: **`storage-class-getter `

`m`

**Procedure: **`storage-class-setter `

`m`

**Procedure: **`storage-class-checker `

`m`

**Procedure: **`storage-class-maker `

`m`

**Procedure: **`storage-class-copier `

`m`

**Procedure: **`storage-class-length `

`m`

**Procedure: **`storage-class-default `

`m`

**Procedure: **`storage-class-data? `

`m`

**Procedure: **`storage-class-data->body `

`m`

Assumes that

is a storage class, created, e.g., by`m`

`(make-storage-class ``getter setter checker maker copier length default data? data->body`)

Then `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``storage-class-default`

returns

, `default``storage-class-data?`

returns

, and `data?``storage-class-data->body`

returns

. Otherwise, it is an error to call any of these procedures.`data->body`

**Variable: **`generic-storage-class`

**Variable: **`char-storage-class`

**Variable: **`s8-storage-class`

**Variable: **`s16-storage-class`

**Variable: **`s32-storage-class`

**Variable: **`s64-storage-class`

**Variable: **`u1-storage-class`

**Variable: **`u8-storage-class`

**Variable: **`u16-storage-class`

**Variable: **`u32-storage-class`

**Variable: **`u64-storage-class`

**Variable: **`f8-storage-class`

**Variable: **`f16-storage-class`

**Variable: **`f32-storage-class`

**Variable: **`f64-storage-class`

**Variable: **`c64-storage-class`

**Variable: **`c128-storage-class`

`generic-storage-class`

is defined as if by

```
(define generic-storage-class
(make-storage-class vector-ref
vector-set!
(lambda (arg) #t)
make-vector
vector-copy!
vector-length
#f
vector?
values))
```

In the sample implementation `char-storage-class`

is defined as

```
(define char-storage-class
(make-storage-class string-ref
string-set!
char?
make-string
string-copy!
string-length
#\0
string?
values))
```

Implementations shall define `s`

for `X`-storage-class

=8, 16, 32, and 64 (which have default values 0 and
manipulate exact integer values between -2`X`^{X-1} and
2^{X-1}-1 inclusive),
`u`

for `X`-storage-class

=1, 8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between 0 and
2`X`^{X}-1 inclusive),
`f`

for `X`-storage-class

= 8, 16, 32, and 64 (which have default value 0.0 and manipulate 8-, 16-, 32-, and 64-bit floating-point numbers), and
`X``c`

for `X`-storage-class

= 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).`X`

Implementations 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`

.

Arrays are a data type distinct from other Scheme data types.

In the examples we use a procedure `array-unveil`

that lists the multi-indices and elements of an array in lexicographical order:

```
(define (array-unveil A)
(let ((D_A (array-domain A))
(A_ (array-getter A)))
(interval-for-each (lambda args
(for-each (lambda (arg)
(display #\space) (display arg))
args)
(for-each display
(list " => " (apply A_ args) #\newline)))
D_A)))
```

For example:

```
(let ((A (make-array (make-interval '#(2 3 2)) list)))
(array-unveil A))
```

displays:

```
0 0 0 => (0 0 0)
0 0 1 => (0 0 1)
0 1 0 => (0 1 0)
0 1 1 => (0 1 1)
0 2 0 => (0 2 0)
0 2 1 => (0 2 1)
1 0 0 => (1 0 0)
1 0 1 => (1 0 1)
1 1 0 => (1 1 0)
1 1 1 => (1 1 1)
1 2 0 => (1 2 0)
1 2 1 => (1 2 1)
```

**Parameter: **`specialized-array-default-safe?`

A parameter as specified in SRFI 39. Initially, `(specialized-array-default-safe?)`

returns `#f`

. It is an error to call `(specialized-array-default-safe? `

if `arg`)

is not a boolean.`arg`

**Parameter: **`specialized-array-default-mutable?`

A parameter as specified in SRFI 39. Initially, `(specialized-array-default-mutable?)`

returns `#t`

. It is an error to call `(specialized-array-default-mutable? `

if `arg`)

is not a boolean.`arg`

**Procedure: **`make-array `

`interval` `getter` [ `setter` ]

Assume first that the optional argument `setter`

is not given.

If

is an interval and `interval`

is a procedure from
`getter`

to Scheme objects, then `interval``make-array`

returns an array with domain

and getter `interval`

.`getter`

It is an error to call `make-array`

if

and `interval`

do not satisfy these conditions.`getter`

If now

is specified, assume that it is a procedure such that getter and setter satisfy: If`setter`

`(`

$\neq$i_{1},...,i_{n})`(`

j_{1},...,j_{n})

are elements of

and `interval`

`(getter ``j`_{1} ... `j`_{n}) => x

then "after"

`(setter v ``i`_{1} ... `i`_{n})

we have

`(getter ``j`_{1} ... `j`_{n}) => x

and

`(getter ``i`_{1},...,`i`_{n}) => v

Then `make-array`

builds a mutable array with domain

, getter `interval`

, and
setter `getter`

. It is an error to call `setter``make-array`

if its arguments do not satisfy these conditions.

**Example: **

```
(define a (make-array (make-interval '#(1 1) '#(11 11))
(lambda (i j)
(if (= i j)
1
0))))
```

defines an array for which `(array-getter a)`

returns 1 when i=j and 0 otherwise.

**Example: **

```
(define a ;; a sparse array
(let ((domain
(make-interval '#(1000000 1000000)))
(sparse-rows
(make-vector 1000000 '())))
(make-array
domain
(lambda (i j)
(cond ((assv j (vector-ref sparse-rows i))
=> cdr)
(else
0.0)))
(lambda (v i j)
(cond ((assv j (vector-ref sparse-rows i))
=> (lambda (pair)
(set-cdr! pair v)))
(else
(vector-set!
sparse-rows
i
(cons (cons j v)
(vector-ref sparse-rows i)))))))))
(define a_ (array-getter a))
(define a! (array-setter a))
(a_ 12345 6789) => 0.
(a_ 0 0) => 0.
(a! 1.0 0 0) => undefined
(a_ 12345 6789) => 0.
(a_ 0 0) => 1.
```

**Example: ** If an array `A`

is empty, e.g., `(make-array (make-interval '#(0 0)) getter setter)`

, then it is an error to call `getter`

or `setter`

. Still, such arrays can usefully exist to simplify limit cases of some algorithms.

**Example: **`(define a (make-array (make-interval '#()) (lambda () 42)))`

makes an array with a zero-dimensional domain whose getter takes no arguments and always returns 42.

**Example: **We can have the following interactive session, which builds a zero-dimensional mutable array:

```
> (define a
(let ((contents (box 42)))
(make-array
(make-interval '#())
(lambda ()
(unbox contents))
(lambda (val)
(set-box! contents val)))))
> (define a_ (array-getter a))
> (define a! (array-setter a))
> (a_)
42
> (a! 23)
> (a_)
23
```

**Procedure: **`array? `

`object`

Returns `#t`

if

is an array and `object``#f`

otherwise.

**Procedure: **`array-domain `

`array`

**Procedure: **`array-getter `

`array`

Assumes that

is an array built, e.g., by`array`

`(make-array ``interval` `getter` [`setter`])

(with or without the optional

argument). Then `setter``array-domain`

returns

and `interval``array-getter`

returns

.
It is an error to call `getter``array-domain`

or `array-getter`

if

is not an array.`array`

**Example: **

```
(define a (make-array (make-interval '#(1 1) '#(11 11))
(lambda (i j)
(if (= i j)
1
0))))
(define a_ (array-getter a))
(a_ 3 3) => 1
(a_ 2 3) => 0
(a_ 11 0) => is an error
```

**Procedure: **`array-dimension `

`array`

Shorthand for `(interval-dimension (array-domain `

. It is an error to call this procedure if `array`))

is not an array.`array`

**Example: **

```
(let ((A (make-array (make-interval '#(3 3)) list))
(B (make-array (make-interval '#()) (lambda () 42))))
(array-dimension A) ;; => 2
(array-dimension B)) ;; => 0
```

**Procedure: **`mutable-array? `

`object`

Returns `#t`

if

is a mutable array and `object``#f`

otherwise.

**Example: **

```
(let ((A (array-copy (make-array (make-interval '#(2 2)) list)
generic-storage-class
#t))
(B (array-copy (make-array (make-interval '#(2 2)) list)
generic-storage-class
#f)))
(mutable-array? A) ;; => #t
(mutable-array? B)) ;; => #f
```

**Procedure: **`array-setter `

`array`

If

is an array built by`array`

`(make-array ``interval` `getter` `setter`)

then `array-setter`

returns

. Other procedures can build mutable arrays, e.g., `setter``array-copy`

. It is an error to call `array-setter`

if

is not a mutable array.`array`

**Procedure: **`array-freeze! `

`array`

Modifies the array

so it is not mutable. Returns the modified argument.`array`

It is an error if

is not an array.`array`

**Example: **

```
(let ((array (array-copy (make-array (make-interval '#(2 2)) list)
generic-storage-class
#t)))
(mutable-array? array) ;; => #t
(array-freeze! array)
(mutable-array? array)) ;; => #f
```

**Procedure: **`array-empty? `

`array`

Assumes

is an array, and returns `array``(interval-empty? (array-domain `

. It is an error if the argument is not an array.`array`))

**Example: **

```
(let ((A (make-array (make-interval '#(2 2)) list))
(B (make-array (make-interval '#(4 0 4)) list)))
(array-empty? A) ;; => #f
(array-empty? B)) ;; => #t
```

**Procedure: **`make-specialized-array `

`interval` [ `storage-class` [ `initial-value` [ `safe?` ] ] ]

Constructs a mutable specialized array from its arguments.

must be given an interval. If given, `interval`

must be a storage class; if it is not given, it defaults to `storage-class``generic-storage-class`

. If given,

must be a value that can be manipulated by `initial-value`

; if it is not given, it defaults to `storage-class``(storage-class-default `

. If given, `storage-class`)

must be a boolean; if it is not given, it defaults to the current value of `safe?``(specialized-array-default-safe?)`

.

The body of the result is constructed as

` ((storage-class-maker ``storage-class`)
(interval-volume `interval`)
`initial-value`)

The indexer of the resulting array is constructed as the lexicographical mapping of

onto the interval `interval``[0,(interval-volume `

.`interval`))

If

is `safe``#t`

, then the arguments of the getter and setter (including the value to be stored) of the resulting array are checked for correctness.

After correctness checking (if needed), `(array-getter `

is defined simply as `array`)

```
(lambda multi-index
((storage-class-getter
````storage-class`)
(array-body `array`)
(apply (array-indexer `array`) multi-index)))

and `(array-setter `

is defined as `array`)

```
(lambda (val . multi-index)
((storage-class-setter
````storage-class`)
(array-body `array`)
(apply (array-indexer `array`) multi-index)
val))

It is an error if the arguments of `make-specialized-array`

do not satisfy these conditions.

**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

```
(define (make-u16-array interval)
(make-specialized-array interval u16-storage-class 0 #f))
```

and then simply call, e.g., `(make-u16-array (make-interval '#(3 3)))`

.

**Example: **

```
(let ((A (make-specialized-array (make-interval '#(2 3)) u8-storage-class 42)))
(array-unveil A))
```

displays:

```
0 0 => 42
0 1 => 42
0 2 => 42
1 0 => 42
1 1 => 42
1 2 => 42
```

**Procedure: **`make-specialized-array-from-data `

`data` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

This routine constructs a new specialized array using

as part of the body of the result without copying.`data`

Any missing optional arguments are assigned the values `generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

This routine exploits the low-level representation of the body of a specialized array of a specific storage class, and as such may not be portable between implementations. Here are several examples.

The sample implementation uses homogeneous vectors to represent the bodies of arrays with storage classes `u8-storage-class`

, `s8-storage-class`

, ..., `s64-storage-class`

, `f32-storage-class`

, and `f64-storage-class`

. Another implementation might use byte-vectors as the bodies of arrays for all these storage classes.

The sample implementation uses homogeneous (f32 and f64) vectors with an even number of elements to represent the bodies of arrays with storage classes `c64-storage-class`

and `c128-storage-class`

. Another implementation with purely inexact complex numbers might make another choice.

Finally, the sample implementation uses `(vector `

to represent the body of an array with a `n` (u16vector ...))`u1-storage-class`

, where

represents the valid number of bits (no more than 16 times the length of the `n``u16vector`

). A Scheme with bitvectors might choose those as the underlying representation of bodies of arrays with `u1-storage-class`

.

This routine assumes that

and `mutable?`

, if given, are booleans, and that `safe?`

, if given, is a storage class. `storage-class`

must be an object for which `data``((storage-class-data? `

returns `storage-class`) `data`)`#t`

.

This routine constructs a new one-dimensional array with storage class

, mutability `storage-class`

, safety `mutable?`

, body `safe?``((storage-class-data->body `

, with domain `storage-class`) `data`)`(make-interval (vector `

, where `N`))

is the greatest number of elements one can fit into `N`

, and indexer `data``(lambda (i) i)`

.

It is an error if the arguments do not satisfy these conditions, or if

is true and `mutable?`

is not mutable.`data`

**Example: **

```
(let ((A (make-specialized-array-from-data '#(dog cat bird))))
(array-unveil A))
```

displays:

```
0 => dog
1 => cat
2 => bird
```

**Discussion:** Correct transformations on specialized arrays *require* that the array's indexer, which maps the domain of the array to exact integers that index elements of the one-dimensional body of the array, be *affine*. The procedure `make-specialized-array-from-data`

provides a structured way to turn externally-provided data into an array with a known, very simple, one-dimensional affine indexer. With this start, the programmer can apply array transforms (e.g., `array-extract`

, `specialized-array-reshape`

, etc.) to massage the data into the shape needed.

**Example: **To build a zero-dimensional array that stores its single element in a pre-existing vector, one could use the code:

```
(pretty-print
(array->list*
(specialized-array-reshape ;; Reshape to a zero-dimensional array
(make-specialized-array-from-data ;; The basic one-dimensional array
(vector 'foo))
(make-interval '#()))))
```

prints simply

`foo`

**Example: **In the sample implementation, if you want to construct a $3\times3$ array with storage class `u1-storage-class`

from a length-one `u16vector`

named

then one could write`board`

```
(let* ((board (u16vector #b111100110111))
(A (specialized-array-reshape ;; Reshape to a 3x3 array
(array-extract ;; Only the first 9 elements
(make-specialized-array-from-data ;; The basic one-dimensional array
board u1-storage-class)
(make-interval '#(9)))
(make-interval '#(3 3))))
(B (list->array (make-interval '#(3 3)) ;; Another array with same elements
'(1 1 1
0 1 1
0 0 1)
u1-storage-class)))
(define (pad n s)
(string-append (make-string (- n (string-length s)) #\0) s))
(for-each display (list "(array-every = A B) => " (array-every = A B) #\newline))
(for-each display (list "(array-body A) => " (array-body A) #\newline))
(for-each display (list "(array-body B) => " (array-body B) #\newline))
(for-each
display
(list "(pad 16 (number->string (u16vector-ref (vector-ref (array-body A) 1) 0) 2)) => "
#\newline
(pad 16 (number->string (u16vector-ref (vector-ref (array-body A) 1) 0) 2))
#\newline))
(for-each
display
(list "(pad 16 (number->string (u16vector-ref (vector-ref (array-body B) 1) 0) 2)) => "
#\newline
(pad 16 (number->string (u16vector-ref (vector-ref (array-body B) 1) 0) 2))
#\newline)))
```

prints

(array-every = A B) => #t (array-body A) => #(16 #u16(3895)) (array-body B) => #(9 #u16(311)) (pad 16 (number->string (u16vector-ref (vector-ref (array-body A) 1) 0) 2)) => 0000111100110111 (pad 16 (number->string (u16vector-ref (vector-ref (array-body B) 1) 0) 2)) => 0000000100110111

The 9 low-order bits of board represent the entries of the array `A`

, ignoring higher order bits, and you can see the bit order that is used to represent a `u1-storage-class-body`

.

**Procedure: **`specialized-array? `

`object`

Returns `#t`

if

is a specialized-array, and `object``#f`

otherwise. A specialized-array is an array.

**Procedure: **`array-storage-class `

`array`

**Procedure: **`array-indexer `

`array`

**Procedure: **`array-body `

`array`

**Procedure: **`array-safe? `

`array`

Assumes that

is a specialized array. `array``array-storage-class`

returns the storage-class of

. `array``array-safe?`

is true if and only if the arguments of `(array-getter `

and `array`)`(array-setter `

(including the value to be stored in the array) are checked for correctness.`array`)

`(array-body `

is a linearly indexed, vector-like object (e.g., a vector, string, u8vector, etc.) indexed from 0.`array`)

`(array-indexer `

is assumed to be a one-to-one, but not necessarily onto, affine mapping from `array`)`(array-domain `

into the indexing domain of `array`)`(array-body `

.`array`)

Please see `make-specialized-array`

for how `(array-body `

, etc., are used.`array`)

It is an error to call any of these procedures if

is not a specialized array.`array`

**Procedure: **`array-packed? `

`array`

Assumes that

is a specialized array, in which case it returns `array``#t`

if the elements of

, taken in lexicographical order, are stored in `array``(array-body `

with increasing and consecutive indices, and `array`)`#f`

otherwise.

It is an error if

is not a specialized array.`array`

**Example: **

```
(let* ((A (make-specialized-array-from-data '#(0 1 2 3)))
(B (array-reverse A))
(C (array-sample A '#(2))))
(array-unveil A)
(display (array-packed? A)) (newline)
(array-unveil B)
(display (array-packed? B)) (newline) ;; adjacent boxes but decreasing index
(array-unveil C)
(display (array-packed? C)) (newline)) ;; non-adjacent boxes
```

displays:

```
0 => 0
1 => 1
2 => 2
3 => 3
#t
0 => 3
1 => 2
2 => 1
3 => 0
#f
0 => 0
1 => 2
#f
```

**Procedure: **`specialized-array-share `

`array` `new-domain` `new-domain->old-domain`

Constructs a new specialized array that shares the body of the specialized array

.
Returns an object that is behaviorally equivalent to a specialized array with the following fields:`array`

```
domain: new-domain
storage-class: (array-storage-class
````array`)
body: (array-body `array`)
indexer: (lambda multi-index
(call-with-values
(lambda ()
(apply `new-domain->old-domain`
multi-index))
(array-indexer `array`)))

The resulting array inherits its safety and mutability from

.`array`

Because

is assumed to be a one-to-one mapping, the volume of `new-domain->old-domain`

must be no greater than the number of elements of `new-domain`

.`array`

Note: It is assumed that the affine structure of the composition of

and `new-domain->old-domain``(array-indexer `

will be used to simplify:`array`)

```
(lambda multi-index
(call-with-values
(lambda ()
(apply
````new-domain->old-domain` multi-index))
(array-indexer `array`)))

It is an error if

is not a specialized array, or if `array`

is not an interval, or if `new-domain`

is not a one-to-one affine mapping from `new-domain->old-domain`

to
`new-domain``(array-domain `

.`array`)

**Example: **One can apply a "shearing" operation to an array as follows:

```
(define a
(array-copy
(make-array (make-interval '#(5 10))
list)))
(define b
(specialized-array-share
a
(make-interval '#(5 5))
(lambda (i j)
(values i (+ i j)))))
;; Print the "rows" of b
(array-for-each (lambda (row)
(pretty-print (array->list row)))
(array-curry b 1))
;; which prints
;; ((0 0) (0 1) (0 2) (0 3) (0 4))
;; ((1 1) (1 2) (1 3) (1 4) (1 5))
;; ((2 2) (2 3) (2 4) (2 5) (2 6))
;; ((3 3) (3 4) (3 5) (3 6) (3 7))
;; ((4 4) (4 5) (4 6) (4 7) (4 8))
```

This "shearing" operation cannot be achieved by combining the procedures `array-extract`

, `array-translate`

, `array-permute`

, `array-translate`

, `array-curry`

, `array-reverse`

, and `array-sample`

.

**Procedure: **`array-copy `

`array` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

**Procedure: **`array-copy! `

`array` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is an array, `array`

is a storage class that can manipulate all the elements of `storage-class`

, and `array`

and `mutable?`

are booleans.`safe?`

The specialized array returned by `array-copy`

can be defined conceptually by:

```
(list->array (array-domain array)
(array->list array)
storage-class
mutable?
safe?)
```

If

is a specialized array, then if any of `array`

, `storage-class`

, `mutable?`

are omitted, their values are assigned `safe?``(array-storage-class `

, `array`)`(mutable-array? `

, and `array`)`(array-safe? `

, respectively.`array`)

Otherwise, omitted arguments are assigned the values `generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if the arguments do not satisfy these conditions.

**Example: ** The example of reading PGM files exploits the fact that `array->list`

, and hence `array-copy`

and `array-copy!`

, evaluates an array's getter in lexicographical order.

**Example: **

```
(let* ((A (make-array (make-interval '#(2 2)) list))
(B (array-copy A)))
(display (specialized-array? A)) (newline)
(array-unveil A)
(display (specialized-array? B)) (newline)
(array-unveil B))
```

displays:

```
#f
0 0 => (0 0)
0 1 => (0 1)
1 0 => (1 0)
1 1 => (1 1)
#t
0 0 => (0 0)
0 1 => (0 1)
1 0 => (1 0)
1 1 => (1 1)
```

**Procedure: **`array-curry `

`array` `inner-dimension`

Assumes that

is an array whose domain is an interval $[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$, and `array`

is an exact integer between $0$ and $d$ (inclusive). Then `inner-dimension``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})$.

For example, if `A`

and `B`

are defined by

```
(define interval (make-interval '#(10 10 10 10)))
(define A (make-array interval list))
(define B (array-curry A 1))
(define A_ (array-getter A))
(define B_ (array-getter B))
```

so

`(A_ i j k l) => (list i j k l)`

then `B`

is an immutable array with domain `(make-interval '#(10 10 10))`

, each
of whose elements is itself an (immutable) array and

```
(equal?
(A_ i j k l)
((array-getter (B_ i j k)) l)) => #t
```

for all multi-indices `i j k l`

in `interval`

.

The subarrays are immutable, mutable, or specialized according to whether the array argument is immutable, mutable, or specialized.

More precisely, if

`0 <= ``inner-dimension` <= (interval-dimension (array-domain `array`))

then `array-curry`

returns a result as follows.

If the input array is specialized, then array-curry returns

```
(call-with-values
(lambda () (interval-projections (array-domain
````array`)
`inner-dimension`))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(specialized-array-share
`array`
inner-interval
(lambda inner-multi-index
(apply values
(append outer-multi-index
inner-multi-index))))))))

Otherwise, if the input array is mutable, then array-curry returns

```
(call-with-values
(lambda () (interval-projections (array-domain
````array`)
`inner-dimension`))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(make-array
inner-interval
(lambda inner-multi-index
(apply (array-getter `array`)
(append outer-multi-index
inner-multi-index)))
(lambda (v . inner-multi-index)
(apply (array-setter `array`)
v
(append outer-multi-index
inner-multi-index))))))))

Otherwise, array-curry returns

```
(call-with-values
(lambda () (interval-projections (array-domain
````array`)
`inner-dimension`))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(make-array
inner-interval
(lambda inner-multi-index
(apply (array-getter `array`)
(append outer-multi-index
inner-multi-index))))))))

It is an error to call `array-curry`

if its arguments do not satisfy these conditions.

If

is a specialized array, the subarrays of the result inherit their safety and mutability from `array`

.`array`

**Note: **Let's denote by

the result of `B``(array-curry `

. While the result of calling `A` `k`)`(array-getter `

is an immutable, mutable, or specialized array according to whether `B`)

itself is immutable, mutable, or specialized, `A`

is always an immutable array, where `B``(array-getter `

, which returns an array, is computed anew for each call. If `B`)`(array-getter `

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.`B`)

Please see the note in the discussion of array-tile.

**Example:**

```
(let* ((A (make-array (make-interval '#(10 10)) list))
(B (array-curry A 1)))
(array-ref A 3 4) ;; => (3 4)
(array-ref (array-ref B 3) 4)) ;; => (3 4)
```

**Example: **NumPy has the operation numpy.squeeze, which can eliminate, or "squeeze" out, all axes of an array with length 1. It can be implemented using `partition`

from SRFI 1 by

```
(define (array-squeeze a)
(call-with-values
(lambda ()
(let ((interval (array-domain a)))
(partition (lambda (k)
(eqv? (interval-width interval k) 1))
(iota (array-dimension a)))))
(lambda (ones rest)
(car (array->list
(array-curry ;; this array has exactly one element
(array-permute
a (list->vector (append ones rest))) ;; put all length-one axes at beginning
(length rest)))))))
(array->list* (array-squeeze (make-array (make-interval '#(1 2 1 2)) list)))
=>
(((0 0 0 0)
(0 0 0 1))
((0 1 0 0)
(0 1 0 1)))
(array->list*
(array-squeeze
(make-array (make-interval '#(1 2 3 4) '#(2 3 4 5))
(lambda args (apply string-append (map number->string args))))))
=>
"1234"
(array-dimension
(array-squeeze
(make-array (make-interval '#(1 2 3 4) '#(2 3 4 5))
(lambda args (apply string-append (map number->string args))))))
=>
0
(array->list*
(array-squeeze
(make-array (make-interval '#(1 2 3 4) '#(3 3 4 5))
(lambda args (apply string-append (map number->string args))))))
=>
("1234" "2234")
(array-dimension
(array-squeeze
(make-array (make-interval '#(1 2 3 4) '#(3 3 4 5))
(lambda args (apply string-append (map number->string args))))))
=>
1
```

**Procedure: **`array-extract `

`array` `interval`

Returns a new array with the same getter (and setter, if appropriate) of the first argument, defined on the second argument.

Assumes that

is an array and `array`

is an interval that is a sub-interval of `interval``(array-domain `

. If `array`)

is a specialized array, then returns `array`

` (specialized-array-share ``array`
`interval`
values)

Otherwise, if

is a mutable array, then `array``array-extract`

returns

` (make-array ``interval`
(array-getter `array`)
(array-setter `array`))

Finally, if

is an immutable array, then `array``array-extract`

returns

` (make-array ``interval`
(array-getter `array`))

It is an error if the arguments of `array-extract`

do not satisfy these conditions.

If

is a specialized array, the resulting array inherits its mutability and safety from `array`

.`array`

**Example: **

```
(let* ((A (make-array (make-interval '#(3 3)) list))
(B (array-extract A (make-interval '#(1 0) '#(3 2)))))
(display "A:\n")
(array-unveil A)
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
0 0 => (0 0)
0 1 => (0 1)
0 2 => (0 2)
1 0 => (1 0)
1 1 => (1 1)
1 2 => (1 2)
2 0 => (2 0)
2 1 => (2 1)
2 2 => (2 2)
B:
1 0 => (1 0)
1 1 => (1 1)
2 0 => (2 0)
2 1 => (2 1)
```

**Procedure: **`array-tile `

`array` `S`

Decomposes the array

into subarrays, or `array`*tiles*, specified by *cuts* perpendicular to the coordinate axes of

, which are specified by the elements second argument, `array`

, and returns an array $T$ whose elements are those tiles.`S`

If the $k$th axis of

has zero width, then the $k$th component of `array`

must be a nonempty vector of exact zeros.`array`

Otherwise, if the $k$th component of

is a positive exact integer $s$, then the cuts perpendicular to the $k$th coordinate axis are evenly spaced, beginning at the lower bound in the $k$th axis, $l_k$, cutting `S`

into slices of uniform width, except possibly for the last slice. If the $k$ component of `array`

is a vector $C$ of nonnegative exact integers that sum to `S``(interval-width (array-domain `

, then the cuts in the $k$th direction create slices with widths $C_0, C_1, \ldots$, beginning at the lower bound $l_k$. These subarrays completely "tile" `array`) k)

, in the sense that every entry in `array`

is an entry of precisely one entry of the result $T$.`array`

More formally, assume the domain of

is the interval $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$; $T$ is an immutable array with all lower bounds zero. We specify the lower and upper bounds of the array comprising each element of $T$ that is extracted from `array`

in the sense of `array``array-extract`

, as follows.

If the $k$th component of

is an exact positive integer $s$, then the elements of $T$ with $k$th coordinates $j_k$ are subarrays of `S`

with $k$th lower and upper bounds given by $l_k+j_k\times s$ and $\min(l_k+(j_k+1)s, u_k)$, respectively. (The "minimum" operator is necessary if $u_k-l_k$ is not divisible by $s$.)`array`

If, on the other hand, the $k$ component of

is a vector of nonnegative exact integers $C$ whose components sum to $u_k-l_k$, then the elements of $T$ with $k$th coordinates $j_k$ are subarrays of `S`

with $k$th lower and upper bounds given by
$$
l_k+\sum_{i<j_k} C_i\quad\text{ and }\quad l_k+\sum_{i\leq j_k} C_i,\quad\text{respectively.}
$$
`array`

It is an error if the arguments of `array-tile`

do not satisfy these conditions.

If

is a specialized array, the subarrays of the result inherit safety and mutability from `array`

.`array`

**Example: **

```
(define T
(list*->array
2
'(( 1 2 3 4 5 6)
( 7 8 9 10 11 12)
(13 14 15 16 17 18)
(19 20 21 22 23 24)
(25 26 27 28 29 30)
(31 32 33 34 35 36))))
(array->list*
(array-map array->list*
(array-tile T '#(#(3 1 2)
3))))
=>
((((1 2 3) ;; upper left corner
(7 8 9)
(13 14 15))
((4 5 6) ;; upper right corner
(10 11 12)
(16 17 18)))
(((19 20 21)) ;; left middle row
((22 23 24))) ;; right middle row
(((25 26 27) ;; lower left corner
(31 32 33))
((28 29 30) ;; lower right corner
(34 35 36))))
```

**Note: **The procedures `array-tile`

and `array-curry`

both decompose an array into subarrays, but in different ways. For example, if

is defined as `A``(make-array (make-interval '#(10 10)) list)`

, then `(array-tile `

returns an array with domain `A` '#(1 10))`(make-interval '#(10 1))`

for which the value at the multi-index `(`

is an array with domain `i` 0)`(make-interval (vector `

(i.e., a two-dimensional array whose elements are two-dimensional arrays), while `i` 0) (vector (+ `i` 1) 10))`(array-curry `

returns an array with domain `A` 1)`(make-interval '#(10))`

, each element of which has domain `(make-interval '#(10))`

(i.e., a one-dimensional array whose elements are one-dimensional arrays).

**Procedure: **`array-translate `

`array` `translation`

Assumes that

is an array, `array`

is a translation, and that the dimensions of the array and the translation are the same. The resulting array will have domain `translation``(interval-translate (array-domain array) translation)`

.

If

is a specialized array, returns a new specialized array`array`

```
(specialized-array-share
````array`
(interval-translate (array-domain `array`)
`translation`)
(lambda multi-index
(apply values
(map -
multi-index
(vector->list `translation`)))))

that shares the body of

, as well as inheriting its safety and mutability.`array`

If

is not a specialized array but is a mutable array, returns a new mutable array`array`

```
(make-array
(interval-translate (array-domain
````array`)
`translation`)
(lambda multi-index
(apply (array-getter `array`)
(map -
multi-index
(vector->list `translation`))))
(lambda (val . multi-index)
(apply (array-setter `array`)
val
(map -
multi-index
(vector->list `translation`)))))

that employs the same getter and setter as the original array argument.

If

is not a mutable array, returns a new array`array`

```
(make-array
(interval-translate (array-domain
````array`)
`translation`)
(lambda multi-index
(apply (array-getter `array`)
(map - multi-index (vector->list `translation`)))))

that employs the same getter as the original array.

It is an error if the arguments do not satisfy these conditions.

**Example: **

```
(let* ((A (make-array (make-interval '#(2 3)) list))
(B (array-translate A '#(1 -3))))
(display "A:\n")
(array-unveil A)
(interval= (array-domain B) (make-interval '#(1 -3) '#(3 0))) ;; => #t
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
0 0 => (0 0)
0 1 => (0 1)
0 2 => (0 2)
1 0 => (1 0)
1 1 => (1 1)
1 2 => (1 2)
B:
1 -3 => (0 0)
1 -2 => (0 1)
1 -1 => (0 2)
2 -3 => (1 0)
2 -2 => (1 1)
2 -1 => (1 2)
```

**Procedure: **`array-permute `

`array` `permutation`

Assumes that

is an array and `array`

is a permutation, and that the dimensions of the array and the permutation are the same. The resulting array will have domain `permutation``(interval-permute (array-domain array) permutation)`

.

We begin with an example. Assume that the domain of

is represented by the interval $[0,4)\times[0,8)\times[0,21)\times [0,16)$, as in the example for `array``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)$.

So 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

as `old-getter``(array-getter `

, the definition of the new array must be in fact`array`)

`(make-array (interval-permute (array-domain ``array`)
'#(3 0 1 2))
(lambda (l i j k)
(old-getter i j k l)))

So you see that if the first argument if the new getter is in $[0,16)$, then indeed the fourth argument of

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 `old-getter``lambda`

, that must be permuted.

Mathematically, 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 `

$\pi^{-1}$`old-getter` (` multi-index)))`

. We have assumed that $\pi^{-1}$ takes a list as an argument and returns a list as a result.

Employing this same pseudo-code, if

is a specialized array and we denote the permutation by $\pi$, then `array``array-permute`

returns the new specialized array

`(specialized-array-share ``array`
(interval-permute (array-domain `array`) π)
(lambda multi-index
(apply values (π^{-1} multi-index))))

The resulting array shares the body of

, as well as its safety and mutability.`array`

Again employing this same pseudo-code, if

is not a specialized array, but is
a mutable array, then `array``array-permute`

returns the new mutable array

`(make-array (interval-permute (array-domain ``array`) π)
(lambda multi-index
(apply (array-getter `array`)
(π^{-1} multi-index)))
(lambda (val . multi-index)
(apply (array-setter `array`)
val
(π^{-1} multi-index))))

which employs the setter and the getter of the argument to `array-permute`

.

Finally, if

is not a mutable array, then `array``array-permute`

returns

`(make-array (interval-permute (array-domain ``array`) π)
(lambda multi-index
(apply (array-getter `array`)
(π^{-1} multi-index))))

The only length-zero permutation is the empty permutation, specified by `'#()`

.

It is an error to call `array-permute`

if its arguments do not satisfy these conditions.

**Example: **

```
(let* ((A (make-array (make-interval '#(1 3 2)) list))
(B (array-permute A '#(2 1 0))))
(display "A:\n")
(array-unveil A)
(interval= (array-domain B) (make-interval '#(2 3 1))) ;; => #t
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
0 0 0 => (0 0 0)
0 0 1 => (0 0 1)
0 1 0 => (0 1 0)
0 1 1 => (0 1 1)
0 2 0 => (0 2 0)
0 2 1 => (0 2 1)
B:
0 0 0 => (0 0 0)
0 1 0 => (0 1 0)
0 2 0 => (0 2 0)
1 0 0 => (0 0 1)
1 1 0 => (0 1 1)
1 2 0 => (0 2 1)
```

**Procedure: **`array-reverse `

`array` #!optional `flip?`

We assume that

is an array and `array`

, if given, is a vector of booleans whose length is the same as the dimension of `flip?`

. If `array`

is not given, it is set to a vector with length the same as the dimension of `flip?`

, all of whose elements are `array``#t`

.

`array-reverse`

returns a new array that is specialized, mutable, or immutable according to whether

is specialized, mutable, or immutable, respectively. Informally, if `array``(vector-ref `

is true, then the ordering of multi-indices in the k'th coordinate direction is reversed, and is left undisturbed otherwise.`flip?` k)

More formally, we introduce the procedure

```
(define flip-multi-index
(let* ((domain (array-domain
````array`

))
(lowers (interval-lower-bounds->list domain))
(uppers (interval-upper-bounds->list domain)))
(lambda (multi-index)
(map (lambda (i_k flip?_k l_k u_k)
(if flip?_k
(- (+ l_k u_k -1) i_k)
i_k))
multi-index
(vector->list `flip?`)
lowers
uppers))))

Then if

is specialized, `array``array-reverse`

returns

```
(specialized-array-share
````array`

domain
(lambda multi-index
(apply values
(flip-multi-index multi-index))))

and the result inherits the safety and mutability of

.`array`

Otherwise, if

is mutable, then `array``array-reverse`

returns

```
(make-array
domain
(lambda multi-index
(apply (array-getter
````array`

)
(flip-multi-index multi-index)))
(lambda (v . multi-index)
(apply (array-setter `array`

)
v
(flip-multi-index multi-index)))))

Finally, if

is immutable, then `array``array-reverse`

returns

```
(make-array
domain
(lambda multi-index
(apply (array-getter
````array`

)
(flip-multi-index multi-index)))))

It is an error if

and `array`

don't satisfy these requirements.`flip?`

**Example: **The following example was motivated by a blog post by Joe Marshall.

```
(define (palindrome? s)
(let* ((n
(string-length s))
(a
;; an array accessing the characters of s
(make-array (make-interval (vector n))
(lambda (i)
(string-ref s i))))
(ra
;; the characters accessed in reverse order
(array-reverse a))
(half-domain
(make-interval (vector (quotient n 2)))))
;; If n is 0 or 1 the following extracted arrays
;; are empty.
(array-every
char=?
;; the first half of s
(array-extract a half-domain)
;; the reversed second half of s
(array-extract ra half-domain))))
(palindrome? "") => #t
(palindrome? "a") => #t
(palindrome? "aa") => #t
(palindrome? "ab") => #f
(palindrome? "aba") => #t
(palindrome? "abc") => #f
(palindrome? "abba") => #t
(palindrome? "abca") => #f
(palindrome? "abbc") => #f
```

**Procedure: **`array-sample `

`array` `scales`

Assumes that

is an array all of whose lower bounds are zero, and `array`

is a vector of positive exact integers whose length is the same as the dimension of `scales`

.`array`

Informally, if we construct a new matrix $S$ with the entries of

on the main diagonal, then the $\vec i$th element of `scales``(array-sample `

is the $S\vec i$th element of `array` `scales`)

.`array`

More formally, if

is specialized, then `array``array-sample`

returns

```
(specialized-array-share
````array`

(interval-scale (array-domain `array`

)
`scales`

)
(lambda multi-index
(apply values
(map * multi-index (vector->list `scales`

)))))

with the result inheriting the safety and mutability of

.`array`

Otherwise, if

is mutable, then `array``array-sample`

returns

```
(make-array
(interval-scale (array-domain
````array`

)
`scales`

)
(lambda multi-index
(apply (array-getter `array`

)
(map * multi-index (vector->list `scales`

))))
(lambda (v . multi-index)
(apply (array-setter `array`

)
v
(map * multi-index (vector->list `scales`

)))))

Finally, if

is immutable, then `array``array-sample`

returns

```
(make-array
(interval-scale (array-domain
````array`

)
`scales`

)
(lambda multi-index
(apply (array-getter `array`

)
(map * multi-index (vector->list `scales`

)))))

It is an error if

and `array`

don't satisfy these requirements.`scales`

**Example: **

```
(let* ((A (make-array (make-interval '#(3 2)) list))
(B (array-sample A '#(2 1))))
(interval= (array-domain B) (make-interval '#(2 2))) ;; => #t
(display "A:\n")
(array-unveil A)
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
0 0 => (0 0)
0 1 => (0 1)
1 0 => (1 0)
1 1 => (1 1)
2 0 => (2 0)
2 1 => (2 1)
B:
0 0 => (0 0)
0 1 => (0 1)
1 0 => (2 0)
1 1 => (2 1)
```

**Procedure: **`array-outer-product `

`operator` `array1` `array2`

Implements the outer product of

and `array1`

with the operator `array2`

, similar to the APL function with the same name.`operator`

Assumes that

and `array1`

are arrays and that `array2`

is a procedure of two arguments. `operator`

returns the immutable array`array-outer-product`

```
(make-array (interval-cartesian-product (array-domain array1)
(array-domain array2))
(lambda args
(operator (apply (array-getter array1) (take args (array-dimension array1)))
(apply (array-getter array2) (drop args (array-dimension array1))))))
```

This operation can be considered a partial inverse to `array-curry`

. It is an error if the arguments do not satisfy these assumptions.

**Note: **You can see from the above definition that if

is `C``(array-outer-product `

, then each call to `operator` `A` `B`)`(array-getter `

will call `C`)

as well as `operator``(array-getter `

and `A`)`(array-getter `

. This means that if all elements of `B`)

are eventually accessed, then `C``(array-getter `

will be called `A`)`(interval-volume (array-domain `

times; similarly `B`))`(array-getter `

will be called `B`)`(interval-volume (array-domain `

times. `A`))

This implies that if `(array-getter `

is expensive to compute (for example, if it's returning an array, as does `A`)`array-curry`

) then the elements of

should be precomputed if necessary and stored in a specialized array, typically using `A``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 `array-inner-product`

applies `array-outer-product`

to curried arrays, so we apply `array-copy`

to the arguments before passage to `array-outer-product`

.

**Example: **

```
(let* ((A (make-array (make-interval '#(4)) (lambda (i) (* i 10))))
(B (make-array (make-interval '#(3)) values))
(C (array-outer-product + A B)))
(interval= (array-domain C) (make-interval '#(4 3))) ;; => #t
(array-unveil C))
```

displays:

```
0 0 => 0
0 1 => 1
0 2 => 2
1 0 => 10
1 1 => 11
1 2 => 12
2 0 => 20
2 1 => 21
2 2 => 22
3 0 => 30
3 1 => 31
3 2 => 32
```

**Procedure: **`array-inner-product `

`A` `f` `g` `B`

Assumes that

and `f`

are procedures of two arguments and `g`

and `A`

are arrays of dimension at least one, with the upper and lower bounds of the last axis of `B`

the same as those of the first axis of `A`

. Computes the equivalent of`B`

`(define (array-inner-product ``A f g B`)
(array-outer-product
(lambda (`a b`)
(array-reduce `f` (array-map `g a b`)))
(array-copy (array-curry `A` 1))
(array-copy (array-curry (array-permute `B` (index-rotate (array-dimension `B`) 1)))))

We precompute and store the curried arrays using `array-copy`

for efficiency reasons, as described in `array-outer-product`

.

It is an error if the arguments do not satisfy these constraints.

See the extended examples below that use `array-inner-product`

.

**Procedure: **`array-map `

`f` `array` . `arrays`

Assumes that

, `array``(car `

, ... are arrays with the same domain and `arrays`)

is a procedure. Then `f``array-map`

returns a new immutable array with the same domain and getter

```
(lambda multi-index
(apply
````f`
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons `array` `arrays`)))))

It is assumed that

is appropriately defined to be evaluated in this context.`f`

It is expected that `array-map`

and `array-for-each`

will specialize the construction of

```
(lambda multi-index
(apply
````f`
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons `array`
`arrays`)))))

It is an error to call `array-map`

if its arguments do not satisfy these conditions.

**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

by `a`

```
(array-map +
a
(make-array (array-domain a)
(lambda (i) i)))
```

or even

```
(make-array (array-domain a)
(let ((a_ (array-getter a)))
(lambda (i)
(+ (a_ i) i))))
```

**Example: **

```
(let* ((A (make-array (make-interval '#(1 1) '#(5 5)) list))
(B (array-map (lambda (arg) (apply * arg)) A)))
(display "A:\n")
(array-unveil A)
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
1 1 => (1 1)
1 2 => (1 2)
1 3 => (1 3)
1 4 => (1 4)
2 1 => (2 1)
2 2 => (2 2)
2 3 => (2 3)
2 4 => (2 4)
3 1 => (3 1)
3 2 => (3 2)
3 3 => (3 3)
3 4 => (3 4)
4 1 => (4 1)
4 2 => (4 2)
4 3 => (4 3)
4 4 => (4 4)
B:
1 1 => 1
1 2 => 2
1 3 => 3
1 4 => 4
2 1 => 2
2 2 => 4
2 3 => 6
2 4 => 8
3 1 => 3
3 2 => 6
3 3 => 9
3 4 => 12
4 1 => 4
4 2 => 8
4 3 => 12
4 4 => 16
```

**Procedure: **`array-for-each `

`f` `array` . `arrays`

Assumes that

, `array``(car `

, ... are arrays with the same domain and `arrays`)

is a procedure. Then `f``array-for-each`

calls

```
(interval-for-each
(lambda multi-index
(apply
````f`
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons `array`
`arrays`)))))
(array-domain `array`))

In particular, `array-for-each`

always walks the indices of the arrays in lexicographical order.

It is expected that `array-map`

and `array-for-each`

will specialize the construction of

```
(lambda multi-index
(apply
````f`
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons `array`
`arrays`)))))

It is an error to call `array-for-each`

if its arguments do not satisfy these conditions.

**Example: **

```
(let ((A (make-array (make-interval '#(3 3)) list)))
(display "A:\n")
(array-unveil A)
(newline)
(array-for-each (lambda (entry)
(pretty-print (apply + entry)))
A))
```

displays:

```
A:
0 0 => (0 0)
0 1 => (0 1)
0 2 => (0 2)
1 0 => (1 0)
1 1 => (1 1)
1 2 => (1 2)
2 0 => (2 0)
2 1 => (2 1)
2 2 => (2 2)
0
1
2
1
2
3
2
3
4
```

**Procedure: **`array-fold-left `

`operator` `identity` `array` . `arrays`

**Procedure: **`array-fold-right `

`operator` `identity` `array` . `arrays`

These procedures assume that

is a procedure and `operator``(cons `

is a list of arrays all with the same domain.`array arrays`)

These procedures can be defined as:

```
(define (array-fold-left operator identity array . arrays)
(interval-fold-left (array-getter (apply array-map list array arrays))
(lambda (identity array-elements)
(apply operator identity array-elements))
identity
(array-domain array)))
(define (array-fold-right operator identity array . arrays)
(interval-fold-right (array-getter (apply array-map list array arrays))
(lambda (array-elements identity)
(apply operator (append array-elements (list identity))))
identity
(array-domain array)))
```

See the notes for `interval-fold-left`

and `interval-fold-right`

.

It is an error if the arguments do not satisfy these assumptions.

**Note: ** One can fold over empty arrays, which returns

, but it is an error to call `identity``array-reduce`

on an empty array, because `array-reduce`

must evaluate at least one element of the argument array.

**Example: ** One can define an APL-style `array-depth`

by:

```
(define (array-depth a)
(if (array? a)
(+ 1 (array-fold-left max 0 (array-map array-depth a)))
0))
```

Here non-arrays have depth 0, and each level of array "nesting" increases the depth by 1.

**Example: ** One can define `(array-foldl-on op id array dims)`

, which does not fold over the entire array, but computes:

```
(define (array-foldl-on op id array dims)
(array-map (lambda (a)
(array-fold-left op id a))
(array-curry a dims)))
```

which folds over only the `dims`

rightmost dimensions and returns an array of results. (Note that this works even if `dims`

is `(array-dimension array)`

, in which case the result is a zero-dimensional array containing the left fold of the entire array.)

**Example: **If

is associative with two-sided identity `op`

, then `id``array-fold-left`

and `array-fold-right`

return the same results, but see:

```
(define a (make-array (make-interval '#(10)) (lambda (i) i)))
(array-fold-left cons '() a)
=> ((((((((((() . 0) . 1) . 2) . 3) . 4) . 5) . 6) . 7) . 8) . 9)
(array-fold-right cons '() a)
=> (0 1 2 3 4 5 6 7 8 9)
(array-fold-left - 0 a)
=> -45
(array-fold-right - 0 a)
=> -5
```

**Procedure: **`array-reduce `

`operator` `array`

Assumes that

is a nonempty array and `array`

is a procedure of two arguments that is associative, i.e., `operator``(`

is the same as `operator` (`operator` `x` `y`) `z`)`(`

.`operator` `x` (`operator` `y` `z`))

Then `(array-reduce `

can be defined as`operator` `array`)

```
(define array-reduce
(let ((reduce-base (list 1))) ;; any unique object
(lambda (sum array)
(array-fold-left (lambda (id entry)
(if (eq? id reduce-base)
entry
(sum id entry)))
reduce-base
array))))
```

The implementation is allowed to use the associativity of

to reorder the computations in `operator``array-reduce`

. It is an error if the arguments do not satisfy these conditions.

**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.

We find with $m=1{,}000{,}000{,}000$:

```
(define A (make-array (make-interval '#(1) '#(1000000001))
(lambda (k)
(fl/ (flsquare (inexact k))))))
(define (block-sum A)
(let ((N (interval-volume (array-domain A))))
(cond ((<= N 1000)
(array-reduce fl+ A))
((<= N (square 1000))
(block-sum (array-map block-sum
(array-tile A (vector (integer-sqrt N))))))
(else
(block-sum (array-map block-sum
(array-tile A (vector (quotient N 1000)))))))))
(array-reduce fl+ A) => 1.644934057834575
(block-sum A) => 1.6449340658482325
```

Since $\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.

**Procedure: **`array-any `

`predicate` `array` . `arrays`

Assumes that `(cons `

is a list of arrays, all with the same domain, which we'll call `array arrays`)`interval`

. Also assumes that

is a procedure that takes as many arguments as there are arrays and returns a single value.`predicate`

`array-any`

first computes `(apply `

to the first element of `predicate` (map (lambda (g_) (apply g_ multi-index)) (map array-getter (cons `array arrays`))))`interval`

in lexicographical order.

If the result of

is not `predicate``#f`

, then that result is returned by `array-any`

. If the result of

is `predicate``#f`

, then `array-any`

continues with the second element of `interval`

, etc., returning the first nonfalse value of

.`predicate`

If

always returns `predicate``#f`

, then `array-any`

returns `#f`

.

If it happens that

is applied to `predicate``(map (lambda (g_) (apply g_ multi-index)) (map array-getter (cons `

with `array arrays`)))`multi-index`

the last element of `interval`

, then this last call to

is in tail position.`predicate`

The procedures `(array-getter `

, etc., are applied only to those values of `array`)`interval`

necessary to determine the result of `array-any`

.

It is an error if the arguments do not satisfy these assumptions.

**Example: **

```
(let ((A (make-array (make-interval '#(240) '#(250)) values))
(B (make-array (make-interval '#(250) '#(300)) values)))
(define (square? n)
(and (exact? (sqrt n)) n)) ;; return the value
(array-any square? A) ;; => #f
(array-any square? B)) ;; => 256, the first nonfalse value
```

**Procedure: **`array-every `

`predicate` `array` . `arrays`

Assumes that`(cons `

is a list arrays, all with the same domain, which we'll call `array arrays`)`interval`

. Also assumes that

is a procedure that takes as many arguments as there are arrays and returns a single value.`predicate`

`array-every`

first computes `(apply `

to the first element of `predicate` (map (lambda (g_) (apply g_ multi-index)) (map array-getter (cons `array arrays`))))`interval`

in lexicographical order.

If the result of

is `predicate``#f`

, then that result is returned by `array-every`

. If the result of

is nonfalse, then `predicate``array-every`

continues with the second element of `interval`

, etc., returning the first value of

that is `predicate``#f`

.

If

always returns a nonfalse value, then the last nonfalse value returned by `predicate`

is also returned by `predicate``array-every`

.

If it happens that

is applied to `predicate``(map (lambda (g_) (apply g_ multi-index)) (map array-getter (cons `

with `array arrays`)))`multi-index`

the last element of `interval`

, then this last call to

is in tail position.`predicate`

The procedures `(array-getter `

, etc., are applied only to those values of `array`)`interval`

necessary to determine the result of `array-every`

.

It is an error if the arguments do not satisfy these assumptions.

For an example, see the palindrome example above.

**Procedure: **`array->list `

`array`

Stores the elements of

into a newly allocated list in lexicographical order. It is an error if `array`

is not an array.`array`

It is guaranteed that `(array-getter `

is called precisely once for each multi-index in `array`)`(array-domain `

in lexicographical order.`array`)

**Example: **

```
(let* ((A (make-specialized-array-from-data '#(2 4 6 8)))
(B (array-reverse A)))
(pretty-print (array->list A))
(pretty-print (array->list B)))
```

displays:

```
(2 4 6 8)
(8 6 4 2)
```

**Procedure: **`list->array `

`interval` `list` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a list, `list`

is an interval with volume the same as the length of `interval`

, `list`

is a storage class that can manipulate all the elements of `storage-class`

, and `list`

and `mutable?`

are booleans.`safe?`

Returns a specialized array with domain

whose elements are the elements of the list `interval`

stored in lexicographical order. The result is mutable or safe depending on the values of `list`

and `mutable?`

.`safe?`

Any missing optional arguments are assigned the values `generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if the arguments do not satisfy these assumptions, or if any element of

cannot be stored in the body of `list`

, and this last error shall be detected and raised.`storage-class`

**Example: **

```
(let* ((l (iota 12))
(A (list->array (make-interval '#(2 2 3)) l))
(B (list->array (make-interval '#(12)) l)))
(display "A:\n")
(array-unveil A)
(display "B:\n")
(array-unveil B))
```

displays

```
A:
0 0 0 => 0
0 0 1 => 1
0 0 2 => 2
0 1 0 => 3
0 1 1 => 4
0 1 2 => 5
1 0 0 => 6
1 0 1 => 7
1 0 2 => 8
1 1 0 => 9
1 1 1 => 10
1 1 2 => 11
B:
0 => 0
1 => 1
2 => 2
3 => 3
4 => 4
5 => 5
6 => 6
7 => 7
8 => 8
9 => 9
10 => 10
11 => 11
```

**Procedure: **`list*->array `

`d` `nested-list` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a nonnegative exact integer and, if given, `d`

is a storage class and `storage-class`

and `mutable?`

are booleans.`safe?`

This routine builds a specialized array of dimension

, storage class `d`

, mutability `storage-class`

, and safety `mutable?`

from `safe?`

. It is assumed that following predicate does not return `nested-list``#f`

when passed

and `nested-list`

as arguments:`d`

```
(define (check-nested-list dimension nested-data)
(or (eqv? dimension 0) ;; anything goes in dimension 0
(and (list? nested-data)
(let ((len (length nested-data)))
(cond ((eqv? len 0)
'())
((eqv? dimension 1)
(list len))
(else
(let* ((sublists
(map (lambda (l)
(check-nested-list (fx- dimension 1) l))
nested-data))
(first
(car sublists)))
(and first
(every (lambda (l)
(equal? first l))
(cdr sublists))
(cons len first)))))))))
```

In this case, `list*->array`

returns an array with domain `(make-interval (list->vector (check-nested-list `

. If we denote the getter of the result by `d nested-list`)))`A_`

, then

```
(A_ i_0 ... i_d-2 i_d-1)
=> (list-ref (list-ref (... (list-ref nested-list i_0) ...) i_d-2) i_d-1)
```

and we assume that this value can be manipulated by

.`storage-class`

Any missing optional arguments are assigned the values `generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

Empty and zero-dimensional lists are treated differently; see the discussion for `array->list*`

. For example

```
(list*->array 0 '()) => An array for which ((array-getter (list*->array 0 '()))) => '()
(list*->array 1 '()) => An empty array with domain (make-interval '#(0))
(list*->array 2 '()) => An empty array with domain (make-interval '#(0 0))
(list*->array 2 '(() ())) => An empty array with domain (make-interval '#(2 0))
```

It is an error if the arguments do not satisfy these assumptions.

**Example: **

```
(let ((A (list*->array 3 '(((1 2 3)
(4 5 6))
((7 8 9)
(10 11 12))))))
(array-unveil A))
```

displays:

```
0 0 0 => 1
0 0 1 => 2
0 0 2 => 3
0 1 0 => 4
0 1 1 => 5
0 1 2 => 6
1 0 0 => 7
1 0 1 => 8
1 0 2 => 9
1 1 0 => 10
1 1 1 => 11
1 1 2 => 12
```

**Procedure: **`array->list* `

`array`

Assumes that

is an array, and returns a newly allocated nested list `array`

. If `nested-list`

is nonempty and has positive dimension and we denote the getter of `array`

by `array``array_`

, then

and `nested-list``array_`

satisfy

```
(array_ i_0 ... i_d-2 i_d-1)
=> (list-ref (list-ref (... (list-ref nested-list i_0) ...) i_d-2) i_d-1)
```

Each element of

is accessed once.`array`

If

is zero dimensional, then `array``array->list*`

returns `((array-getter `

. If the argument is an empty array, then the nested lists of the result match the first nonzero dimensions (if any). For example:`array`))

```
(array->list* (make-array (make-interval '#()) (lambda () 2))) => 2 ;; no list
(array->list* (make-array (make-interval '#(0)) error)) => '()
(array->list* (make-array (make-interval '#(0 0)) error)) => '()
(array->list* (make-array (make-interval '#(2 0)) error)) => '(() ())
(array->list* (make-array (make-interval '#(0 2)) error)) => '()
```

It is an error if

is not an array.`array`

**Example: **

```
(let ((B (array->list* (make-array (make-interval '#(6 6)) (lambda (i j) (/ (+ 1 i j)))))))
(pretty-print B))
```

displays:

```
((1 1/2 1/3 1/4 1/5 1/6)
(1/2 1/3 1/4 1/5 1/6 1/7)
(1/3 1/4 1/5 1/6 1/7 1/8)
(1/4 1/5 1/6 1/7 1/8 1/9)
(1/5 1/6 1/7 1/8 1/9 1/10)
(1/6 1/7 1/8 1/9 1/10 1/11))
```

**Procedure: **`array->vector `

`array`

Stores the elements of

into a newly allocated vector in lexicographical order. It is an error if `array`

is not an array.`array`

It is guaranteed that `(array-getter `

is called precisely once for each multi-index in `array`)`(array-domain `

in lexicographical order.`array`)

**Example: **

```
(let* ((A (make-specialized-array-from-data '#(2 4 6 8)))
(B (array-reverse A)))
(pretty-print (array->vector A))
(pretty-print (array->vector B)))
```

displays:

```
#(2 4 6 8)
#(8 6 4 2)
```

**Procedure: **`vector->array `

`interval` `vector` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a vector, `vector`

is an interval with volume the same as the length of `interval`

, `v`

is a storage class that can manipulate all the elements of `storage-class`

, and `vector`

and `mutable?`

are booleans.`safe?`

Returns a specialized array with domain

whose elements are the elements of the vector `interval`

stored in lexicographical order. The result is mutable or safe depending on the values of `vector`

and `mutable?`

.`safe?`

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if the arguments do not satisfy these assumptions, or if any element of

cannot be stored in the body of `vector`

, and this last error shall be detected and raised.`storage-class`

**Example: **

```
(let* ((v (list->vector (iota 12)))
(A (vector->array (make-interval '#(2 2 3)) v))
(B (vector->array (make-interval '#(12)) v)))
(display "A:\n")
(array-unveil A)
(display "B:\n")
(array-unveil B))
```

displays:

```
A:
0 0 0 => 0
0 0 1 => 1
0 0 2 => 2
0 1 0 => 3
0 1 1 => 4
0 1 2 => 5
1 0 0 => 6
1 0 1 => 7
1 0 2 => 8
1 1 0 => 9
1 1 1 => 10
1 1 2 => 11
B:
0 => 0
1 => 1
2 => 2
3 => 3
4 => 4
5 => 5
6 => 6
7 => 7
8 => 8
9 => 9
10 => 10
11 => 11
```

**Procedure: **`vector*->array `

`d` `nested-vector` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a nonnegative exact integer and, if given, `d`

is a storage class and `storage-class`

and `mutable?`

are booleans.`safe?`

This routine builds a specialized array of dimension

, storage class `d`

, mutability `storage-class`

, and safety `mutable?`

from `safe?`

. It is assumed that following predicate does not return `nested-vector``#f`

when passed

and `nested-vector`

as arguments:`d`

```
(define (check-nested-vector dimension nested-data)
(or (eqv? dimension 0) ;; anything goes in dimension 0
(and (vector? nested-data)
(let ((len (vector-length nested-data)))
(cond ((eqv? len 0)
'())
((eqv? dimension 1)
(list len))
(else
(let* ((sublists
(vector-map (lambda (l)
(check-nested-vector (fx- dimension 1) l))
nested-data))
(first
(vector-ref sublists 0)))
(and first
(vector-every (lambda (l)
(equal? first l))
sublists)
(cons len first)))))))))
```

In this case, `vector*->array`

returns an array with domain `(make-interval (list->vector (check-nested-vector `

. If we denote the getter of the result by `d nested-vector`)))`A_`

, then

```
(A_ i_0 ... i_d-2 i_d-1)
=> (vector-ref (vector-ref (... (vector-ref nested-vector i_0) ...) i_d-2) i_d-1)
```

and we assume that this value can be manipulated by

.`storage-class`

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

If the resulting array would be empty or have dimension zero, see the examples for `list*->array`

.

It is an error if the arguments do not satisfy these assumptions.

**Example: **

```
(let ((A (vector*->array 3 '#(#(#(1 2 3)
#(4 5 6))
#(#(7 8 9)
#(10 11 12))))))
(array-unveil A))
```

displays:

```
0 0 0 => 1
0 0 1 => 2
0 0 2 => 3
0 1 0 => 4
0 1 1 => 5
0 1 2 => 6
1 0 0 => 7
1 0 1 => 8
1 0 2 => 9
1 1 0 => 10
1 1 1 => 11
1 1 2 => 12
```

**Procedure: **`array->vector* `

`array`

Assumes that

is an array, and returns a newly allocated nested vector `array`

. If we denote the getter of `nested-vector`

by `array``array_`

, then

and `nested-vector``array_`

satisfy

```
(array_ i_0 ... i_d-2 i_d-1)
=> (vector-ref (vector-ref (... (vector-ref nested-vector i_0) ...) i_d-2) i_d-1)
```

If

is empty or zero dimensional, then see the examples for `array``array->list*`

.:

Each element of

is accessed once.`array`

It is an error if

is not an array.`array`

**Example: **

```
(let ((B (array->vector* (make-array (make-interval '#(6 6)) (lambda (i j) (/ (+ 1 i j)))))))
(pretty-print B))
```

displays:

```
#(#(1 1/2 1/3 1/4 1/5 1/6)
#(1/2 1/3 1/4 1/5 1/6 1/7)
#(1/3 1/4 1/5 1/6 1/7 1/8)
#(1/4 1/5 1/6 1/7 1/8 1/9)
#(1/5 1/6 1/7 1/8 1/9 1/10)
#(1/6 1/7 1/8 1/9 1/10 1/11))
```

**Procedure: **`array-assign! `

`destination` `source`

Assumes that

is a mutable array and `destination`

is an array with the same domain, and that the elements of `source`

can be stored into `source`

.`destination`

Evaluates `(array-getter `

on the multi-indices in `source`)`(array-domain `

in lexicographical order, and assigns each value to the multi-index in `source`)

in the same lexicographical order.`destination`

It is an error if the arguments don't satisfy these assumptions.

If assigning any element of

affects the value of any element of `destination`

, then the result is undefined.`source`

**Example: **

```
(let* ((A (array-copy
(make-array (make-interval '#(5 5))
(lambda (i j) (* i j)))
generic-storage-class
#t)) ;; ensure A is mutable
(D_B (make-interval '#(2 2) '#(5 5)))
(B (make-array D_B (lambda (i j) 100))))
(display "A before assignment:\n")
(pretty-print (array->list* A))
(array-assign! (array-extract A D_B)
B)
(display "A after assignment:\n")
(pretty-print (array->list* A)))
```

displays:

```
A before assignment:
((0 0 0 0 0)
(0 1 2 3 4)
(0 2 4 6 8)
(0 3 6 9 12)
(0 4 8 12 16))
A after assignment:
((0 0 0 0 0)
(0 1 2 3 4)
(0 2 100 100 100)
(0 3 100 100 100)
(0 4 100 100 100))
```

**Procedure: **`array-stack `

`k` `arrays` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

**Procedure: **`array-stack! `

`k` `arrays` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a nonempty list of arrays with identical domains, `arrays`

is an exact integer between 0 (inclusive) and the dimension of the array domains (inclusive), and, if given, `k`

is a storage class, `storage-class`

is a boolean, and `mutable?`

is a boolean.`safe?`

Returns a specialized array equivalent to

```
(array-copy
(make-array
(let ((
````lowers` (interval-lower-bounds->list (array-domain (car `arrays`))))
(`uppers` (interval-upper-bounds->list (array-domain (car `arrays`))))
(`N` (length `arrays`)))
(make-interval (list->vector (append (take `lowers k`) (cons 0 (drop `lowers k`))))
(list->vector (append (take `uppers k`) (cons `N` (drop `uppers k`))))))
(let ((`getters` (map array-getter `arrays`)))
(lambda indices
(let ((`i` (list-ref `indices k`)))
(apply (list-ref `getters i`)
(append (take `indices k`)
(drop `indices` (+ `k` 1)))))))))

In other words we "stack" the argument arrays along a new

'th axis, the lower bound of which is set to 0.`k`

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

Each element of any of the

is accessed once.`arrays`

It is an error if the arguments do not satisfy these constraints.

**Example: **Let's say we have a spreadsheet

and we want to make a new spreadsheet `A`

with the same rows but with the data from only columns 1, 2, 5, and 8. Using the routine `B``array-display`

we define below, code to do this can look like:

```
(let* ((A
(make-array
(make-interval '#(4 10))
list))
(column_
(array-getter ;; the getter of ...
(array-curry ;; a 1-D array of the columns of A
(array-permute A '#(1 0))
1)))
(B
(array-stack ;; stack into a new 2-D array ...
1 ;; along axis 1 (i.e., columns) ...
(map column_ '(1 2 5 8))))) ;; the columns of A you want
(array-display B))
;;; Displays
(0 1) (0 2) (0 5) (0 8)
(1 1) (1 2) (1 5) (1 8)
(2 1) (2 2) (2 5) (2 8)
(3 1) (3 2) (3 5) (3 8)
```

In fact, because

is a generalized array, the only elements of `A`

that are generated are the ones that are assigned as elements of `A`

. The result could also be computed in one line:`B`

`(array-stack 1 (map (array-getter (array-curry (array-permute A '#(1 0)) 1)) '(1 2 5 8)))`

**Procedure: **`array-decurry `

`AofA` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

**Procedure: **`array-decurry! `

`AofA` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a nonempty array of arrays; the elements of `AofA`

are assumed to all have the same (possibly empty) domain. Also assumes that, if given, `AofA`

is a storage class and `storage-class`

and `mutable?`

are booleans.`safe?`

`array-decurry`

evaluates each array element of

once, and evaluates each element of `AofA``AofA`

's array elements once. `array-decurry`

returns a specialized array containing the elements of

's array elements; ignoring optional arguments, the result is equivalent to: `AofA`

```
(let* ((A
(array-copy AofA)) ;; evaluate all elements of A once
(A_dim
(array-dimension A))
(A_
(array-getter A))
(A_D
(array-domain A))
(element-domain
(array-domain (apply A_ (interval-lower-bounds->list A_D))))
(result-domain
(interval-cartesian-product A_D element-domain))
(result
(make-specialized-array result-domain))
(curried-result
(array-curry result (interval-dimension element-domain))))
(array-for-each array-assign! curried-result A)
result-array)
```

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if any of these assumptions are not met, or if the given storage class cannot manipulate the elements of

's array elements.`AofA`

**Example: **

```
(let* ((A (list*->array 1 '(1 2 3)))
(B (list*->array 1 '(4 5 6)))
(C (list*->array 1 '(7 8 9)))
(D (list*->array 1 '(10 11 12)))
(E (list*->array 1 (list A B C D)))
(F (array-decurry E)))
(array-every array? E) ;; => #t
(array-every
(lambda (a)
(interval=
(array-domain a)
(make-interval '#(3)))) ;; => #t
E)
(interval= (array-domain F)
(make-interval '#(4 3))) ;; => #t
(array-unveil F))
```

displays:

```
0 0 => 1
0 1 => 2
0 2 => 3
1 0 => 4
1 1 => 5
1 2 => 6
2 0 => 7
2 1 => 8
2 2 => 9
3 0 => 10
3 1 => 11
3 2 => 12
```

**Procedure: **`array-append `

`k` `arrays` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

**Procedure: **`array-append! `

`k` `arrays` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

Assumes that

is a nonempty list of arrays with domains that differ at most in the `arrays`

'th axis, `k`

is an exact integer between 0 (inclusive) and the dimension of the array domains (exclusive), and, if given, `k`

is a storage class, `storage-class`

is a boolean, and `mutable?`

is a boolean.`safe?`

This routine appends, or concatenates, the argument arrays along the `k`'th axis, with the lower bound of this axis set to 0.

Returns a specialized array equivalent to the result of

```
(define (array-append k arrays)
(let*-values (((axis-subdividers kth-size)
;; compute lower and upper bounds of where along the
;; k'th axis we'll copy each array argument, plus
;; the total size of the kth axis of the result array
(let loop ((result '(0))
(arrays arrays))
(if (null? arrays)
(values (reverse result) (car result))
(let ((interval (array-domain (car arrays))))
(loop (cons (+ (car result)
(- (interval-upper-bound interval k)
(interval-lower-bound interval k)))
result)
(cdr arrays))))))
((lowers)
;; the domains of the arrays differ only in the kth axis
(interval-lower-bounds->vector (array-domain array)))
((uppers)
(interval-upper-bounds->vector (array-domain array)))
((result)
;; the result array
(make-specialized-array
(let ()
(vector-set! lowers k 0)
(vector-set! uppers k kth-size)
(make-interval lowers uppers))))
((translation)
;; a vector we'll use to align each argument
;; array into the proper subarray of the result
(make-vector (array-dimension array) 0)))
(let loop ((arrays arrays)
(subdividers axis-subdividers))
(if (null? arrays)
;; we've assigned every array to the appropriate subarray of result
result
(let ((array (car arrays)))
;; the lower and upper bounds in the kth axis of the result where we copy the
;; next array
(vector-set! lowers k (car subdividers))
(vector-set! uppers k (cadr subdividers))
;; the translation that aligns the next array with the subarray of the result
(vector-set! translation k (- (car subdividers)
(interval-lower-bound (array-domain array) k)))
(array-assign!
(array-extract result (make-interval lowers uppers))
(array-translate array translation))
(loop (cdr arrays)
(cdr subdividers)))))))
```

Each element of any of the

is accessed once.`arrays`

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if the arguments do not satisfy these constraints.

**Example:** Given a two-dimensional array $a$ interpreted as a spreadsheet, with the rows and columns indexed starting at 0, one might want to make a new array with row $k$ moved to be the top row. Then one could do:

```
(let* ((a (make-array (make-interval '#(4 6)) list))
(k 2)
(m (interval-upper-bound (array-domain a) 0))
(n (interval-upper-bound (array-domain a) 1)))
(pretty-print
(array->list* a))
(newline)
(pretty-print
(array->list*
(array-append
0
(list (array-extract a (make-interval (vector k 0) (vector (+ k 1) n)))
(array-extract a (make-interval (vector k n)))
(array-extract a (make-interval (vector (+ k 1) 0) (vector m n))))))))
```

This prints:

```
(((0 0) (0 1) (0 2) (0 3) (0 4) (0 5))
((1 0) (1 1) (1 2) (1 3) (1 4) (1 5))
((2 0) (2 1) (2 2) (2 3) (2 4) (2 5))
((3 0) (3 1) (3 2) (3 3) (3 4) (3 5)))
(((2 0) (2 1) (2 2) (2 3) (2 4) (2 5))
((0 0) (0 1) (0 2) (0 3) (0 4) (0 5))
((1 0) (1 1) (1 2) (1 3) (1 4) (1 5))
((3 0) (3 1) (3 2) (3 3) (3 4) (3 5)))
```

Because this SRFI supports empty arrays, the same code works when $k=0$ (when the second extracted array is empty) or $k=m-1$ (when the third extracted array is empty).

**Procedure: **`array-block `

`AofA` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

**Procedure: **`array-block! `

`AofA` [ `storage-class` [ `mutable?` [ `safe?` ] ] ]

This procedure is an inverse to `array-tile`

. It assumes that

is a nonempty array of arrays, all of which have the same dimension as `AofA`

itself. It also assumes that, if given, `AofA`

is a storage class and `storage-class`

and `mutable?`

are booleans.`safe?`

While ignoring the lower and upper bounds of the element arrays, it assumes that those element arrays have widths (as defined by `interval-widths`

) that allow them to be packed together in the configuration given by their indices in

. We can always do this when `AofA``(array-dimension `

is 1. Otherwise, assuming that the lower bounds of `AofA`)

are zero, we require: `AofA`

```
(every
(lambda (k) ;; for each coordinate direction
(let ((slices ;; the "slices"
;; perpendicular to that direction
(array-curry (array-permute AofA (index-first (array-dimension AofA) k))
(- (array-dimension AofA) 1))))
(array-every
(lambda (slice) ;; for every slice perpendicular
;; to direction k
(let ((slice-kth-width ;; the kth interval width of the
;; "corner" element
(interval-width
(array-domain
(apply (array-getter slice)
(make-list (- (array-dimension AofA) 1) 0)))
k)))
(array-every
(lambda (a) ;; all arrays within that slice
(= (interval-width (array-domain a) k) ;; have the same width in the kth
;; direction
slice-kth-width))
slice)))
slices)))
(iota (array-dimension AofA)))
```

This procedure then returns a specialized array with lower bounds all zero and with the specified storage class, mutability, and safety, whose elements are taken from the array elements of

itself. In principle, one could compute the result by appending all the array elements of `AofA`

successively along each coordinate axis of `AofA`

, in any order of the axes. Each element of `AofA`

is accessed once, and each element of `AofA`

's array elements is accessed once.`AofA`

Each element of

is itself an array; one can copy the contents of each array element of `AofA`

to the result array with the following algorithm:`AofA`

```
(let* ((A_dim
(array-dimension AofA))
(ks
(list->vector (iota A_dim)))
(corner-multi-index
(make-list (fx- A_dim 1) 0))
(slice-offsets ;; the indices in each direction where the "cuts" are
(vector-map
(lambda (k) ;; the direction
(let* ((pencil ;; a pencil in that direction
(apply (array-getter (array-curry (array-permute AofA (index-last A_dim k)) 1))
corner-multi-index))
(pencil_
(array-getter pencil))
(pencil-size
(interval-width (array-domain pencil) 0))
(result ;; include sum of all kth interval-widths in pencil
(make-vector (fx+ pencil-size 1) 0)))
(do ((i 0 (fx+ i 1)))
((fx= i pencil-size) result)
(vector-set! result
(fx+ i 1)
(fx+ (vector-ref result i)
(interval-width (array-domain (pencil_ i)) k))))))
ks))
(result
(make-specialized-array
(make-interval
(vector-map (lambda (v)
(vector-ref v (fx- (vector-length v) 1)))
slice-offsets))
storage-class
(storage-class-default storage-class)
safe?)))
;; We copy the elements from each input array block to the corresponding block
;; in the result array.
(interval-for-each
(lambda multi-index
(let* ((vector-multi-index
(list->vector multi-index))
(corner ;; where the subarray will sit in the result array
(vector-map (lambda (i k)
(vector-ref (vector-ref slice-offsets k) i))
vector-multi-index
ks))
(subarray
(apply array-ref AofA multi-index))
(translated-subarray ;; translate the subarray to corner
(array-translate
subarray
(vector-map -
corner
(interval-lower-bounds (array-domain subarray))))))
(array-assign! (array-extract result (array-domain translated-subarray))
translated-subarray)))
(array-domain AofA))
(if (not mutable?)
(array-freeze! result)
result))
```

`generic-storage-class`

, `(specialized-array-default-mutable?)`

, and `(specialized-array-default-safe?)`

, respectively.

It is an error if the arguments do not satisfy these assumptions, or if all elements of the result cannot by manipulated by the given storage class.

**Examples: **

```
(array->vector*
(array-block (list*->array
2
(list (list (list*->array 2 '((0 1)
(2 3)))
(list*->array 2 '((4)
(5)))
(list*->array 2 '((6 7 8)
(9 10 11))))
(list (list*->array 2 '((12 13)))
(list*->array 2 '((14)))
(list*->array 2 '((15 16 17))))))))
=>
#(#(0 1 4 6 7 8)
#(2 3 5 9 10 11)
#(12 13 14 15 16 17))
(array-block (list*->array
2
(list (list (list*->array 2 '((0 1)
(2 3)))
(list*->array 2 '((4)
(5)))
(list*->array 2 '((6 7) ;; these should each have ...
(9 10)))) ;; three elements ...
(list (list*->array 2 '((12 13)))
(list*->array 2 '((14)))
(list*->array 2 '((15 16 17))))))) ;; to match this array
=> error
```

**Procedure: **`array-ref `

`array` . `multi-index`

Assumes that

is an array, and `array`

is a sequence of exact integers.`multi-index`

Returns `(apply (array-getter `

.`array`) `multi-index`)

It is an error if

is not an array, if the number of elements in `array`

is not the the dimension of `multi-index`

, or if `array`

is not in the domain of `multi-index`

, so, in particular, if `array`

is empty.`array`

**Example: **

```
(let ((A (make-array (make-interval '#(10000 10000)) expt)))
(array-ref A 5 37) ;; => 72759576141834259033203125
(array-ref A 37 5)) ;; => 69343957
```

**Procedure: **`array-set! `

`array` `object` . `multi-index`

Assumes that

is a mutable array, that `array`

is a value that can be stored within that array, and that `object`

is a sequence of exact integers.`multi-index`

Returns `(apply (array-setter `

.`array`) `object multi-index`)

It is an error if

is not a mutable array, if `array`

is not an appropriate value to be stored in that array, if the number of elements in `object`

is not the the dimension of `multi-index`

, or if `array`

is not in the domain of `multi-index`

, so, in particular, if `array`

is empty.`array`

**Example: **

```
(let ((A (array-copy
(list*->array 1 (iota 1000))
generic-storage-class
#t))) ;; ensure array is mutable
(array-ref A 500) ;; => 500
(array-set! A 'grok 500)
(array-ref A 500)) ;; => grok
```

**Note: **In the sample implementation, because `array-ref`

and `array-set!`

take a variable number of arguments and they must check that

is an array of the appropriate type, programs written in a style using these procedures, rather than the style in which `A``1D-Haar-loop`

is coded below, can take up to three times as long runtime.

**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.

**Procedure: **`specialized-array-reshape `

`array` `interval` [ `copy-on-failure?` #f ]

Assumes that

is a specialized array, `array`

is an interval with the same volume as `interval``(array-domain `

, and `array`)

, if given, is a boolean.`copy-on-failure?`

If there is an affine map that takes the multi-indices in

to the cells in `interval``(array-body `

storing the elements of `array`)

in lexicographical order, `array``specialized-array-reshape`

returns a new specialized array, with the same body and elements as

and domain `array`

. The result inherits its mutability and safety from `interval`

.`array`

If there is not an affine map that takes the multi-indices in

to the cells storing the elements of `interval`

in lexicographical order and `array`

is `copy-on-failure?``#t`

, then returns a specialized array copy of

with domain `array`

, storage class `interval``(array-storage-class `

, mutability `array`)`(mutable-array? `

, and safety `array`)`(array-safe? `

.`array`)

It is an error if these conditions on the arguments are not met.

**Note: **The code in the sample implementation to determine whether there exists an affine map from

to the multi-indices of the elements of `interval`

in lexicographical order is modeled on the corresponding code in the Python library NumPy.`array`

**Note: **In the sample implementation, if an array cannot be reshaped and

is `copy-on-failure?``#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.

**Examples: **Reshaping an array is not a Bawden-type array transform. For example, we use `array-display`

defined below to see:

```
;;; The entries of A are the multi-indices of the locations
(define A (array-copy (make-array (make-interval '#(3 4)) list)))
(array-display A)
;;; Displays
;;; (0 0) (0 1) (0 2) (0 3)
;;; (1 0) (1 1) (1 2) (1 3)
;;; (2 0) (2 1) (2 2) (2 3)
(array-display (array-permute A '#(1 0)))
;;; Displays
;;; (0 0) (1 0) (2 0)
;;; (0 1) (1 1) (2 1)
;;; (0 2) (1 2) (2 2)
;;; (0 3) (1 3) (2 3)
(array-display (specialized-array-reshape A (make-interval '#(4 3))))
;;; Displays
;;; (0 0) (0 1) (0 2)
;;; (0 3) (1 0) (1 1)
;;; (1 2) (1 3) (2 0)
;;; (2 1) (2 2) (2 3)
(define B (array-sample A '#(2 1)))
(array-display B)
;;; Displays
;;; (0 0) (0 1) (0 2) (0 3)
;;; (2 0) (2 1) (2 2) (2 3)
(array-display (specialized-array-reshape B (make-interval '#(8)))) => fails
(array-display (specialized-array-reshape B (make-interval '#(8)) #t))
;;; Displays
;;; (0 0) (0 1) (0 2) (0 3) (2 0) (2 1) (2 2) (2 3)
```

The following examples succeed:

```
(specialized-array-reshape
(array-copy (make-array (make-interval '#(2 1 3 1)) list))
(make-interval '#(6)))
(specialized-array-reshape
(array-copy (make-array (make-interval '#(2 1 3 1)) list))
(make-interval '#(3 2)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
(make-interval '#(6)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
(make-interval '#(3 2)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
(make-interval '#(3 2)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
(make-interval '#(3 1 2 1)))
(specialized-array-reshape
(array-sample
(array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list))
'#(#f #f #f #t))
'#(1 1 2 1))
(make-interval '#(4)))
(specialized-array-reshape
(array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list))
'#(#t #f #t #t))
'#(1 1 2 1))
(make-interval '#(4)))
```

The following examples raise an exception:

```
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
(make-interval '#(6)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
(make-interval '#(3 2)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #f))
(make-interval '#(6)))
(specialized-array-reshape
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #t))
(make-interval '#(3 2)))
(specialized-array-reshape
(array-sample
(array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list))
'#(#f #f #f #t))
'#(1 1 2 1))
(make-interval '#(4)) )
(specialized-array-reshape
(array-sample
(array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list))
'#(#f #f #t #t))
'#(1 1 2 1))
(make-interval '#(4)))
```

In 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.

```
(define interval-flat (make-interval '#(100 100 4)))
(define interval-2x2 (make-interval '#(100 100 2 2)))
(define A (array-copy (make-array interval-flat (lambda args (random-integer 5)))))
(define B (array-copy (make-array interval-flat (lambda args (random-integer 5)))))
(define C (array-copy (make-array interval-flat (lambda args 0))))
(define (2x2-matrix-multiply-into! A B C)
(let ((C! (array-setter C))
(A_ (array-getter A))
(B_ (array-getter B)))
(C! (+ (* (A_ 0 0) (B_ 0 0))
(* (A_ 0 1) (B_ 1 0)))
0 0)
(C! (+ (* (A_ 0 0) (B_ 0 1))
(* (A_ 0 1) (B_ 1 1)))
0 1)
(C! (+ (* (A_ 1 0) (B_ 0 0))
(* (A_ 1 1) (B_ 1 0)))
1 0)
(C! (+ (* (A_ 1 0) (B_ 0 1))
(* (A_ 1 1) (B_ 1 1)))
1 1)))
;;; Reshape A, B, and C to change all the 4-vectors to 2x2 matrices
(time
(array-for-each 2x2-matrix-multiply-into!
(array-curry (specialized-array-reshape A interval-2x2) 2)
(array-curry (specialized-array-reshape B interval-2x2) 2)
(array-curry (specialized-array-reshape C interval-2x2) 2)))
;;; Displays
;;; 0.015186 secs real time
;;; 0.015186 secs cpu time (0.015186 user, 0.000000 system)
;;; 4 collections accounting for 0.004735 secs real time (0.004732 user, 0.000000 system)
;;; 46089024 bytes allocated
;;; no minor faults
;;; no major faults
;;; Reshape each 4-vector to a 2x2 matrix individually
(time
(array-for-each (lambda (A B C)
(2x2-matrix-multiply-into!
(specialized-array-reshape A (make-interval '#(2 2)))
(specialized-array-reshape B (make-interval '#(2 2)))
(specialized-array-reshape C (make-interval '#(2 2)))))
(array-curry A 1)
(array-curry B 1)
(array-curry C 1)))
;;; Displays
;;; 0.039193 secs real time
;;; 0.039193 secs cpu time (0.039191 user, 0.000002 system)
;;; 6 collections accounting for 0.006855 secs real time (0.006851 user, 0.000001 system)
;;; 71043024 bytes allocated
;;; no minor faults
;;; no major faults
```

We provide a sample implementation in Gambit Scheme; the nonstandard techniques used
in the implementation are `define-structure`

, `define-macro`

, and DSSSL optional arguments. The sample implementation open codes specialized versions of algorithms for intervals and arrays of dimension no greater than 4, but a simple implementation could employ only the general algorithms that are used for dimensions greater than 4.

There is a git repository of this document, a sample implementation, a test file, and other materials.

NumPy has a popular array library, upon which PyTorch, a machine-learning tensor library, is built.

We list some of the correspondences between the routines of this SRFI and the array routines of NumPy. As we are in no way experienced NumPy programmers, we do not claim any level of completeness in this correspondence.

SRFI 231 | NumPy |
---|---|

array-assign! | copyto |

array-domain | shape |

specialized-array-reshape | reshape |

array-permute with index-first, index-last, index-rotate, index-swap | move-axis, transpose, swap-axes |

array-tile | split, array-split, dsplit, hsplit, vsplit, array slicing notation |

array-reverse | flip, fliplr, flipud |

Combine array-reverse and array-permute | rot90 |

array-append | concatenate |

array-stack, array-decurry | stack, vstack, hstack, dstack, column_stack, row_stack |

array-block | block |

array-map with array-copy, array-copy!, or array-assign! | Elementwise array operations |

array-translate | No correspondence (array indices always start at 0) |

array-sample | array slicing notation |

array-extract | array slicing notation |

array-copy and array-copy! | Use the copy method |

array-reduce | Use the reduce method |

array-outer-product | Use the outer method |

Combine array-curry, array-permute, and array-map | Operations specifying the axes along which they will be performed |

make-array and array-copy or array-copy! | fromfunction |

make-specialized-array-from-data | frombuffer, fromstring |

list->array, list*->array, vector->array, vector*->array | asarray, asanyarray |

array->list* | tolist |

Final SRFIs 25, 47, 58, and 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 procedures in these previous SRFIs have corresponding forms in this SRFI. For example, from SRFI 63, we can translate:

`(array? object)`

`(array? object)`

`(array-rank A)`

`(array-dimension A)`

`(make-array prototype k1 ...)`

`(make-specialized-array (make-interval (vector k1 ...)))`

.`(make-shared-array A mapper k1 ...)`

`(specialized-array-share A (make-interval (vector k1 ...)) mapper)`

`(array-in-bounds? A index1 ...)`

`(interval-contains-multi-index? (array-domain A) index1 ...)`

`(array-ref A k1 ...)`

`(let ((A_ (array-getter A))) ... (A_ k1 ...) ... )`

or`(array-ref A k1 ...)`

`(array-set! A object k1 ...)`

`(let ((A! (array-setter A))) ... (A! object k1 ...) ...)`

or`(array-set! A object k1 ...)`

Racket has an extensive array library, written by Neil Toronto, as part of its "Math Library". I do not claim to have a complete understanding of Racket's array library, but attempt here to give a superficial comparison of some aspects of Racket's library with this proposal:

- Racket's library has what it calls broadcasting and slicing; this proposal lacks these features as primitives.
- Racket's nonstrict arrays correspond to our "generalized arrays".
- Racket's arrays axes are indexed from zero; this SRFI allows nonzero lower bounds. Thus Racket's library has no need for array-translate.
- Racket's array-map is similar to the one in this proposal, except that ours always returns a generalized array.
- Racket's array-axis-ref can be implemented in this SRFI with
`array-permute`

and`array-curry`

. Racket's array-axis-permute is similar to our`array-permute`

. Both have`array-reshape`

and`array-append`

. Racket's array-flatten is similar to`(array->vector array)`

in this proposal. - Racket's array-axis-fold can be implemented in this SRFI as

If one wants what Racket calls a "strict" array as a result, apply`(define (array-axis-fold arr k f init) (array-map (lambda (pencil) (array-fold-left f init pencil)) (array-curry (array-permute arr (index-last (array-dimension arr) k)) 1)))`

`array-copy`

to the result. One can define Racket's "*-axis-*" procedures similarly. - Racket's library has specialized mathematical array operations for many math procedures; this library does not.
- Racket's library has flonum and complex flonum arrays; this library has similar features, including for various other homogeneous storage types, and is extendable.
- Racket has many procedures to select and recombine data from various axes of arrays, some of which can be simulated in this SRFI with
`array-permute`

,`array-curry`

, and`array-stack`

. - I don't see procedures in Racket's library corresponding to
`array-curry`

,`array-decurry`

,`array-tile`

,`array-block`

,`array-reverse`

, or`array-sample`

. - I don't see a procedure in Racket's library that corresponds to specialized-array-share in this SRFI.

Image processing applications provided significant motivation for this SRFI.

**Manipulating images in PGM format. **On a system with eight-bit chars, one
can write procedures 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:

```
(define make-pgm cons)
(define pgm-greys car)
(define pgm-pixels cdr)
(define (read-pgm file)
(define (read-pgm-object port)
(skip-white-space port)
(let ((o (read port)))
;; to skip the newline or next whitespace
(read-char port)
(if (eof-object? o)
(error "reached end of pgm file")
o)))
(define (skip-to-end-of-line port)
(let loop ((ch (read-char port)))
(if (not (eqv? ch #\newline))
(loop (read-char port)))))
(define (white-space? ch)
(case ch
((#\newline #\space #\tab) #t)
(else #f)))
(define (skip-white-space port)
(let ((ch (peek-char port)))
(cond ((white-space? ch)
(read-char port)
(skip-white-space port))
((eqv? ch #\#)
(skip-to-end-of-line port)
(skip-white-space port))
(else #f))))
;; The image file formats defined in netpbm
;; are problematical because they read the data
;; in the header as variable-length ISO-8859-1 text,
;; including arbitrary whitespace and comments,
;; and then they may read the rest of the file
;; as binary data.
;; So we give here a solution of how to deal
;; with these subtleties in Gambit Scheme.
(call-with-input-file
(list path: file
char-encoding: 'ISO-8859-1
eol-encoding: 'lf)
(lambda (port)
;; We're going to read text for a while,
;; then switch to binary.
;; So we need to turn off buffering until
;; we switch to binary.
(port-settings-set! port '(buffering: #f))
(let* ((header (read-pgm-object port))
(columns (read-pgm-object port))
(rows (read-pgm-object port))
(greys (read-pgm-object port)))
;; Now we switch back to buffering
;; to speed things up.
(port-settings-set! port '(buffering: #t))
(make-pgm
greys
(array-copy
(make-array
(make-interval (vector rows columns))
(cond ((or (eq? header 'p5)
(eq? header 'P5))
;; pgm binary
(if (< greys 256)
;; one byte/pixel
(lambda (i j)
(char->integer
(read-char port)))
;; two bytes/pixel,
;;little-endian
(lambda (i j)
(let* ((first-byte
(char->integer
(read-char port)))
(second-byte
(char->integer
(read-char port))))
(+ (* second-byte 256)
first-byte)))))
;; pgm ascii
((or (eq? header 'p2)
(eq? header 'P2))
(lambda (i j)
(read port)))
(else
(error "not a pgm file"))))
(if (< greys 256)
u8-storage-class
u16-storage-class)))))))
(define (write-pgm pgm-data file #!optional force-ascii)
(call-with-output-file
(list path: file
char-encoding: 'ISO-8859-1
eol-encoding: 'lf)
(lambda (port)
(let* ((greys
(pgm-greys pgm-data))
(pgm-array
(pgm-pixels pgm-data))
(domain
(array-domain pgm-array))
(rows
(fx- (interval-upper-bound domain 0)
(interval-lower-bound domain 0)))
(columns
(fx- (interval-upper-bound domain 1)
(interval-lower-bound domain 1))))
(if force-ascii
(display "P2" port)
(display "P5" port))
(newline port)
(display columns port) (display port)
(display rows port) (newline port)
(display greys port) (newline port)
(array-for-each
(if force-ascii
(let ((next-pixel-in-line 1))
(lambda (p)
(write p port)
(if (fxzero? (fxand next-pixel-in-line 15))
(begin
(newline port)
(set! next-pixel-in-line 1))
(begin
(display port)
(set! next-pixel-in-line
(fx+ 1 next-pixel-in-line))))))
(if (fx< greys 256)
(lambda (p)
(write-u8 p port))
(lambda (p)
(write-u8 (fxand p 255) port)
(write-u8 (fxarithmetic-shift-right p 8)
port))))
pgm-array)))))
```

One can write a procedure to convolve an image with a filter as follows:

```
(define (array-convolve source filter)
(let* ((source-domain
(array-domain source))
(S_
(array-getter source))
(filter-domain
(array-domain filter))
(F_
(array-getter filter))
(result-domain
(interval-dilate
source-domain
;; the left bound of an interval is an equality,
;; the right bound is an inequality, hence the
;; the difference in the following two expressions
(vector-map -
(interval-lower-bounds->vector filter-domain))
(vector-map (lambda (x)
(- 1 x))
(interval-upper-bounds->vector filter-domain)))))
(make-array result-domain
(lambda (i j)
(array-fold-left
(lambda (p q)
(+ p q))
0
(make-array
filter-domain
(lambda (k l)
(* (S_ (+ i k)
(+ j l))
(F_ k l))))))
)))
```

together with some filters

```
(define sharpen-filter
(list->array
(make-interval '#(-1 -1) '#(2 2))
'(0 -1 0
-1 5 -1
0 -1 0)))
(define edge-filter
(list->array
(make-interval '#(-1 -1) '#(2 2))
'(0 -1 0
-1 4 -1
0 -1 0)))
```

Our computations might results in pixel values outside the valid range, so we define

```
(define (round-and-clip pixel max-grey)
(max 0 (min (exact (round pixel)) max-grey)))
```

We can then compute edges and sharpen a test image as follows:

```
(define test-pgm (read-pgm "girl.pgm"))
(let ((greys (pgm-greys test-pgm)))
(write-pgm
(make-pgm
greys
(array-map (lambda (p)
(round-and-clip p greys))
(array-convolve
(pgm-pixels test-pgm)
sharpen-filter)))
"sharper-test.pgm"))
(let* ((greys (pgm-greys test-pgm))
(edge-array
(array-copy
(array-map
abs
(array-convolve
(pgm-pixels test-pgm)
edge-filter))))
(max-pixel
(array-fold-left max 0 edge-array))
(normalizer
(inexact (/ greys max-pixel))))
(write-pgm
(make-pgm
greys
(array-map (lambda (p)
(- greys
(round-and-clip (* p normalizer) greys)))
edge-array))
"edge-test.pgm"))
```

**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

, then one could get 1024 axial views, each $512\times512$, of this 3D body by `body``(array-curry `

; or 512 median views, each $1024\times512$, by `body` 2)`(array-curry (array-permute `

; or finally 512 frontal views, each again $1024\times512$ pixels, by `body` (index-first 3 1)) 2)`(array-curry (array-permute `

; see Anatomical plane. Note that you want to have the head up in both the median and frontal views, so we use `body` (index-first 3 2)) 2)`index-first`

to provide the appropriate permutations.

**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 these notes on wavelets and approximation theory for more details.

Using 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:

```
(define (all-second-differences image direction)
(let ((image-domain (array-domain image)))
(let loop ((i 1)
(result '()))
(let ((negative-scaled-direction
(vector-map (lambda (j)
(* -1 j i))
direction))
(twice-negative-scaled-direction
(vector-map (lambda (j)
(* -2 j i))
direction)))
(cond ((interval-intersect
image-domain
(interval-translate
image-domain
negative-scaled-direction)
(interval-translate
image-domain
twice-negative-scaled-direction))
=>
(lambda (subdomain)
(loop
(+ i 1)
(cons
(array-copy
(array-map
(lambda (f_i f_i+d f_i+2d)
(+ f_i+2d
(* -2. f_i+d)
f_i))
(array-extract
image
subdomain)
(array-extract
(array-translate
image
negative-scaled-direction)
subdomain)
(array-extract
(array-translate
image
twice-negative-scaled-direction)
subdomain)))
result))))
(else
(reverse result)))))))
```

We can define a small synthetic image of size 8x8 pixels and compute its second differences in various directions:

```
(define image
(array-copy
(make-array (make-interval '#(8 8))
(lambda (i j)
(exact->inexact (+ (* i i) (* j j)))))))
(define (expose difference-images)
(pretty-print (map (lambda (difference-image)
(list (array-domain difference-image)
(array->list* difference-image)))
difference-images)))
(begin
(display "\nOriginal image:\n")
(pretty-print (list (array-domain image)
(array->list* image)))
(display
"\nSecond-differences in the direction $k\times (1,0)$:\n")
(expose (all-second-differences image '#(1 0)))
(display
"\nSecond-differences in the direction $k\times (1,1)$:\n")
(expose (all-second-differences image '#(1 1)))
(display
"\nSecond-differences in the direction $k\times (1,-1)$:\n")
(expose (all-second-differences image '#(1 -1))))
```

On Gambit 4.8.5, this yields (after some hand editing):

Original image: (#<%%interval #53 dimension: 2 %%volume: 64 lower-bounds: #(0 0) upper-bounds: #(8 8)> ((0. 1. 4. 9. 16. 25. 36. 49.) (1. 2. 5. 10. 17. 26. 37. 50.) (4. 5. 8. 13. 20. 29. 40. 53.) (9. 10. 13. 18. 25. 34. 45. 58.) (16. 17. 20. 25. 32. 41. 52. 65.) (25. 26. 29. 34. 41. 50. 61. 74.) (36. 37. 40. 45. 52. 61. 72. 85.) (49. 50. 53. 58. 65. 74. 85. 98.))) Second-differences in the direction $k\times (1,0)$: ((#<##interval #2 lower-bounds: #(0 0) upper-bounds: #(6 8)> ((2. 2. 2. 2. 2. 2. 2. 2.) (2. 2. 2. 2. 2. 2. 2. 2.) (2. 2. 2. 2. 2. 2. 2. 2.) (2. 2. 2. 2. 2. 2. 2. 2.) (2. 2. 2. 2. 2. 2. 2. 2.) (2. 2. 2. 2. 2. 2. 2. 2.))) (#<##interval #3 lower-bounds: #(0 0) upper-bounds: #(4 8)> ((8. 8. 8. 8. 8. 8. 8. 8.) (8. 8. 8. 8. 8. 8. 8. 8.) (8. 8. 8. 8. 8. 8. 8. 8.) (8. 8. 8. 8. 8. 8. 8. 8.))) (#<##interval #4 lower-bounds: #(0 0) upper-bounds: #(2 8)> ((18. 18. 18. 18. 18. 18. 18. 18.) (18. 18. 18. 18. 18. 18. 18. 18.)))) Second-differences in the direction $k\times (1,1)$: ((#<##interval #5 lower-bounds: #(0 0) upper-bounds: #(6 6)> ((4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.))) (#<##interval #6 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((16. 16. 16. 16.) (16. 16. 16. 16.) (16. 16. 16. 16.) (16. 16. 16. 16.))) (#<##interval #7 lower-bounds: #(0 0) upper-bounds: #(2 2)> ((36. 36.) (36. 36.)))) Second-differences in the direction $k\times (1,-1)$: ((#<##interval #8 lower-bounds: #(0 2) upper-bounds: #(6 8)> ((4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.) (4. 4. 4. 4. 4. 4.))) (#<##interval #9 lower-bounds: #(0 4) upper-bounds: #(4 8)> ((16. 16. 16. 16.) (16. 16. 16. 16.) (16. 16. 16. 16.) (16. 16. 16. 16.))) (#<##interval #10 lower-bounds: #(0 6) upper-bounds: #(2 8)> ((36. 36.) (36. 36.))))

You 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.

**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 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:

```
(define (make-separable-transform 1D-transform)
(lambda (a)
(let ((n (array-dimension a)))
(do ((d 0 (fx+ d 1)))
((fx= d n))
(array-for-each
1D-transform
(array-curry (array-permute a (index-last n d)) 1))))))
```

Here we put each axis in turn at the end and then apply `1D-transform`

to each of the pencils along that axis.

Wavelet 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:

```
(define (recursively-apply-transform-and-downsample transform)
(lambda (a)
(let ((sample-vector (make-vector (array-dimension a) 2)))
(define (helper a)
(if (fx< 1 (interval-upper-bound (array-domain a) 0))
(begin
(transform a)
(helper (array-sample a sample-vector)))))
(helper a))))
(define (recursively-downsample-and-apply-transform transform)
(lambda (a)
(let ((sample-vector (make-vector (array-dimension a) 2)))
(define (helper a)
(if (fx< 1 (interval-upper-bound (array-domain a) 0))
(begin
(helper (array-sample a sample-vector))
(transform a))))
(helper a))))
```

By adding a single loop that calculates scaled sums and differences of adjacent elements in a one-dimensional array, we can define various Haar wavelet transforms as follows:

```
(define (1D-Haar-loop a)
(let ((a_ (array-getter a))
(a! (array-setter a))
(n (interval-upper-bound (array-domain a) 0)))
(do ((i 0 (fx+ i 2)))
((fx= i n))
(let* ((a_i (a_ i))
(a_i+1 (a_ (fx+ i 1)))
(scaled-sum (fl/ (fl+ a_i a_i+1) (flsqrt 2.0)))
(scaled-difference (fl/ (fl- a_i a_i+1) (flsqrt 2.0))))
(a! scaled-sum i)
(a! scaled-difference (fx+ i 1))))))
(define 1D-Haar-transform
(recursively-apply-transform-and-downsample 1D-Haar-loop))
(define 1D-Haar-inverse-transform
(recursively-downsample-and-apply-transform 1D-Haar-loop))
(define hyperbolic-Haar-transform
(make-separable-transform 1D-Haar-transform))
(define hyperbolic-Haar-inverse-transform
(make-separable-transform 1D-Haar-inverse-transform))
(define Haar-transform
(recursively-apply-transform-and-downsample
(make-separable-transform 1D-Haar-loop)))
(define Haar-inverse-transform
(recursively-downsample-and-apply-transform
(make-separable-transform 1D-Haar-loop)))
```

We 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:

```
(let ((image
(array-copy
(make-array (make-interval '#(4 4))
(lambda (i j)
(case i
((0) 1.)
((1) -1.)
(else 0.)))))))
(display "
Initial image:
")
(pretty-print (list (array-domain image)
(array->list* image)))
(hyperbolic-Haar-transform image)
(display "\nArray of hyperbolic Haar wavelet coefficients: \n")
(pretty-print (list (array-domain image)
(array->list* image)))
(hyperbolic-Haar-inverse-transform image)
(display "\nReconstructed image: \n")
(pretty-print (list (array-domain image)
(array->list* image))))
```

This yields:

Initial image: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((1. 1. 1. 1.) (-1. -1. -1. -1.) (0. 0. 0. 0.) (0. 0. 0. 0.))) Array of hyperbolic Haar wavelet coefficients: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((0. 0. 0. 0.) (2.8284271247461894 0. 0. 0.) (0. 0. 0. 0.) (0. 0. 0. 0.))) Reconstructed image: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((.9999999999999996 .9999999999999996 .9999999999999996 .9999999999999996) (-.9999999999999996 -.9999999999999996 -.9999999999999996 -.9999999999999996) (0. 0. 0. 0.) (0. 0. 0. 0.)))

In 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.

We can apply the (nonhyperbolic) Haar transform to the same image as follows:

```
(let ((image
(array-copy
(make-array (make-interval '#(4 4))
(lambda (i j)
(case i
((0) 1.)
((1) -1.)
(else 0.)))))))
(display "\nInitial image:\n")
(pretty-print (list (array-domain image)
(array->list* image)))
(Haar-transform image)
(display "\nArray of Haar wavelet coefficients: \n")
(pretty-print (list (array-domain image)
(array->list* image)))
(Haar-inverse-transform image)
(display "\nReconstructed image: \n")
(pretty-print (list (array-domain image)
(array->list* image))))
```

This yields:

Initial image: (#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((1. 1. 1. 1.) (-1. -1. -1. -1.) (0. 0. 0. 0.) (0. 0. 0. 0.))) Array of Haar wavelet coefficients: (#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((0. 0. 0. 0.) (1.9999999999999998 0. 1.9999999999999998 0.) (0. 0. 0. 0.) (0. 0. 0. 0.))) Reconstructed image: (#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)> ((.9999999999999997 .9999999999999997 .9999999999999997 .9999999999999997) (-.9999999999999997 -.9999999999999997 -.9999999999999997 -.9999999999999997) (0. 0. 0. 0.) (0. 0. 0. 0.)))

You 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.

**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.

Given 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-map`

, `array-assign!`

, `specialized-array-share`

, `array-extract`

, and `array-outer-product`

.

```
(define (LU-decomposition A)
;; Assumes the domain of A is [0,n)\times [0,n)
;; and that Gaussian elimination can be applied
;; without pivoting.
(let ((n
(interval-upper-bound (array-domain A) 0))
(A_
(array-getter A)))
(do ((i 0 (fx+ i 1)))
((= i (fx- n 1)) A)
(let* ((pivot
(A_ i i))
(column/row-domain
;; both will be one-dimensional
(make-interval (vector (+ i 1))
(vector n)))
(column
;; the column below the (i,i) entry
(specialized-array-share A
column/row-domain
(lambda (k)
(values k i))))
(row
;; the row to the right of the (i,i) entry
(specialized-array-share A
column/row-domain
(lambda (k)
(values i k))))
;; the subarray to the right and
;; below the (i,i) entry
(subarray
(array-extract
A (make-interval
(vector (fx+ i 1) (fx+ i 1))
(vector n n)))))
;; Compute multipliers.
(array-assign!
column
(array-map (lambda (x)
(/ x pivot))
column))
;; Subtract the outer product of i'th
;; row and column from the subarray.
(array-assign!
subarray
(array-map -
subarray
(array-outer-product * column row)))))))
```

We then define a $4\times 4$ Hilbert matrix:

```
(define A
(array-copy
(make-array (make-interval '#(4 4))
(lambda (i j)
(/ (+ 1 i j))))))
(define (array-display A)
(define (display-item x)
(display x) (display "\t"))
(newline)
(case (array-dimension A)
((1) (array-for-each display-item A) (newline))
((2) (array-for-each (lambda (row)
(array-for-each display-item row)
(newline))
(array-curry A 1)))
(else
(error "array-display can't handle > 2 dimensions: " A))))
(display "\nHilbert matrix:\n\n")
(array-display A)
;;; which displays:
;;; 1 1/2 1/3 1/4
;;; 1/2 1/3 1/4 1/5
;;; 1/3 1/4 1/5 1/6
;;; 1/4 1/5 1/6 1/7
(LU-decomposition A)
(display "\nLU decomposition of Hilbert matrix:\n\n")
(array-display A)
;;; which displays:
;;; 1 1/2 1/3 1/4
;;; 1/2 1/12 1/12 3/40
;;; 1/3 1 1/180 1/120
;;; 1/4 9/10 3/2 1/2800
```

We can now define matrix multiplication as follows to check our result:

```
;;; Functions to extract the lower- and upper-triangular
;;; matrices of the LU decomposition of A.
(define (L a)
(let ((a_ (array-getter a))
(d (array-domain a)))
(make-array
d
(lambda (i j)
(cond ((= i j) 1) ;; diagonal
((> i j) (a_ i j)) ;; below diagonal
(else 0)))))) ;; above diagonal
(define (U a)
(let ((a_ (array-getter a))
(d (array-domain a)))
(make-array
d
(lambda (i j)
(cond ((<= i j) (a_ i j)) ;; diagonal and above
(else 0)))))) ;; below diagonal
(display "\nLower triangular matrix of decomposition of Hilbert matrix:\n\n")
(array-display (L A))
;;; which displays:
;;; 1 0 0 0
;;; 1/2 1 0 0
;;; 1/3 1 1 0
;;; 1/4 9/10 3/2 1
(display "\nUpper triangular matrix of decomposition of Hilbert matrix:\n\n")
(array-display (U A))
;;; which displays:
;;; 1 1/2 1/3 1/4
;;; 0 1/12 1/12 3/40
;;; 0 0 1/180 1/120
;;; 0 0 0 1/2800
;;; We'll define a brief, not-very-efficient matrix multiply procedure.
(define (matrix-multiply A B)
(array-inner-product A + * B))
;;; We'll check that the product of the result of LU
;;; decomposition of A is again A.
(define product (matrix-multiply (L A) (U A)))
(display "\nProduct of lower and upper triangular matrices \n")
(display "of LU decomposition of Hilbert matrix:\n\n")
(array-display product)
;;; which displays:
;;; 1 1/2 1/3 1/4
;;; 1/2 1/3 1/4 1/5
;;; 1/3 1/4 1/5 1/6
;;; 1/4 1/5 1/6 1/7
```

**Conway's Game of Life. **Alex Harsányi implemented Conway's Game of Life using Racket's array library; here we implement the game using this SRFI.

Our strategy is to extend the original array periodically to an array dilated by one row and column above and below, left and right:

```
(define (array-pad-periodically a N)
;; Pad a periodically with N rows and columns top and bottom, left and right.
;; Assumes that the domain of a has zero lower bounds.
;; Returns a generalized array.
(let* ((domain (array-domain a))
(m (interval-upper-bound domain 0))
(n (interval-upper-bound domain 1))
(a_ (array-getter a)))
(make-array (interval-dilate domain (vector (- N) (- N)) (vector N N))
(lambda (i j)
(a_ (modulo i m) (modulo j n))))))
(define (neighbor-count a)
(let* ((big-a (array-copy (array-pad-periodically a 1)
(array-storage-class a)))
(domain (array-domain a))
(translates (map (lambda (translation)
(array-extract (array-translate big-a translation)
domain))
'(#(1 0) #(0 1) #(-1 0) #(0 -1)
#(1 1) #(1 -1) #(-1 1) #(-1 -1)))))
;; Returns a generalized array that contains the number
;; of 1s in the 8 cells surrounding each cell in the original array.
(apply array-map + translates)))
(define (game-rules a neighbor-count)
;; a is a single cell, neighbor-count is the count of 1s in
;; its 8 neighboring cells.
(if (= a 1)
(if (or (= neighbor-count 2)
(= neighbor-count 3))
1 0)
;; (= a 0)
(if (= neighbor-count 3)
1 0)))
(define (advance a)
;; Returns a specialized array
(array-copy
(array-map game-rules a (neighbor-count a))
(array-storage-class a)))
(define glider
(list*->array
2
'((0 0 0 0 0 0 0 0 0 0)
(0 0 1 0 0 0 0 0 0 0)
(0 0 0 1 0 0 0 0 0 0)
(0 1 1 1 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
u1-storage-class))
(define (generations a N)
(do ((i 0 (fx+ i 1))
(a a (advance a)))
((fx= i N))
(newline)
(pretty-print (array->list* a))))
(generations glider 5)
```

which prints

```
((0 0 0 0 0 0 0 0 0 0)
(0 0 1 0 0 0 0 0 0 0)
(0 0 0 1 0 0 0 0 0 0)
(0 1 1 1 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
((0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 1 0 1 0 0 0 0 0 0)
(0 0 1 1 0 0 0 0 0 0)
(0 0 1 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
((0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 1 0 0 0 0 0 0)
(0 1 0 1 0 0 0 0 0 0)
(0 0 1 1 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
((0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 1 0 0 0 0 0 0 0)
(0 0 0 1 1 0 0 0 0 0)
(0 0 1 1 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
((0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 1 0 0 0 0 0 0)
(0 0 0 0 1 0 0 0 0 0)
(0 0 1 1 1 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0)
(0 0 0 0 0 0 0 0 0 0))
```

**Inner products. **Our `array-inner-product`

procedure differs from that found in APL in a number of ways, including that the result is always an array, never a scalar.

We take some examples from the APLX Language Reference:

```
;; Examples from
;; http://microapl.com/apl_help/ch_020_020_880.htm
(define TABLE1
(list->array
(make-interval '#(3 2))
'(1 2
5 4
3 0)))
(define TABLE2
(list->array
(make-interval '#(2 4))
'(6 2 3 4
7 0 1 8)))
(array-display (inner-product TABLE1 + * TABLE2))
;;; Displays
;;; 20 2 5 20
;;; 58 10 19 52
;;; 18 6 9 12
(define X (list*->array 1 '(1 3 5 7)))
(define Y (list*->array 1 '(2 3 6 7)))
(array->list* (array-inner-product X + (lambda (x y) (if (= x y) 1 0)) Y))
=>
2
```

The SRFI author thanks Edinah K Gnang, John Cowan, Sudarshan S Chawathe, Jamison Hope, Per Bothner, Alex Shinn, Jens Axel Søgaard, and Marc Nieper-Wißkirchen for their comments and suggestions, and Arthur A. Gleckler, SRFI Editor, for his guidance and patience.

- "multi-dimensional arrays in R5RS?", by Alan Bawden.
- SRFI 4: Homogeneous Numeric Vector Datatypes, by Marc Feeley.
- SRFI 25: Multi-dimensional Array Primitives, by Jussi Piitulainen.
- SRFI 47: Array, by Aubrey Jaffer.
- SRFI 58: Array Notation, by Aubrey Jaffer.
- SRFI 63: Homogeneous and Heterogeneous Arrays, by Aubrey Jaffer.
- SRFI 164: Enhanced multi-dimensional Arrays, by Per Bothner.

© 2016, 2018, 2020, 2022 Bradley J Lucier. All Rights Reserved.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice (including the next paragraph) shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Editor: Arthur A. Gleckler