// This takes a simplex from each brane, chooses a
// random point in each, and returns distance^2 between them.
#include <stdlib.h>
#include <stdio.h>
#include "l.h"
R smp(R* w, int a, int b){int k = a+b+1; R*tr = w+k*(k-1);
  if(0){int i, j; for(i=0; i<k; ++i) {for(j=0; j<k; ++j) printf("%8.5f ", w[i*k+j]); printf("kt\n");}}
  R x[k]; {int j=k; while(j--) x[j] = tr[j];}
  void ss(R* a, int ac){R const f = 0x1p-31;
    if(ac) {ui sia[ac]; {int i=ac; while(i--) sia[i] = random();}
      rdx(sia, sia+ac-1, 1<<30);
      R W[ac]; {int j=ac-1; while(j--) W[j+1] = f*(sia[j+1]-sia[j]); W[0] = f*sia[0];}
      if(0){int h=ac; while(h--) printf("%6.3f ", W[h]); printf("httk\n");}
      if(1){R s=0; {int j=ac; while(j--) s+= W[j];} if(s>1) {int j=ac;
         while(j--) printf("%d %x %6.2f, ", j, sia[j], W[j]); printf("vjoe %7.3f\n", s);}}
      {int j = k; while(j--) {R s=0; {int i=ac; while(i--) s+= W[i]*a[i*k+j];} x[j] += s;}}}}
   if(0) {ss(w+a*k, b); {int j=k; while(j--) x[j]=-x[j];} ss(w, a);}
   else {ss(w, a);{int j=k; while(j--) x[j]=-x[j];} ss(w+a*k, b);}
  {int j=k; R s=0; while(j--) s+= x[j]*x[j];
  if (0) if(s>4) {int j=k; while(j--) printf("%6.2f  ", x[j] );
         pm((mat){w, k, k}, "wjot"); exit(printf("whog %e\n", s));}
  if(0) printf("jje s=%e a=%d b=%d\n", s, a, b);
  return s;}}

/*
#include <math.h>
static const int sz = 10000;
int main(){R t[][4] =
     {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
//  {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
//ave = 1.282212e+00 mx=1.686407e+00 mn=1.000017e+00
//  {{1, 0, 0}, {0, 0, 1}, {3, 0, 0}};
//ave = 2.568355e+00 mx=3.153968e+00 mn=2.003977e+00
  if(0){int k=10; while(k--) printf("%8.4f mk\n", sqrt(smp((R*) t, 2, 1)));}
  else{int k=sz; R s=0, mx=0, mn=1e9; while(k--) {R u = sqrt(smp((R*) t, 2, 1));
     s += u; mx=mx<u?u:mx; mn=mn<u?mn:u;}
     printf("ave = %e mx=%e mn=%e\n", s/sz, mx, mn);}
  return 0;}
*/