定义
通用数组切片可以通过引用每个数组通过一个dope向量或描述符来实现（无论是否内置于语言中）- 包含第一个数组元素的地址，然后是每个索引的范围和相应的索引公式中的系数。该技术还允许立即进行数组转置、索引反转、子采样等操作。对于像C这样的语言，其中索引始终从零开始，具有d个索引的数组的dope向量至少具有1 + 2d个参数。
http://en.wikipedia.org/wiki/Array_slicing#Details
这是一个密集的段落，但实际上所有内容都在里面。因此我们需要这样的数据结构：
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。为了我的目的，我已经将各种类型编码成整数句柄，因此int可以胜任我的工作，但你的情况可能有所不同。
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
所有指针成员都可以被分配到与结构本身相同的内存块中（我们称之为头部）。但是通过替换早期使用的偏移量和结构体技巧，算法的实现可以独立于块内或块外的实际内存布局。
自包含数组对象的基本内存布局如下：
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!
如果一个间接数组包含沿另一个轴的正交切片，则其布局与基本间接数组相同，但dims和weight会适当修改。
具有索引（i0 i1 ... iN）的元素的访问公式为
a->data[ i0*a->weight[0] + i1*a->weight[1] + ...
+ iN*a->weight[N] ]
假设每个索引 i[j] 都在 [0 ... dims[j]) 范围内。
对于一个按行主序(normal 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} }
构建
因此，要创建其中一个，我们需要来自之前问题的相同辅助函数来从维度列表计算大小。
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}
要分配内存，只需调用
malloc函数分配足够的内存，并将指针设置到相应的位置即可。
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;
}
使用前面答案中的同样技巧，我们可以创建一个可变参数接口，使得使用更加容易。
void loaddimsv(int rank, int dims[], va_list ap){
int i;
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
}
arr (array)(int rank, ...){
va_list ap;
int dims[rank];
int i;
int x;
arr z;
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
z = arraya(rank,dims);
return z;
}
甚至可以利用 ppnarg 的强大功能通过计算其他参数来自动产生 rank 参数。
#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__)
现在构建其中之一非常容易。
arr a = array(2,3,4); // create a dynamic [2][3][4] array
访问元素
通过调用elema函数来检索一个元素，该函数将每个索引乘以相应的权重，求和后索引data指针。我们返回一个指向元素的指针，以便调用者可以读取或修改它。
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函数可以将
rankint作为参数，并将-1解释为选择整个维度，0..
dimsi-1解释为选择单个索引。
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风格多维数组的动态数组结构，从而进行转换。
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;
}
转置
“dope vector”还可以用于实现转置。可以改变维度（和索引）的顺序。
请记住，这不是像其他人一样的普通“转置”。我们根本不重新排列数据。更准确地说，这应该称为“dope-vector伪转置”。
我们只是改变访问公式中的常数，重新排列多项式的系数。在实质上，这只是加法的可交换性和结合律的应用。
因此，为了具体起见，假设数据按顺序排列，从假想地址500开始。
500: 0
501: 1
502: 2
503: 3
504: 4
505: 5
506: 6
如果a的秩为2，维度为{1,7}，权重为(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)
因此，dope-vector 伪转置重新排列权重和维度以匹配新的索引顺序，但总和保持不变。初始指针保持不变。数据不会移动。
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)
内积（也称矩阵乘法）
因此，通过使用一般切片或转置+“行”切片（更容易），可以实现广义内积。
首先，我们需要两个辅助函数，将二元操作应用于产生向量结果的两个向量，并使用二元操作减少向量以产生标量结果。
如在上一个问题中所述，我们将传递运算符，因此同一函数可与许多不同的运算符一起使用。对于这里的样式，我将运算符作为单个字符传递，因此从C运算符到我们函数的运算符已经存在间接映射。这是一个x宏表。
#define OPERATORS(_) \
\
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
#define binop(X,F,Y) (binop)(X,*#F,Y)
arr (binop)(arr x, char f, arr y);
表格中的额外元素是为了处理空向量参数的
reduce恒等元素+*函数。在这种情况下，reduce应该返回运算符的，对于来说是0，对于来说是1。
#define reduce(F,X) (reduce)(*#F,X)
int (reduce)(char f, arr a);
因此，
binop在运算符字符上执行循环和开关操作。
#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函数会进行反向循环，并将初始值设为最后一个元素（如果存在），同时预设初始值为运算符的身份元素。
#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 提供了一个可变参数接口。
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将1D索引转换为nD索引向量。
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文件。其他支持函数也在Github文件中。
这个问答对是一系列相关问题的一部分，这些问题出现在实现C语言编写的APL语言解释器inca中。其他问题: 如何处理动态分配的任意维数组？ 和 如何将C数学运算符（+-*/%）传递到函数result=mathfunc(x,+,y);中？。其中一些材料以前发布在comp.lang.c上。更多背景资料请参见comp.lang.apl。
