Move the outer function into a separate file.
This commit is contained in:
George Weigt
parent
3fa9fe52b3
commit
27a7c264f9
6 changed files with 122 additions and 167 deletions
2
Makefile
2
Makefile
|
|
@ -39,7 +39,7 @@ sin.o cos.o tan.o arcsin.o arccos.o arctan.o \
|
|||
sinh.o cosh.o tanh.o arcsinh.o arccosh.o arctanh.o \
|
||||
abs.o mod.o roots.o eigen.o simplify.o for.o isprime.o index.o wedge.o \
|
||||
rationalize.o prog.o lcm.o floor.o ceiling.o condense.o userfunc.o find.o \
|
||||
init.o primetab.o bignum.o symbol.o run.o atomize.o pollard.o \
|
||||
init.o primetab.o bignum.o symbol.o run.o atomize.o pollard.o outer.o \
|
||||
coeff.o \
|
||||
main.o \
|
||||
misc.o \
|
||||
|
|
|
|||
15
eval.cpp
15
eval.cpp
|
|
@ -27,6 +27,7 @@ extern void eval_filter(void);
|
|||
extern void eval_floor(void);
|
||||
extern void eval_isprime(void);
|
||||
extern void eval_mod(void);
|
||||
extern void eval_outer(void);
|
||||
extern void eval_product(void);
|
||||
extern void eval_roots(void);
|
||||
extern void eval_sample(void);
|
||||
|
|
@ -638,20 +639,6 @@ eval_operator(void)
|
|||
list(tos - h);
|
||||
}
|
||||
|
||||
static void
|
||||
eval_outer(void)
|
||||
{
|
||||
push(cadr(p1));
|
||||
eval();
|
||||
p1 = cddr(p1);
|
||||
while (iscons(p1)) {
|
||||
push(car(p1));
|
||||
eval();
|
||||
outer();
|
||||
p1 = cdr(p1);
|
||||
}
|
||||
}
|
||||
|
||||
static void
|
||||
eval_power(void)
|
||||
{
|
||||
|
|
|
|||
118
outer.cpp
Normal file
118
outer.cpp
Normal file
|
|
@ -0,0 +1,118 @@
|
|||
// outer product of tensors
|
||||
|
||||
#include "stdafx.h"
|
||||
#include "defs.h"
|
||||
|
||||
static void outer_f(void);
|
||||
|
||||
void
|
||||
eval_outer(void)
|
||||
{
|
||||
p1 = cdr(p1);
|
||||
push(car(p1));
|
||||
eval();
|
||||
p1 = cdr(p1);
|
||||
while (iscons(p1)) {
|
||||
push(car(p1));
|
||||
eval();
|
||||
outer();
|
||||
p1 = cdr(p1);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
outer(void)
|
||||
{
|
||||
save();
|
||||
p2 = pop();
|
||||
p1 = pop();
|
||||
if (istensor(p1) && istensor(p2))
|
||||
outer_f();
|
||||
else {
|
||||
push(p1);
|
||||
push(p2);
|
||||
if (istensor(p1))
|
||||
tensor_times_scalar();
|
||||
else if (istensor(p2))
|
||||
scalar_times_tensor();
|
||||
else
|
||||
multiply();
|
||||
}
|
||||
restore();
|
||||
}
|
||||
|
||||
void
|
||||
outer_f(void)
|
||||
{
|
||||
int i, j, k, ndim, nelem;
|
||||
|
||||
ndim = p1->u.tensor->ndim + p2->u.tensor->ndim;
|
||||
|
||||
if (ndim > MAXDIM)
|
||||
maxdim_error();
|
||||
|
||||
nelem = p1->u.tensor->nelem * p2->u.tensor->nelem;
|
||||
|
||||
p3 = alloc_tensor(nelem);
|
||||
|
||||
p3->u.tensor->ndim = ndim;
|
||||
|
||||
for (i = 0; i < p1->u.tensor->ndim; i++)
|
||||
p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
|
||||
|
||||
j = i;
|
||||
|
||||
for (i = 0; i < p2->u.tensor->ndim; i++)
|
||||
p3->u.tensor->dim[j + i] = p2->u.tensor->dim[i];
|
||||
|
||||
k = 0;
|
||||
|
||||
for (i = 0; i < p1->u.tensor->nelem; i++)
|
||||
for (j = 0; j < p2->u.tensor->nelem; j++) {
|
||||
push(p1->u.tensor->elem[i]);
|
||||
push(p2->u.tensor->elem[j]);
|
||||
multiply();
|
||||
p3->u.tensor->elem[k++] = pop();
|
||||
}
|
||||
|
||||
push(p3);
|
||||
}
|
||||
|
||||
static char *s[] = {
|
||||
|
||||
"outer(a,b)",
|
||||
"a*b",
|
||||
|
||||
"outer(a,(b1,b2))",
|
||||
"(a*b1,a*b2)",
|
||||
|
||||
"outer((a1,a2),b)",
|
||||
"(a1*b,a2*b)",
|
||||
|
||||
"H33=hilbert(3)",
|
||||
"",
|
||||
|
||||
"H44=hilbert(4)",
|
||||
"",
|
||||
|
||||
"H55=hilbert(5)",
|
||||
"",
|
||||
|
||||
"H3344=outer(H33,H44)",
|
||||
"",
|
||||
|
||||
"H4455=outer(H44,H55)",
|
||||
"",
|
||||
|
||||
"H33444455=outer(H33,H44,H44,H55)",
|
||||
"",
|
||||
|
||||
"inner(H3344,H4455)-contract(H33444455,4,5)",
|
||||
"0",
|
||||
};
|
||||
|
||||
void
|
||||
test_outer(void)
|
||||
{
|
||||
test(__FILE__, s, sizeof s / sizeof (char *));
|
||||
}
|
||||
|
|
@ -38,6 +38,7 @@ extern void test_mprime(void);
|
|||
extern void test_mroot(void);
|
||||
extern void test_msub(void);
|
||||
extern void test_multiply(void);
|
||||
extern void test_outer(void);
|
||||
extern void test_power(void);
|
||||
extern void test_product(void);
|
||||
extern void test_prog(void);
|
||||
|
|
@ -112,6 +113,7 @@ selftest(void)
|
|||
test_for();
|
||||
test_lcm();
|
||||
test_mod();
|
||||
test_outer();
|
||||
test_product();
|
||||
test_prog();
|
||||
test_rationalize();
|
||||
|
|
|
|||
110
tensor.cpp
110
tensor.cpp
|
|
@ -276,116 +276,6 @@ cofactor(U *p, int n, int row, int col)
|
|||
negate();
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
//
|
||||
// outer product of tos and tos-1
|
||||
//
|
||||
// the args themselves may also be outer products
|
||||
//
|
||||
//-----------------------------------------------------------------------------
|
||||
|
||||
static void __outer(void);
|
||||
|
||||
void
|
||||
outer(void)
|
||||
{
|
||||
int h, n;
|
||||
|
||||
save();
|
||||
|
||||
p2 = pop();
|
||||
p1 = pop();
|
||||
|
||||
h = tos;
|
||||
|
||||
if (car(p1) == symbol(OUTER)) {
|
||||
push_cars(cdr(p1));
|
||||
p1 = pop();
|
||||
}
|
||||
|
||||
p3 = p2;
|
||||
|
||||
if (car(p2) == symbol(OUTER))
|
||||
p2 = cadr(p2);
|
||||
|
||||
push(p1);
|
||||
push(p2);
|
||||
|
||||
if (istensor(p1) && istensor(p2))
|
||||
__outer();
|
||||
else if (isnum(p1) && istensor(p2))
|
||||
scalar_times_tensor();
|
||||
else if (istensor(p1) && isnum(p2))
|
||||
tensor_times_scalar();
|
||||
|
||||
if (car(p3) == symbol(OUTER))
|
||||
push_cars(cddr(p3));
|
||||
|
||||
n = tos - h;
|
||||
|
||||
if (n > 1) {
|
||||
list(n);
|
||||
p1 = pop();
|
||||
push_symbol(OUTER);
|
||||
push(p1);
|
||||
cons();
|
||||
}
|
||||
|
||||
restore();
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
//
|
||||
// compute the outer product to tos-1 and tos
|
||||
//
|
||||
// both args are tensors
|
||||
//
|
||||
//-----------------------------------------------------------------------------
|
||||
|
||||
static void
|
||||
__outer(void)
|
||||
{
|
||||
int i, j, k, ndim, nelem;
|
||||
|
||||
save();
|
||||
|
||||
p2 = pop();
|
||||
p1 = pop();
|
||||
|
||||
ndim = p1->u.tensor->ndim + p2->u.tensor->ndim;
|
||||
|
||||
if (ndim > MAXDIM)
|
||||
maxdim_error();
|
||||
|
||||
nelem = p1->u.tensor->nelem * p2->u.tensor->nelem;
|
||||
|
||||
p3 = alloc_tensor(nelem);
|
||||
|
||||
p3->u.tensor->ndim = ndim;
|
||||
|
||||
for (i = 0; i < p1->u.tensor->ndim; i++)
|
||||
p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
|
||||
|
||||
j = i;
|
||||
|
||||
for (i = 0; i < p2->u.tensor->ndim; i++)
|
||||
p3->u.tensor->dim[j + i] = p2->u.tensor->dim[i];
|
||||
|
||||
k = 0;
|
||||
|
||||
for (i = 0; i < p1->u.tensor->nelem; i++)
|
||||
for (j = 0; j < p2->u.tensor->nelem; j++) {
|
||||
push(p1->u.tensor->elem[i]);
|
||||
push(p2->u.tensor->elem[j]);
|
||||
multiply();
|
||||
p3->u.tensor->elem[k++] = pop();
|
||||
}
|
||||
|
||||
push(p3);
|
||||
|
||||
restore();
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
//
|
||||
// Inner product (aka dot product)
|
||||
|
|
|
|||
42
test.cpp
42
test.cpp
|
|
@ -1043,48 +1043,6 @@ char *script[] = {
|
|||
"B=quote(B)",
|
||||
"",
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
//
|
||||
// outer
|
||||
//
|
||||
//-----------------------------------------------------------------------------
|
||||
|
||||
"outer(a,b)",
|
||||
"outer(a,b)",
|
||||
|
||||
"A=quote(A)",
|
||||
"A",
|
||||
|
||||
"B=quote(B)",
|
||||
"B",
|
||||
|
||||
"C=quote(C)",
|
||||
"C",
|
||||
|
||||
"outer(C,outer(B,A))",
|
||||
"outer(C,B,A)",
|
||||
|
||||
"H33=hilbert(3)",
|
||||
"H33",
|
||||
|
||||
"H44=hilbert(4)",
|
||||
"H44",
|
||||
|
||||
"H55=hilbert(5)",
|
||||
"H55",
|
||||
|
||||
"H3344=outer(H33,H44)",
|
||||
"H3344",
|
||||
|
||||
"H4455=outer(H44,H55)",
|
||||
"H4455",
|
||||
|
||||
"H33444455=outer(H33,H44,H44,H55)",
|
||||
"H33444455",
|
||||
|
||||
"inner(H3344,H4455)-contract(H33444455,4,5)",
|
||||
"0",
|
||||
|
||||
// rank
|
||||
|
||||
"rank(A)",
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue