/* 資源配分問題を解く増分法のアルゴリズム */

#include <stdio.h>
#include <math.h>
#define M 100         /* 配列Aの最大サイズ */
struct value          /* 構造体valueの定義 */
{
 float d;
 int j;
};
/* 関数の宣言 */
void resource(int n, int N);
float fn(float x, int j);
struct value deletemin(struct value *A, int n);
void downmin(int i, struct value *A, int n);
void insert(struct value dj, struct value *A, int n); 
void upmin(int i, struct value *A, int n);
void swap(int i, int j, struct value *A);

main()
/* 資源配分問題のテストプログラム */
{
 int n, N;
 n=5; N=6;      /* 変数の数nと右辺の値Nの設定 */
 resource(n, N);
}

void resource(int n, int N)
/* 増分法を用いて資源配分問題を解く. 変数の個数n, 右辺の値N. */
{
 struct value dj, min, inc[M];  /* incは増分量を貯えるヒープ */
 int j, i, x[M];                /* xは解 */
 float sum, fjx, fjx1;
 sum=0.0;
 for(j=0; j<n; j++)             /* 初期化 */
   {
    x[j]=0;                     /* 初期解 */
    fjx=fn(0.0, j); fjx1=fn(1.0, j);
    sum=sum+fjx;
    dj.d=fjx1-fjx; dj.j=j;
    insert(dj, inc, j);         /* 初期ヒープ */
   }
 for(i=1; i<=N; i++)            /* 増分法の反復 */
   {
    min=deletemin(inc, n);
    x[min.j]=x[min.j]+1;        /* 増分量最小のxjを1増加 */
    sum=sum+min.d;
    fjx=fn(x[min.j], min.j);    
    fjx1=fn(x[min.j]+1.0, min.j);
    dj.d=fjx1-fjx; dj.j=min.j;
    insert(dj, inc, n-1);       /* ヒープの更新 */
   }
 printf("Solution x =");        /* 最適解の出力 */
 for(j=0; j<n; j++) printf("%d ", x[j]);
 printf("\n");
 printf("Objective value = %f\n", sum);
 }

float fn(float x, int j)
/* 目的関数値fj(x)を与えるプログラム例 */
{
 float f;
 if(j==0) f=1.0/(x+1.0);
 if(j==1) f=(x-2.5)*(x-2.5);
 if(j==2) f=1.0+0.1*x*x*x;
 if(j==3) f=-1.0+(x-2.0)*(x-2.0);
 if(j==4) f=3.0-sqrt(x);
 if(j<0 || j>4)                 /* jの誤りチェック */
   {printf("Error: j = %d is out of the range.", j); exit();}
 return(f);
}

void insert(struct value dj, struct value *A, int n)
/* ヒープ A[0], ...,A[n-1]に対し要素 dj を挿入 */
{
 A[n]=dj;
 upmin(n, A, n+1);              /* ヒープの条件の回復 */
}

struct value deletemin(struct value *A, int n)
/* ヒープ A[0], ...,A[n-1]から最小要素を消去 */
{
 struct value min;
 min=A[0]; A[0]=A[n-1];
 if(n>1) downmin(0, A, n-1);    /* ヒープの条件の回復 */
 return(min);
}

void upmin(int i, struct value *A, int n)
/* 要素 A[i] から上へ向かってヒープ条件の回復 */
{
 int j;
 if(i<0 || i>=n)                /* 位置 i のチェック */
   { 
    printf("Illegal element i = %d for n = %d\n", i, n); 
    exit();
   }
 if(i==0) return;
 j=(i+1)/2-1;                   /* A[i] の親 A[j] */
 if(A[j].d>A[i].d)
   {
    swap(i, j, A); 
    upmin(j, A, n);             /* A[j] から上へ */
   }
 return;
}


void downmin(int i, struct value *A, int n)
/* A[i] から下へ向かってヒープ条件の回復; 最終要素は A[n-1] */
{
 int j;
 if(i<0 || i>=n)                /* 位置 i のチェック */
   { 
    printf("Illegal element i = %d for n = %d\n", i, n); 
    exit();
   }
 j=2*i+1;                       /* A[i] の子 A[j],A[j+1] */
 if(j>=n) return;
 if(j+1<n && A[j].d>A[j+1].d) j=j+1; 
 if(A[j].d<A[i].d) 
   {
    swap(i, j, A);
    downmin(j, A, n);           /* A[j] から下へ */
   }
 return;
}


void swap(int i, int j, struct value *A)
/* 構造体A[i]とA[j]の交換 */
{
 struct value temp;
 temp=A[i]; A[i]=A[j]; A[j]=temp;
 return;
}