/************************************************/
/*    Tomas Gudauskas ir Jolita Svambaryte      */
/*    VDU MIF 1  (Informatika, Matematika       */
/*                                              */
/*  Nuotakos uzdavinys - kompiuterio pirkimas   */
/*                                              */
/* Losimu teorija      Dest. Jonas Mockus       */
/*                                              */
/*   1997 Gruodis  VDU  Kaunas                  */
/*                                              */
/* Kompiliuota vytautas.vdu.lt                  */
/*     gcc nuotaka.c -o bride -lm               */
/*                                              */
/************************************************/  



#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const double pi = 3.14159265359;
/* Nepavyko padaryti nepriklausomus masyvus nuo konstantes */
/* todel apsibreziu maximalu iteraciju kieki */
const maxIter = 150;

double pw(double, double, double);
double psw(double, double, double);
double pws(double, double, double[], int, double, double, double);
double uns(double, double[], int, double, double, double);

/*****************************************************************
  PAGRINDINE PROGRAMA
******************************************************************/

int main(void)
{
   /* masyvai apibreziami galimam max kiekiui */
   double wix[maxIter], six[maxIter];
   double skr[maxIter];
   double un1s[maxIter];
   int dn1[maxIter];

   int N;
   double miu, sigma0, sigma, c; 	
   double aaa, bbb, ccc;
   double pradzia, galas, K, h;
   int KK;
   int i, a, b, x;
   double un;
   double sumA, sumB;
/*****************************************************************
   DUOMENU IVEDIMAS 
******************************************************************/
   printf("Hello, this is BRYDE's problem solution!!!\n");
   printf("Co-workers: Tomas Gudauskas and Jolita Svambaryte\n\n");
   printf("Enter the quantity of grooms N=");
   scanf("%d", &N);
   printf("Enter the average miu=");
   scanf("%lf", &miu);
   printf("Enter the scatter  sigma0=");
   scanf("%lf", &sigma0);
   printf("Enter the brides slip  sigma=");
   scanf("%lf", &sigma);
   printf("Enter the waiting looses   c=");
   scanf("%lf", &c);

/****************************************************************
   PRADINIU REIKSMIU UZSIDAVIMAS
****************************************************************/
   h = 0.1;
   K = 3*sigma0;
   pradzia = miu - K;
   galas = miu + K;
/* iteraciju kiekis */
   KK = ((galas-pradzia)/h);

/* apsidrauskime nuo masyvu perpildymo */
   if (KK>maxIter)
   {
      h = (galas-pradzia)/maxIter;
      KK = (galas-pradzia)/h;
   }

   printf("Gerumu intervalas [%lf, %lf], zingsnelis %lf, iteraciju kiekis %d \n", pradzia, galas, h, KK); 

   for (i=0; i<=KK; i++)
   {
      wix[i] = i*h+pradzia;
      six[i] = wix[i];
      /*printf(" %d  wix=%lf  \n", i, wix[i]);*/
   }
   printf("Skaiciuojame....\n");

/*****************************************************************
   SKAICIUOJAM
******************************************************************/
  printf("   \n");

   skr[N] = wix[0];  /* kritinis = blogiausiam */

   sumA = 0;
   for (a=0; a<=KK; a++)
   {
      sumB = 0;
      for (b=0; b<=KK; b++)
      {
         sumB = sumB + pw(wix[b], sigma0, miu) * psw(six[a], wix[b], sigma);
      }
/*      printf("SumB=%lf   ", sumB);*/
      sumA = sumA + uns(six[a], wix, KK, sigma0, sigma, miu) * sumB / (KK+1);
/*      printf(" sumA=%lf uns=%lf dalyba = %lf \n", sumA, uns(six[a], wix, KK, sigma0, sigma, miu), (1/ (KK+1))); */
   }
   un = sumA/ (KK+1);
/*   printf("\n  un=%lf.\n", un); */


/************************************/

for (x=N-1; x>0; x--)
{

/*************************************
iteracijos pradzia
*************************************/

/*   skr[N-1] = wix[0]; /* imam pradine reiksme */

/**** irasom kritini ******/
   i = 0;
   while ((i<=KK) && (uns(six[i], wix, KK, sigma0, sigma, miu) < un)) 
   {
      i++;
   }
   skr[x] = wix[i]; /* buvo N-1 */

/*   printf(" %d Kritinis yra %d -as ir lygus %lf.\n", x, i, skr[N-1]); */

/***** susirasome d(N-x) ********/
   for (i=0; i<=KK; i++)
   {
      dn1[i] = (uns(six[i], wix, KK, sigma0, sigma, miu) < un) ? 0 : 1;
   }

/***** suskaiciuojame u(N-x) nuo s *******/
   for (i=0; i<=KK; i++)
   {
      un1s[i] =  dn1[i] * uns(six[i], wix, KK, sigma0, sigma, miu) + (1-dn1[i]) * (un - c * (N-x) );
   }

/***** suskaiciuojam u(N-x)  ********/
   sumA = 0;
   for (a=0; a<=KK; a++)
   {
      sumB = 0;
      for (b=0; b<=KK; b++)
      {
         sumB = sumB + pw(wix[b], sigma0, miu) * psw(six[a], wix[b], sigma);
      }
/* printf("SumA=%lf, SumB=%lf, un1s[a]=%lf\n", sumA, sumB, un1s[a]); */
      sumA = sumA + un1s[a] * sumB / (KK+1);
   }
/* printf("sumA=%lf\n", sumA); */
   un = sumA/ (KK+1);
/*   printf("\n  un=%lf.\n", un); */

/*******************************/

} 

/*********************************
iteracijos pabaiga
*********************************/

printf("*********************************************************\n");
printf("*       REZULTATU (gerumu) LENTELE                      *\n");
printf("*********************************************************\n");

for (i=1; i<=N; i++)
{
   printf("* Uz %d-o\t teka jei geresnis nei %lf\t*\n", i, skr[i]);
} 

printf("*********************************************************\n");
printf("\nTai ivyko prie tokiu duomenu: \n");
   printf("Ivedete N=%d, miu=%f, sigma0=%f, sigma=%f, c=%f\n", N, miu, sigma0, sigma, c);

   return 0;
}

/**********************************************************************
      FUNKCIJOS
**********************************************************************/
double pw(double ww, double sigma0, double miu)
{
   double aa;
   aa = (1/(sqrt(2*pi*sigma0)))*exp(-0.5*((ww-miu)/sigma0)*((ww-miu)/sigma0));
   return aa;   
}

double psw(double ss, double ww, double sigma)
{
   double aa;
   aa = exp(-0.5*((ss-ww)/sigma)*((ss-ww)/sigma)) / (sqrt(2*pi*sigma));
   return aa;
}
 
double pws(double ww, double ss, double wix[], int KK, double sigma0, double sigma, double miu)
{
   double ats, sum;
   int b;
   sum = 0;
   for (b=0; b<KK; b++)
   {
      sum = sum + pw(wix[b], sigma0, miu) * psw(ss, wix[b], sigma);
   }
   ats = (pw(ww, sigma0, miu) * psw(ss, ww, sigma)) / (sum/KK);
   return ats;
}

double uns(double ss, double wix[], int KK, double sigma0, double sigma, double miu)
{
   double ats, sum;
   int b;
   sum = 0;
   for (b=0; b<KK; b++)
   {
      sum = sum + wix[b] * pws(wix[b], ss, wix, KK, sigma0, sigma, miu);
   }
   ats = sum/KK; 
   return ats;
}