使用涂料向量访问多维数组的任意轴向切片？
Use a dope vector to access arbitrary axial slices of a multidimensional array?
我正在构建一套函数来处理 ，我希望能够定义数组的任意 切片 ，这样我就可以实现一个两个任意矩阵（又名 张量 或 n-d 数组）的广义内积。
我读过的一篇 APL 论文（老实说我找不到——我读了很多）在左矩阵 X 上定义了矩阵乘积，维度为 A;B;C;D;E;F，右矩阵为 A;B;C;D;E;F维度 G;H;I;J;K 的矩阵 Y 其中 F==G 为
Z <- X +.× Y
Z[A;B;C;D;E;H;I;J;K] <- +/ X[A;B;C;D;E;*] × Y[*;H;I;J;K]
其中 +/ 是 的和，× 将逐个元素应用于两个相同长度的向量。
所以我需要左侧 "row" 片和右侧 "column" 片。我当然可以进行转置，然后进行 "row" 切片来模拟 "column" 切片，但我宁愿做得更优雅。
维基百科关于 slicing leads to a stub about dope vectors 的文章似乎是我正在寻找的灵丹妙药，但没有太多内容可以继续。
如何使用涂料矢量实现任意切片？
（很久以后我才注意到 Stride of an array 里面有一些细节。）
定义
General array slicing can be implemented (whether or not built into the language) by referencing every array through a dope vector or descriptor — a record that contains the address of the first array element, and then the range of each index and the corresponding coefficient in the indexing formula. This technique also allows immediate array transposition, index reversal, subsampling, etc. For languages like C, where the indices always start at zero, the dope vector of an array with d indices has at least 1 + 2d parameters.
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。出于我自己的目的，我设计了一种将各种类型编码为 integer handles 的编码，因此 int 为 me 完成了工作，YMMV。
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
所有指针成员都可以在
与结构本身相同的分配内存块（我们将
调用 header)。但是通过替换早期使用的偏移量
和struct-hackery，可以实现算法
独立于内部（或外部）的实际内存布局
堵塞。
self-contained 数组 object 的基本内存布局是
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 个或多个 row-slices）
将具有以下内存布局
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] ) 之间。
在通常 laid-out 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;
}
并使用与上一个答案相同的技巧，我们可以制作一个 variable-argument 界面以简化使用。
/* 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 指针。我们 return 指向元素的指针，以便调用者可以读取或修改它。
/* 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)。
轴向切片
要进行切片，首先我们需要一种方法来指定要提取哪些维度以及折叠哪些维度。如果我们只需要 select 单个索引或维度中的所有元素，那么 slice 函数可以将 rank ints 作为参数并且将 -1 解释为 selecting 整个维度或 0..dimsi-1 作为 selecting 单个维度指数。
/* 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对应的所有维度和权重，并用于创建新数组header。所有 >= 0 的参数都乘以它们的相关权重并添加到 data 指针，offsetting 指向正确元素的指针。
就数组访问公式而言，我们将其视为多项式。
offset = constant + sum_i=0,n( weight[i] * index[i] )
因此，对于我们从中 select 获取单个元素（+所有较低维度）的任何维度，我们 factor-out now-constant 索引并将其添加到常数项在公式中（在我们的 C 表示中是 data 指针本身）。
辅助函数 casta 创建新数组 header 并共享 data。 slicea 当然改变了权重值，但是通过计算权重本身，casta 变成了一个更普遍可用的函数。它甚至可以用来构造一个动态数组结构，直接在一个常规的 C-style 多维数组上运行，因此 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;
}
转置
掺杂向量也可用于实现转置。可以更改维度（和索引）的顺序。
请记住，这与其他人不同 'transpose'
做。我们根本不重新排列数据。这是更多
正确地称为 'dope-vector pseudo-transpose'。
我们不改变数据，而是改变
访问公式中的常量，重新排列
多项式的系数。真正意义上，这
只是交换律的一个应用，并且
加法的结合性。
所以为了具体起见，假设数据已经安排好
按顺序从假设地址 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)
所以 dope-vectorpseudo-transpose 重新排列
重量和尺寸以匹配新的订购
指数，但总和保持不变。最初的
指针保持不变。数据不会移动。
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)
内积（又名矩阵乘法）
因此，通过使用通用切片或转置+"row"-切片（更容易），可以实现广义内积。
首先，我们需要两个辅助函数，用于将二元运算应用于两个向量以产生向量结果，并使用二元运算减少向量以产生标量结果。
如 previous question we'll pass in the operator, so the same function can be used with many different operators. For the style here, I'm passing operators as single characters, so there's already an indirect mapping from C operators to our function's operators. This is an 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 */
table 中的额外元素用于空向量参数情况下的 reduce 函数。在这种情况下，reduce 应该 return 运算符的 标识元素 ，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
使用这些工具，内积可以higher-level方式实现。
/* 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 将 1D 索引转换为 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 an interpreter for the APL language written in C. Others: How can I work with dynamically-allocated arbitrary-dimensional arrays? , and How to pass a C math operator (+-*/%) into a function result=mathfunc(x,+,y);? . Some of this material was previously posted to comp.lang.c. More background in comp.lang.apl.
过程中出现的一系列相关问题的一部分
我正在构建一套函数来处理
我读过的一篇 APL 论文（老实说我找不到——我读了很多）在左矩阵 X 上定义了矩阵乘积，维度为 A;B;C;D;E;F，右矩阵为 A;B;C;D;E;F维度 G;H;I;J;K 的矩阵 Y 其中 F==G 为
Z <- X +.× Y
Z[A;B;C;D;E;H;I;J;K] <- +/ X[A;B;C;D;E;*] × Y[*;H;I;J;K]
其中 +/ 是 的和，× 将逐个元素应用于两个相同长度的向量。
所以我需要左侧 "row" 片和右侧 "column" 片。我当然可以进行转置，然后进行 "row" 切片来模拟 "column" 切片，但我宁愿做得更优雅。
维基百科关于 slicing leads to a stub about dope vectors 的文章似乎是我正在寻找的灵丹妙药，但没有太多内容可以继续。
如何使用涂料矢量实现任意切片？
（很久以后我才注意到 Stride of an array 里面有一些细节。）
定义
General array slicing can be implemented (whether or not built into the language) by referencing every array through a dope vector or descriptor — a record that contains the address of the first array element, and then the range of each index and the corresponding coefficient in the indexing formula. This technique also allows immediate array transposition, index reversal, subsampling, etc. For languages like C, where the indices always start at zero, the dope vector of an array with d indices has at least 1 + 2d parameters.
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。出于我自己的目的，我设计了一种将各种类型编码为 integer handles 的编码，因此 int 为 me 完成了工作，YMMV。
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
所有指针成员都可以在 与结构本身相同的分配内存块（我们将 调用 header)。但是通过替换早期使用的偏移量 和struct-hackery，可以实现算法 独立于内部（或外部）的实际内存布局 堵塞。
self-contained 数组 object 的基本内存布局是
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 个或多个 row-slices） 将具有以下内存布局
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] ) 之间。
在通常 laid-out 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;
}
并使用与上一个答案相同的技巧，我们可以制作一个 variable-argument 界面以简化使用。
/* 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 指针。我们 return 指向元素的指针，以便调用者可以读取或修改它。
/* 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)。
轴向切片
要进行切片，首先我们需要一种方法来指定要提取哪些维度以及折叠哪些维度。如果我们只需要 select 单个索引或维度中的所有元素，那么 slice 函数可以将 rank ints 作为参数并且将 -1 解释为 selecting 整个维度或 0..dimsi-1 作为 selecting 单个维度指数。
/* 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对应的所有维度和权重，并用于创建新数组header。所有 >= 0 的参数都乘以它们的相关权重并添加到 data 指针，offsetting 指向正确元素的指针。
就数组访问公式而言，我们将其视为多项式。
offset = constant + sum_i=0,n( weight[i] * index[i] )
因此，对于我们从中 select 获取单个元素（+所有较低维度）的任何维度，我们 factor-out now-constant 索引并将其添加到常数项在公式中（在我们的 C 表示中是 data 指针本身）。
辅助函数 casta 创建新数组 header 并共享 data。 slicea 当然改变了权重值，但是通过计算权重本身，casta 变成了一个更普遍可用的函数。它甚至可以用来构造一个动态数组结构，直接在一个常规的 C-style 多维数组上运行，因此 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;
}
转置
掺杂向量也可用于实现转置。可以更改维度（和索引）的顺序。
请记住，这与其他人不同 'transpose' 做。我们根本不重新排列数据。这是更多 正确地称为 'dope-vector pseudo-transpose'。 我们不改变数据，而是改变 访问公式中的常量，重新排列 多项式的系数。真正意义上，这 只是交换律的一个应用，并且 加法的结合性。
所以为了具体起见，假设数据已经安排好 按顺序从假设地址 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)
所以 dope-vectorpseudo-transpose 重新排列 重量和尺寸以匹配新的订购 指数，但总和保持不变。最初的 指针保持不变。数据不会移动。
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)
内积（又名矩阵乘法）
因此，通过使用通用切片或转置+"row"-切片（更容易），可以实现广义内积。
首先，我们需要两个辅助函数，用于将二元运算应用于两个向量以产生向量结果，并使用二元运算减少向量以产生标量结果。
如 previous question we'll pass in the operator, so the same function can be used with many different operators. For the style here, I'm passing operators as single characters, so there's already an indirect mapping from C operators to our function's operators. This is an 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 */
table 中的额外元素用于空向量参数情况下的 reduce 函数。在这种情况下，reduce 应该 return 运算符的 标识元素 ，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
使用这些工具，内积可以higher-level方式实现。
/* 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 将 1D 索引转换为 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 an interpreter for the APL language written in C. Others: How can I work with dynamically-allocated arbitrary-dimensional arrays? , and How to pass a C math operator (+-*/%) into a function result=mathfunc(x,+,y);? . Some of this material was previously posted to comp.lang.c. More background in comp.lang.apl.
过程中出现的一系列相关问题的一部分