/* tabauss.c 
 * 
 * table de la loi de Laplace-Gauss (loi normale N(0 ; 1))
 * produit un fichier au format .tex de LaTeX
 * 
 * compilation :
 * -------------
 * gcc -o tabauss tabauss.c -lm
 *
 * utilisation :
 * -------------
 * tabgauss -d 4 -b 4 > laplgauss.tex
 * latex laplgauss.tex
 * dvips -t a4 laplgauss.dvi -o
 * pdflatex laplgauss.tex
 * 
 *           ------------------ / ------------------
 * tabauss.c, Copyright (C) 1998 Jean-Paul Davalan.
 *
 * 29900 Concarneau (France) le : Mon Oct 19 23:33:17 GMT+0200 1998
 *
 * $Modified:
 *
 * Tabgauss is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * Tabgauss is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Tabgauss; see the file COPYING.  If not, write to
 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 * Jean-Paul Davalan <jeanpaul.dvl@wanadoo.fr>
 */

/* f(t) = 1/sqrt(2*PI) * e^(-t^2/2) en utilisant la fonction erf()
  erf(x) = 2/sqrt(pi)* integral from 0 to x of exp(-t*t) dt
  pour calculer p(x) = 1/sqrt(2*pi)* integral from -infty  to x of exp(-t*t/2) dt
  Pour x >=0 on a  p(x) = 0.5 + 1/sqrt(2*pi)* integral from 0 to x of exp(-t*t/2) dt
  le changement de variable u = t/sqrt(2) donne dt = sqrt(2) * du et les bornes
  deviennent 0 et x/sqrt(2) donc
  p(x) = 0.5 + 1/sqrt(2*pi)* integral from 0 to x/sqrt(2) of exp(-u*u) * sqrt(2) * du
       = 0.5 + 1/sqrt(pi)* integral from 0 to x/sqrt(2) of exp(-u*u)  du
       = 0.5 + erf(x/sqrt(2))
 */
       

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>

#define VERSION "version 0.0"

#define max(U,V) ((U>V)?U:V)
#define min(U,V) ((U<V)?U:V)




/* Gestion des durees */
struct tm *tm;
time_t date_begin, date_now, date_end;
char *my_date();
unsigned long cpu_begin, cpu_now, cpu_end;
unsigned long cputime();
char *cpu_past(unsigned long d);
/* Options */
int     opt_;

static char *helparr[] =
{
   "Options :",
   "  --help         -h",
   "  --version      -v",
   "  --copyright    -c",
   "  --precision=n  -p n  (n positif,  default: 4)",
   "  --decimals=n   -d n (same as precision)",
   "  --bound=n      -b n (default 5)",
   "  --maximum=n    -m n (same as bound)",
   0
};

static const struct option long_options[] =
{
   {"help", no_argument, 0, 'h'},
   {"version", no_argument, 0, 'v'},
   {"copyright", no_argument, 0, 'c'},
   {"precision", required_argument, 0, 'p'},
   {"decimals", required_argument, 0, 'd'},
   {"bound", required_argument, 0, 'b'},
   {"maximum",required_argument, 0, 'm'},
   {0, 0, 0}
};

void    options();
void    help();
void    copyright();
void    init_rand(void);
void    init();
void    printtab();
void    mainloop();
void    end_application();
double  my_rand(void);

int decimals = 4;
float limit = 5.0;


int
main(int argc, char *argv[])
{
  options(argc, argv);
  init(argc, argv);
  mainloop(argc, argv);
  end_application(argc, argv);
  return 0;
}

