定义
通用数组切片可以通过一个 dope 向量或描述符引用每个数组来实现（无论是否内置在语言中）——一个包含第一个数组元素地址的记录，然后是每个索引的范围和指数公式中的相应系数。这种技术还允许立即数组转置、索引反转、子采样等。对于像 C 这样的语言，索引总是从零开始，具有 d 个索引的数组的涂料向量至少有 1 + 2d 个参数。
@987654321 @
这是一个密集的段落，但实际上都在里面。所以我们需要这样的数据结构：
struct {
TYPE *data; //address of first array element
int rank; //number of dimensions
int *dims; //size of each dimension
int *weight; //corresponding coefficient in the indexing formula
};
其中TYPE 是元素类型，即矩阵的字段。为简单起见，我们只使用int。出于我自己的目的，我设计了一种将各种类型的编码转换为integer handles，所以int 为me，YMMV 完成了这项工作。
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
所有指针成员都可以在
与结构本身相同的分配内存块（我们将
调用标题）。但是通过替换早期使用的偏移量
和struct-hackery，可以实现算法的实现
独立于内部（或外部）的实际内存布局
堵塞。
自包含数组对象的基本内存布局是
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
data[0] data[1] ... data[ product(dims)-1 ]
共享数据的间接数组（整个数组或 1 个或多个行切片）
将有以下内存布局
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
//no data! it's somewhere else!
还有一个包含正交切片的间接数组
另一个轴将具有与基本间接数组相同的布局，
但适当修改了尺寸和重量。
具有索引 (i0 i1 ... iN) 的元素的访问公式
是
a->data[ i0*a->weight[0] + i1*a->weight[1] + ...
+ iN*a->weight[N] ]
，假设每个索引 i[j] 在 [ 0 ... dims[j] 之间）。
在一个正常布局的row-major 数组的weight 向量中，每个元素都是所有低维的乘积。
for (int i=0; i<rank; i++)
weight[i] = product(dims[i+1 .. rank-1]);
所以对于一个 3×4×5 的数组，元数据应该是
{ .rank=3, .dims=(int[]){3,4,5}, .weight=(int[]){4*5, 5, 1} }
或者对于一个 7×6×5×4×3×2 的数组，元数据是
{ .rank=6, .dims={7,6,5,4,3,2}, .weight={720, 120, 24, 6, 2, 1} }
建设
因此，要创建其中之一，我们需要来自 previous question 的相同辅助函数来计算维度列表的大小。
/* multiply together rank integers in dims array */
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}
要分配，只需malloc 足够的内存并将指针设置到适当的位置。
/* create array given rank and int[] dims */
arr arraya(int rank, int dims[]){
int datasz;
int i;
int x;
arr z;
datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));
z->rank = rank;
z->dims = z + 1;
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
使用与上一个答案相同的技巧，我们可以制作一个可变参数接口以简化使用。
/* load rank integers from va_list into int[] dims */
void loaddimsv(int rank, int dims[], va_list ap){
int i;
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
}
/* create a new array with specified rank and dimensions */
arr (array)(int rank, ...){
va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];
int i;
int x;
arr z;
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
z = arraya(rank,dims);
//free(dims);
return z;
}
甚至通过使用 ppnarg 的强大功能计算其他参数来自动生成 rank 参数。
#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__) /* create a new array with specified dimensions */
现在构建其中之一非常容易。
arr a = array(2,3,4); // create a dynamic [2][3][4] array
访问元素
通过对elema 的函数调用检索元素，该函数将每个索引乘以相应的权重，将它们相加，并索引data 指针。我们返回一个指向元素的指针，以便调用者可以读取或修改它。
/* access element of a indexed by int[] */
int *elema(arr a, int *ind){
int idx = 0;
int i;
for (i=0; i<a->rank; i++){
idx += ind[i] * a->weight[i];
}
return a->data + idx;
}
同样的可变参数技巧可以使界面更简单elem(a,i,j,k)。
轴向切片
要进行切片，首先我们需要一种方法来指定要提取哪些维度以及要折叠哪些维度。如果我们只需要从一个维度中选择单个索引或所有元素，那么slice 函数可以将 rank ints 作为参数，并将 -1 解释为选择整个维度或 0。 .dimsi-1 作为选择单个索引。
/* take a computed slice of a following spec[] instructions
if spec[i] >= 0 and spec[i] < a->rank, then spec[i] selects
that index from dimension i.
if spec[i] == -1, then spec[i] selects the entire dimension i.
*/
arr slicea(arr a, int spec[]){
int i,j;
int rank;
for (i=0,rank=0; i<a->rank; i++)
rank+=spec[i]==-1;
int dims[rank];
int weight[rank];
for (i=0,j=0; i<rank; i++,j++){
while (spec[j]!=-1) j++;
if (j>=a->rank) break;
dims[i] = a->dims[j];
weight[i] = a->weight[j];
}
arr z = casta(a->data, rank, dims);
memcpy(z->weight,weight,rank*sizeof(int));
for (j=0; j<a->rank; j++){
if (spec[j]!=-1)
z->data += spec[j] * a->weight[j];
}
return z;
}
所有与参数数组中的 -1 对应的维度和权重都被收集并用于创建新的数组头。所有 >= 0 的参数都乘以它们的相关权重并添加到 data 指针，偏移指向正确元素的指针。
就数组访问公式而言，我们将其视为多项式。
offset = constant + sum_i=0,n( weight[i] * index[i] )
因此，对于我们从中选择单个元素的任何维度（+ 所有较低维度），我们会分解现在的常数索引并将其添加到公式中的常数项（在我们的 C 表示中是data 指针本身）。
辅助函数casta 使用共享data 创建新的数组头。 slicea 当然会改变权重值，但是通过计算权重本身，casta 成为更普遍可用的功能。它甚至可以用来构造一个动态数组结构，直接在一个常规的 C 风格的多维数组上操作，从而casting。
/* create an array header to access existing data in multidimensional layout */
arr casta(int *data, int rank, int dims[]){
int i,x;
arr z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = z + 1;
z->weight = z->dims + rank;
z->data = data;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
转置
涂料向量也可用于实现转置。维度（和索引）的顺序可以更改。
请记住，这不是像其他人一样的正常“转置”
做。我们根本不重新排列数据。这更
恰当地称为“涂料矢量伪转置”。
我们不改变数据，而是改变
访问公式中的常量，重新排列
多项式的系数。在真正意义上，这
只是交换律的一个应用
加法的结合性。
因此，为了具体起见，假设数据已排列
依次从假设地址 500 开始。
500: 0
501: 1
502: 2
503: 3
504: 4
505: 5
506: 6
如果 a 是 rank 2，dims {1, 7)，weight (7, 1)，那么
指数总和乘以相关权重
添加到初始指针 (500) 产生适当的
每个元素的地址
a[0][0] == *(500+0*7+0*1)
a[0][1] == *(500+0*7+1*1)
a[0][2] == *(500+0*7+2*1)
a[0][3] == *(500+0*7+3*1)
a[0][4] == *(500+0*7+4*1)
a[0][5] == *(500+0*7+5*1)
a[0][6] == *(500+0*7+6*1)
所以涂料向量伪转置重新排列
重量和尺寸以匹配新的排序
指数，但总和保持不变。最初的
指针保持不变。数据不移动。
b[0][0] == *(500+0*1+0*7)
b[1][0] == *(500+1*1+0*7)
b[2][0] == *(500+2*1+0*7)
b[3][0] == *(500+3*1+0*7)
b[4][0] == *(500+4*1+0*7)
b[5][0] == *(500+5*1+0*7)
b[6][0] == *(500+6*1+0*7)
内积（又称矩阵乘法）
因此，通过使用通用切片或转置+“行”-切片（更容易），可以实现广义内积。
首先，我们需要两个辅助函数，用于对两个向量应用二元运算以产生向量结果，并通过二元运算对向量进行归约以产生标量结果。
就像在previous question 中一样，我们将传入运算符，因此同一个函数可以与许多不同的运算符一起使用。对于这里的风格，我将运算符作为单个字符传递，因此已经存在从 C 运算符到 我们函数的 运算符的间接映射。这是x-macro table。
#define OPERATORS(_) \
/* f F id */ \
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
/**/
#define binop(X,F,Y) (binop)(X,*#F,Y)
arr (binop)(arr x, char f, arr y); /* perform binary operation F upon corresponding elements of vectors X and Y */
表中的额外元素是针对空向量参数情况下的reduce 函数。在这种情况下，reduce 应该返回操作符的身份元素，+ 为 0，* 为 1。
#define reduce(F,X) (reduce)(*#F,X)
int (reduce)(char f, arr a); /* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
所以binop 在操作符上执行循环和开关。
/* perform binary operation F upon corresponding elements of vectors X and Y */
#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break;
arr (binop)(arr x, char f, arr y){
arr z=copy(x);
int n=x->dims[0];
int i;
for (i=0; i<n; i++){
switch(f){
OPERATORS(BINOP)
}
}
return z;
}
#undef BINOP
如果有足够的元素，reduce 函数会执行向后循环，如果有，则将初始值设置为最后一个元素，并将初始值预设为运算符的标识元素。
/* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
#define REDID(f,F,id) case f: x = id; break;
#define REDOP(f,F,id) case f: x = *elem(a,i) F x; break;
int (reduce)(char f, arr a){
int n = a->dims[0];
int x;
int i;
switch(f){
OPERATORS(REDID)
}
if (n) {
x=*elem(a,n-1);
for (i=n-2;i>=0;i--){
switch(f){
OPERATORS(REDOP)
}
}
}
return x;
}
#undef REDID
#undef REDOP
使用这些工具，可以以更高级别的方式实现内部产品。
/* perform a (2D) matrix multiplication upon rows of x and columns of y
using operations F and G.
Z = X F.G Y
Z[i,j] = F/ X[i,*] G Y'[j,*]
more generally,
perform an inner product on arguments of compatible dimension.
Z = X[A;B;C;D;E;F] +.* Y[G;H;I;J;K] |(F = G)
Z[A;B;C;D;E;H;I;J;K] = +/ X[A;B;C;D;E;*] * Y[*;H;I;J;K]
*/
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr xdims = cast(x->dims,1,x->rank);
arr ydims = cast(y->dims,1,y->rank);
xdims->dims[0]--;
ydims->dims[0]--;
ydims->data++;
arr z=arraya(x->rank+y->rank-2,catv(xdims,ydims)->data);
int datasz = productdims(z->rank,z->dims);
int k=y->dims[0];
arr xs = NULL;
arr ys = NULL;
for (i=0; i<datasz; i++){
int idx[x->rank+y->rank];
vector_index(i,z->dims,z->rank,idx);
int *xdex=idx;
int *ydex=idx+x->rank-1;
memmove(ydex+1,ydex,y->rank);
xdex[x->rank-1] = -1;
free(xs);
free(ys);
xs = slicea(x,xdex);
ys = slicea(y,ydex);
z->data[i] = (reduce)(f,(binop)(xs,g,ys));
}
free(xs);
free(ys);
free(xdims);
free(ydims);
return z;
}
此函数还使用函数cast，它为casta 提供可变参数接口。
/* create an array header to access existing data in multidimensional layout */
arr cast(int *data, int rank, ...){
va_list ap;
int dims[rank];
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
return casta(data, rank, dims);
}
它还使用vector_index 将一维索引转换为索引的nD 向量。
/* compute vector index list for ravel index ind */
int *vector_index(int ind, int *dims, int n, int *vec){
int i,t=ind, *z=vec;
for (i=0; i<n; i++){
z[n-1-i] = t % dims[n-1-i];
t /= dims[n-1-i];
}
return z;
}
github file。其他支持功能也在 github 文件中。
这个 Q/A 对是在实现 inca 一个用 C 编写的 APL 语言的解释器时出现的一系列相关问题的一部分。其他：How can I work with dynamically-allocated arbitrary-dimensional arrays? 和 How to pass a C math operator (+-*/%) into a function result=mathfunc(x,+,y);?。其中一些材料之前已发布到comp.lang.c。 comp.lang.apl 中的更多背景信息。