// Simplicial Array
#include <stdio.h>
#include "l.h"
void rep(R* w){int j; for(j=0; j<4; ++j) printf("%6.3f ", w[j]); printf("\n");}
void r(R* w, char* m){printf("ws %s", m); rep(w);}

void sma(mat g, int e, dp rout, R* x){
int dim = g.m,  k = g.n; if(1) printf("k=%d dim=%d\n", k, dim);
R del[dim][k];
{R b = 1./e; int i=dim; while(i--) {int j=k;
   while(j--)  del[i][j] = b*g.a[i*k+j];}}
if(0){rep(del[0]); rep(del[1]);}
// inner recursive routine
void sm(int l, R * x, int e){
   R y[k];
   printf("l=%d e=%d\n", l, e);
   {int j=k; while(j--) y[j] = x[j];}
   if(l+1<dim)
    {int j=e; while(j--) {sm(l+1, y, e--); {int i=k; while(i--) y[i] += del[l][i];}}}
    else {int j=e; while(j--) {rout(y); {int i=k; while(i--) y[i] += del[l][i];}}}}
sm(0, x, e);}

int main(){
if(0){R d[][3] = {{-3, 3, 0}, {-3, 0, 3}};
// Consider a 3 by 3 by 3 cube 0 < x, y, z < 3.
// The points in the place x+y+z = 3 span an equilateral triangle, our simplex.
// We take the vertex (3, 0, 0) as bc origin.
R wh[] = {3, 0, 0};
sma((mat){(R*)d, 2, 3}, 2, rep, wh);// -> 1, 1, 1 the center.
}
else {R d[][4] = {{-4, 4, 0, 0}, {-4, 0, 4, 0}, {-4, 0, 0, 4}};
R wh[] = {4, 0, 0, 0};
sma((mat){(R*)d, 3, 4}, 4, rep, wh);// -> 1, 1, 1 the center.
}
return 0;}