void
options(int argc, char *argv[])
{
  int     c, longoptindex = 0, errflag = 0;
/*  if (argc == 1) {
    help(argc, argv);
    exit(2);
  }
  */
  while ((c = getopt_long_only(argc, argv,
           "hvci:o:p:d:b:m:",
           long_options, &longoptindex)) != EOF) {
    switch (c) {
      case 0:
        switch (longoptindex) {
          case 1:
            break;
          default:
            break;
        }
        break;
      case 'h':
        help(argc, argv);
        exit(1);
        break;
      case 'v':
        printf("%s\n", VERSION);
        break;
      case 'c':
        copyright(argc, argv);
        exit(1);
        break;
      case 'd':
      case 'p':
        if(optarg) {
          decimals = atoi(optarg);
          if(decimals < 0 ) {
            fprintf(stderr, "precision may be positive\n");
            exit(1);
          }
        }
        break;
      case 'b':
      case 'm':
        if(optarg){
          limit = atof(optarg);
          if(limit <= 0) {
            fprintf(stderr, "%s may be positive\n",optarg);
            exit(1);
          }
        }
        break;
      default:
        errflag++;
        break;
    }
  }
  if (errflag) {
    help(argc, argv);
    exit(2);
  }
}
/*---------------------------------------------------------------*/
unsigned long cputime ()
{
  return ((clock ()* 1000) / CLOCKS_PER_SEC);
}
/*---------------------------------------------------------------*/
void
help(int argc, char *argv[])
{
  char  **pt;
  fprintf(stderr, "usage: %s options\n", argv[0]);
  for (pt = helparr; *pt; pt++)
    fprintf(stderr, "%s\n", *pt);
}
/*---------------------------------------------------------------*/
void
copyright(int argc, char *argv[])
{
 fprintf(stdout,
    " Application %s (Gauss), Copyright (C) 1998 Jean-Paul Davalan.\n\n"
    " 29900 Concarneau (France) le : Mon Oct 19 23:33:17 GMT+0200 1998\n\n"
    " %s (Gauss) is free software; you can redistribute it and/or modify\n"
    " it under the terms of the GNU General Public License as published by\n"
    " the Free Software Foundation; either version 2, or (at your option)\n"
    " any later version.\n\n"
    " %s (Gauss) is distributed in the hope that it will be useful,\n"
    " but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
    " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
    " GNU General Public License for more details.\n"
    " You should have received a copy of the GNU General Public License\n"
    " along with %s (Gauss); see the file COPYING.  If not, write to\n"
    " the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.\n\n"
    " Jean-Paul Davalan <jeanpaul.davalan@wanadoo.fr>\n"
    ,argv[0],argv[0],argv[0],argv[0]);
}

/*---------------------------------------------------------------*/
void
mainloop(int argc, char *argv[])
{
  printtab() ;

}
/*---------------------------------------------------------------*/
void
init(int argc, char *argv[])
{
  char *s;
  date_begin = time(NULL);
  cpu_begin = cputime();
  fprintf(stdout,"%% %s running\n",argv[0]);
  s = my_date(date_begin);
  fprintf(stdout,"%% %s\n",s);
  free(s);
}
/*---------------------------------------------------------------*/
void
end_application(int argc, char *argv[])
{
  char *s;
  unsigned long cpudiff;
  cpudiff = (cpu_end = cputime()) - cpu_begin;
  s = cpu_past(cpudiff);
  fprintf(stdout,"%% %s done, cpu time = %s\n",argv[0],s);
}
/*---------------------------------------------------------------*/
char *
cpu_past(unsigned long d)
{
  unsigned long tmp;
  /* unsigned long mo = 2592000000L; */
  char chn[100];
  char * str = (char *)malloc(100 * sizeof(char));
  *str = '\x0';
/*
  if (d> 31104000000) {
    tmp = d/31104000000;
    d = d%31104000000;
    chn[0] = '\x0';
    sprintf(chn,"%ld y ",tmp);
    strcat(str,chn);
  }
  if (d>(unsigned long) mo) {
    tmp = d/ mo;
    d = d%mo;
    chn[0] = '\x0';
    sprintf(chn,"%ld m ",tmp);
    strcat(str,chn);
  }
*/
  if (d>86400000 ) {
    tmp = d/86400000;
    d = d%86400000;
    chn[0] = '\x0';
    sprintf(chn,"%ld d ",tmp);
    strcat(str,chn);
  }
  if (d>3600000 ) {
    tmp = d/3600000;
    d = d%3600000;
    chn[0] = '\x0';
    sprintf(chn,"%ld h ",tmp);
    strcat(str,chn);
  }
  if (d>60000 ) {
    tmp = d/60000;
    d = d%60000;
    chn[0] = '\x0';
    sprintf(chn,"%ld min ",tmp);
    strcat(str,chn);
  }
  if (d>1000 ) {
    tmp = d/1000;
    d = d%1000;
    chn[0] = '\x0';
    sprintf(chn,"%ld s ",tmp);
    strcat(str,chn);
  }
  if (d>0 ) {
    tmp = d;
    chn[0] = '\x0';
    sprintf(chn,"%ld ms ",tmp);
    strcat(str,chn);
  }
  return str;
}
/*---------------------------------------------------------------*/
char *
my_date(time_t t)
{
  struct tm *tm;
  char *str, chn[100];
  chn[0] = '\x0';
  tm = gmtime(&t);
  strftime(chn,60,"%A %B %d - %H h %M min %S s GMT - %Y", tm);
  str = (char *)malloc((strlen(chn)+1)*sizeof(char));
  str = strcpy(str,chn);
  return str;
}
/*---------------------------------------------------------------*/

