/*
  Copyright (C) 2003, Guillaume-Latapy - LIAFA (University of Paris 7)

  Disclaimer: this code is provided 'as-is', without any express or
  implied warranty. In no event will the authors be held liable for
  any damages arising from the use of this software.

  Permission is granted to anyone to use this code for any purpose,
  to alter it and redistribute it without restrictions. Just feel
  free to inform us if you plan to use it. If you face any problems
  or find any bug, please contact us using our webpage:
    http://www.liafa.jussieu.fr/~guillaume/

  --------------------------------------------------

  Various degree distribution generators.
*/

#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>

using namespace std;
using std::cout;
using std::cin;
using std::cerr;

#define UNIFORM  0
#define POWERLAW 1
#define POISSON  2

#define PI 3.141592654

void
usage(char *program_name, int status) {
  if (status == EXIT_SUCCESS)
    {
      cout << "Usage: " << endl
	   << program_name << " -u -m min -M max -n balls -s sum" << endl
	   << program_name << " -p exponent -m min -M max -n balls -s sum" << endl
	   << program_name << " -e mean -m min -M max -n balls -s sum" << endl
	   << "  Generates a random degree distribution for various laws." << endl
           << "  Graph is printed on stdout with one link (betwen top and bottom) per line." << endl
	   << "  Laws:" << endl
	   << "    -u: uniform law between m and M" << endl
	   << "    -p exponent: power law of given exponent (default: 1) between m and M" << endl
	   << "    -e mean: Poisson law of given mean (default: 1)" << endl
	   << "  Options:" << endl
           << "    -m: miminal value" << endl
           << "    -M: maximal value" << endl
           << "    -n: number of balls" << endl
           << "    -s: sum of ball values" << endl
           << "    -h, this usage" << endl
           << "Remarks:" << endl
           << "  If given sum of values is not consistent with the expected sum of values," << endl
	   << "  program might not end." << endl ;
    }
  else
    {
      cerr << "Try '" << program_name << " -h' for usage information." << endl;
    }
  exit(status);
}

void
parse_args(int argc, char **argv, int &law, int &m, int &M, int &n, int &s, float &p) {
  extern char *optarg;
  extern int optind, opterr, optopt;
  char c;

  while ((c = getopt(argc, argv, "m:M:n:s:p:e:uh")) != EOF){
    switch (c) {
    case 'h':       /* help */
      usage(argv[0], EXIT_SUCCESS);
      break;
    case 'm':       /* nb vertices */
      m=atoi(optarg);
      break;
    case 'M':       /* nb links */
      M=atoi(optarg);
      break;
    case 'n':       /* nb links */
      n=atoi(optarg);
      break;
    case 's':       /* nb links */
      s=atoi(optarg);
      break;
    case 'u' :
      law=UNIFORM;
      break;
    case 'p' :
      p=atof(optarg);
      law=POWERLAW;
      break;
    case 'e' :
      p=atof(optarg);
      law=POISSON;
      break;
    case ':':       /* missing operand */
      cerr << argv[0] << ": Option -" << optopt << " requires an operand." << endl;
      usage(argv[0], EXIT_FAILURE);
      break;
    case '?':       /* unknown option */
      cerr << argv[0] << ": Unrecognized option -" << optopt << endl;
      usage(argv[0], EXIT_FAILURE);
    }
  }
}

float gammaln(float xx) {
  double x,y , tmp, ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
			24.01409824083091,-1.231739572450155,
			0.1208650973866179e-2,-0.5395239384953e-5};
  
  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (int j=0;j<5;j++)
    ser+=cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);
}

float
poisson(float xm) {
  static float sq, alxm, g, oldm=(-1.0);
  float em,t,y;

  if(xm<12.0) {
    if(xm!=oldm) {
      oldm=xm;
      g=exp(-xm);
    }
    em=-1;
    t=1.0;
    do {
      ++em;
      t*=(random()*1.)/(RAND_MAX*1.);
    } while (t>g);
  } else {
    if(xm!=oldm) {
      oldm=xm;
      sq=sqrt(2.0*xm);
      alxm=log(xm);
      g=xm*alxm-gammaln(xm+1.0);
    }
    do {
      do {
	y=tan(PI+(random()*1.)/(RAND_MAX*1.));
	em=sq*y+xm;
      } while (em<0.0);
      em=floor(em);
      t=0.9*(1.0+y*y)*exp(em*alxm-gammaln(em+1.0)-g);
    } while ((random()*1.)/(RAND_MAX*1.) > t);
  }
  return em;
}


int
myrandom(int law, int m, int M, float p) {
  if (law==UNIFORM) {
    return random()%(M-m+1)+m;
  } else   if (law==POWERLAW) {
    double r=(random()*1.)/(RAND_MAX*1.);
    double nexp=p+1.;
    double norm=1./(pow(M*1.0,nexp)-pow(m*1.0,nexp));
    double expo=log10(r/norm+pow(m*1.0, nexp))/nexp;
    return (int)floor(pow(10.,expo));
  } else { //  case POISSON:
    int val;
    do {
      val=poisson(p);
    } while (val<m || val>M);
    return val;
  }
}

vector<int> distrib(vector<int> v) {
  int m=0;
  for (int i=0; i<v.size() ; i++)
    m=max(m, v[i]);

  vector<int> d(m+1);
  for (int i=0; i<m+1 ; i++)
    d[i]=0;
  for (int i=0; i<v.size() ; i++)
    if (v[i]>=0)
      d[v[i]]++;
  return d;
}

int
main(int argc, char **argv) {
  srandom(time(NULL));

  int m=-1, M=-1, n=-1, s=-1, law=-1;
  float p=1.;

  parse_args(argc, argv, law, m, M, n, s, p);
  if (n<0 || m<0 || M<0)
    usage(argv[0], EXIT_FAILURE);
  
  vector<int> numbers(n);

  // generates numbers
  int current_sum=0;
  for(int i=0 ; i<n ; i++) {
    numbers[i]=myrandom(law, m, M, p);
    current_sum+=numbers[i];
  }


  int t=0;

  // match sum of balls with prescribed sum
  if(s>=0) {
    while(current_sum!=s) {
      t++;
      int new_ball = random()%n;
      current_sum -= numbers[new_ball];
      numbers[new_ball] = myrandom(law, m, M, p);
      current_sum += numbers[new_ball];
    }
  }

  //print distribution
  vector<int> degree_dist = distrib(numbers);
  for(int i=0;i<degree_dist.size();i++)
    cout << i << " " << degree_dist[i] << endl;
}