#include <imsl.h>
#include <math.h>
#define LD 50000
void rob1(mtot, wt, g, nd, dpar, q, amu, asd)
int mtot, nd; double *amu, *asd, *wt, *g, *dpar, *q;
{
   int nv1[LD], m, id, l, lstep;
   double ql, w, wg, w1, wg1, wg2, ak, t[LD], d[LD];

   for (m=0; m<=mtot-1; m++) 
      nv1[m]=m;
   imsl_d_sort(mtot,g,IMSL_RETURN_USER,g,IMSL_PERMUTATION_USER, nv1, 0);
   w1=0; wg1=0; wg2=0;
   for (m=0; m<=mtot-1; m++) {
      w=wt[nv1[m]];
      wg=w*g[m];   
      w1=w1+w;
      wg1=wg1+wg;
      wg2=wg2+wg*g[m];
      d[m]=w1;
      t[m]=wg1;
   }
   *amu=wg1/w1;
   *asd=sqrt(wg2/w1-(*amu)*(*amu));
   for (id=0; id<=nd-1; id++) {
      ak=exp(dpar[id]);
      l=mtot/2;
      lstep=mtot/2;
      do {
         lstep=(lstep-1)/2+1;
         ql=(t[l-1]+ak*(wg1-t[l-1]))/(d[l-1]+ak*(w1-d[l-1]));
         if (ql > g[l])   
            l=l+lstep;
         else if (ql < g[l-1]) 
            l=l-lstep;
      } while ((ql < g[l-1]) || (ql > g[l]));
      lstep=(lstep-1)/2+1;
      ql=(t[l-1]+ak*(wg1-t[l-1]))/(d[l-1]+ak*(w1-d[l-1]));
      q[id]=ql;
   }
   return;
}