void
init_rand(void)
{
  srand((unsigned int) (time(NULL) * getpid()));
}

/*---------------------------------------------------------------
  retourne une valeur entre 0 et 1
*/

double
my_rand()
{
  return ((double) rand()) / (RAND_MAX + 1.0);
}

/*---------------------------------------------------------------*/
char  formatstring[20]="";

void printtab(void)
{
double i, j, v, rac2 = sqrt(2)/2;
int k, r;
double bound = 10.0*limit;
sprintf(formatstring,"& $%%.%df$ \x0",decimals);
fprintf(stdout,
  "\\documentclass[10pt,a4,openright,color]{report}\n"
  "\\usepackage{epic,eepic}\n"
  "\\usepackage{aeguill}\n"
  "\\usepackage{times}\n"
  "\\usepackage{fancybox,psfig,array}\n"
  "\\usepackage[sumlimits,namelimits]{amsmath}\n"
  "\\usepackage{amsfonts,amsxtra,eucal,multicol,indentfirst}\n"
  "\\usepackage[latin1]{inputenc}\n"
  "\\usepackage[francais]{babel}\n"
  "\\usepackage[dvips]{graphicx}\n"
  "\\usepackage[dvips]{color}\n"
  "\\voffset -0.9in\n"
  "\\hoffset -1in\n"
  "\\oddsidemargin=0.4cm\n"
  "\\evensidemargin=0.4cm\n"
  "\\topmargin=1.0cm\n"
  "\\textwidth=20.2cm\n"
  "\\textheight=25.0cm\n"
  "\\footskip=1.5cm\n"
  "%%\\makeindex\n"
  "\\begin{document}\n"
/* ---------------------------------------------------------------------------------
  "\\begin{center}\n"
  "\\begin{tabular}{| c | c | c | c | c | c | c | c | c | c | c |}\n"
  "\\hline\n"
  " t & $0.00$ & $0.01$ & $0.02$ & $0.03$ & $0.04$ & "
  "$0.05$ & $0.06$ & $0.07$ & $0.08$ & $0.09$  \\\\\n"
  "\\hline"
  "\n");
  for(i = 0, r = 0; r < bound ; i += 0.1, r++) {
    fprintf(stdout,"$%.1f$ ",i);
    for(j = 0, k = 0 ; k<= 9 ; k++, j += 0.01){
      v = 0.5 + erf((i+j)*rac2);
      fprintf(stdout,formatstring , v;
    }
    fprintf(stdout,"\\\\\n");
    if(r%10 == 9) {
//    fprintf(stdout,"&&&&&&&&&&&\n");
      fprintf(stdout,"\\hline\n");
    }
  }

 fprintf(stdout,
  "\\hline\n"
  "\\end{tabular}\n"
  "\\end{center}\n"
 --------------------------------------------------------------------------------- */
  "\\begin{center}\n"
  "\\fbox{\\bf Intégrale $\\Pi(t)$ de la Loi Normale Centrée Réduite ${\\cal N}(0;\\,1)$.}\\\\\n"
  "\n"
  "$$\\Pi(t)=P(X \\leq t) = \\int_{-\\infty}^t\\frac{1}{\\sqrt{2\\pi}}\\ e^{-\\frac{x^2}{2}}dx"
  "\\ \\ \\mbox{ et }\\ \\ \\Pi(-t)=1 - \\Pi(t).$$\n"
  "\\begin{tabular}{| c | c | c | c | c | c | c | c | c | c | c |}\n"
  "\\hline\n"
  " t & $0.00$ & $0.01$ & $0.02$ & $0.03$ & $0.04$ & "
  "$0.05$ & $0.06$ & $0.07$ & $0.08$ & $0.09$  \\\\\n"
  "\\hline"
  "\n");
  for(i = 0, r = 0; r < bound ; i += 0.1, r++) {
    fprintf(stdout,"$%.1f$ ",i);
    for(j = 0, k = 0 ; k<= 9 ; k++, j += 0.01)
      fprintf(stdout,formatstring ,0.5 + erf((i+j)*rac2)/2);
    fprintf(stdout,"\\\\\n");
    if(r%10 == 9) {
//    fprintf(stdout,"&&&&&&&&&&&\n");
      fprintf(stdout,"\\hline\n");
    }
  }
 fprintf(stdout,
  "\\hline\n"
  "\\end{tabular}\n"
  "\\end{center}\n"

  "\\end{document}\n" );
}