2. projekt IZP2008

/*
 * Soubor:  proj2.c
 * Datum:   25.11.2008
 * Autor:   kolemjdouci
 * Projekt: Iteracni vypocty, projekt è. 2 pro pøedmìt IZP
 * Popis:   vyocty funkci ln a tangens a take prevadi zrychleni
 *          na polohu a naopak
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdbool.h>
#include <math.h>
#include <errno.h>
#include <limits.h>


#define absolute(val)  ( ( (val) >=0)? (val) : -(val) )
#define MIN_ACC 1e-15
#define TIME_SPAN 0.5
#define TIME_SPAN_SQ_MUL_HALF 0.125 // "(time span sqared)*1/2"
#define DIMENSION 3 // 3 rozmery nekdy nestaci :D

/** konstanty */
const double IZP_E = 2.7182818284590452354;        // e
const double IZP_PI = 3.14159265358979323846;      // pi
const double IZP_2PI = 6.28318530717958647692;     // 2*pi
const double IZP_PI_2 = 1.57079632679489661923;    // pi/2
const double IZP_PI_4 = 0.78539816339744830962;    // pi/4
//

/** hlavicky fci */
int tan_converg(double *);
int sin_converg(double *);
unsigned int nat_log_converg(double *, int *);
void printCode(int);
int signum(double);
double sinus(double, double);
double cosin(double, double);
double nat_log(double, double);
struct param getParams(int, char *[]);
//double accel2position(double, TRADAR *); jak se vlastne delaji prototypy se
//double position2accel(double, TRADAR *); strukturama v argumentu?
void math_cycle(int,double );
void rad_cykl(int);

/** struktura pro hodnoty pozice a rychlost radaru */
typedef struct radar{
  double position;
  double speed;
} TRADAR;

/** struktura pro prenos parametru shellu */
typedef struct param{
  unsigned int mode;
  double acc;
  unsigned int state;
} TPARAM;


/** kody jednotlivych modu programu */
enum tparams
{ //nazvy jsou snad pochopitelne
  RAD1 = 0,
  RAD2,
  LN,
  TAN,
  SIN,
  COS,
};

/** Kódy chyb a stavù programu */
enum tcodes
{
  EOK = 0,     /**< vse v poradku */
  ECLWRONG,    /**< Chybný pøíkazový øádek */
  EACCWRONG,   /**< Chybna presnost */
  EUNKNOWN,    /**< Neznámá chyba */
  CHELP,       /**< Tisk napovedy */
  ENDCODE,     /**< zará¾ka */
};


/** Pole textových øetìzcù vypisovaných funkcí printCode. pouzite z proj1 */
const char *CODEMSG[] =
{
  /* EOK */
  "V¹e v poøádku.\n",
  /* ECLWRONG */
  "Chybné parametry pøíkazového øádku!\n",
  /* EACCWRONG */
  "Jako presnost vypoctu muze figurovat pouze kladne realne cislo\n",
  /* EUNKNOWN */
  "Nastala nepøedvídaná chyba! Vypnìte poèítaè a rychle uteète.\n",
  /* CHELP */
  "Program Iteracni vypocty\n"
  "Autor: kolemjdouci - 2008\n"
  "Program provádí ...\n"
  "Pou¾ití: proj1 -h\n"
  "         proj1 -ln double (presnost)\n"
  "         proj1 -tan double (presnost)\n"
  "         proj1 -sin double (presnost)\n"
  "         proj1 -cos double (presnost)\n"
  "         proj1 -radar1\n"
  "         proj1 -radar2\n"
  "Popis parametrù:\n"
  "  [...]\n"
};

/**
 * "extended_strtod" - nastavba strtod, pri jakemkoliv problemu vraci NAN
 * @param *string - jaky retezec se pravadi na double
 * @return pri ne uplne korektnim vstupu vraci NAN, jinak stejne jako strtod
 */
double ext_strtod(char *string)
{
  char *endptr;
  double value;
  errno=0;
  value = strtod(string, &endptr);
  if (errno != 0 || *endptr != '\0' )
  {
    value=0.0/0.0 ;
  }
  return value ;
}


/**
 * Zpracuje argumenty pøíkazového øádku a vrátí je ve struktuøe TParams.
 * Pokud je formát argumentù chybný, ukonèí program s chybovým kódem.
 * @param argc Poèet argumentù.
 * @param argv Pole textových øetìzcù s argumenty.
 * @return Vrací analyzované argumenty pøíkazového øádku.
 */
TPARAM getParams(int argc, char *argv[])
{
  TPARAM result =
  { // inicializace struktury
    .mode = 0,
    .acc = 0,
    .state = EOK,
  };
  if (argc == 2)
  { // jeden parametr
    if (strcmp("-h", argv[1]) == 0)
    {// tisk nápovìdy
      result.state = CHELP;
    }
    else if (strcmp("-radar1", argv[1]) == 0)
    {
      result.mode = RAD1 ;
    }
    else if (strcmp("-radar2", argv[1]) == 0)
    {
      result.mode = RAD2 ;
    }
    else
    {
      result.state = ECLWRONG ;
    }
  }
  else if (argc == 3)
  { // dva parametry

    if ( (strcmp("-ln", argv[1]) == 0))
    {
      result.mode = LN ;
    }
    else if ( (strcmp("-tan", argv[1]) == 0))
    {
      result.mode = TAN ;
    }
    else if ( (strcmp("-sin", argv[1]) == 0))
    {
      result.mode = SIN ;
    }
    else if ( (strcmp("-cos", argv[1]) == 0))
    {
      result.mode = COS ;
    }
    else
    {
      result.state = ECLWRONG ;
    }
    /** kontrola zadane presnosti */
    if (result.state!=ECLWRONG)
    {
      char *str;
      double val;
      str=argv[2];
      val=ext_strtod(str);
      if (val==0.0/0.0 || val<=0 || !(isfinite(val)) )
      {
        result.state = EACCWRONG ;
      }
      else
      {
        result.acc = val ;
      }
    }
  }
  else
  { // chybný poèet parametrù
    result.state = ECLWRONG;
  }
  return result;
}

/**
 * premistuje tangens do intervalu nejlepsi konvergence,
 * skoro to samy jako sinus, jenom o jedne presun min
 * @param *x - hodnota s kterou pocitame v radianech
 * @return - vysledne znamenko
 */
int tan_converg(double *x)
{
  int sign;
  sign = signum(*x);
  *x = absolute(*x);
  if ( *x >= IZP_PI)
  {
    double help;
    help = floor(*x / IZP_PI);
    *x = *x - ( help * IZP_PI );
  }
  if ( *x >= IZP_PI_2)
  {
      *x = IZP_PI - *x;
      sign = -sign ;
  }
  if (*x == 0) // problemy se zapornou nulou
  {
      return 0;
  }
  return sign;
}

/**
 * premistuje sinus do intervalu nejlepsi konvergence
 * nejdrive premisti na interval <0;inf) (l. 253-254)
 * pak na interval <0;2pi) (l. 255-260)
 * pak na interval <0;pi) (l. 261-265)
 * a nakonec na interval <0;pi/2) (l. 266-269)
 * @param *x - hodnota s kterou pocitame v radianech
 * @return - vysledne znamenko
 */
int sin_converg(double *x)
{
  int sign;
  sign = signum(*x);
  *x = absolute(*x);
  if ( *x >= IZP_2PI)
  {
    double help;
    help = floor(*x / IZP_2PI);
    *x = *x - ( help * IZP_2PI );
  }
  if ( *x >= IZP_PI)
  {
    sign = -sign;
    *x = *x - IZP_PI ;
  }
  if ( *x >= IZP_PI_2)
  {
    *x = IZP_PI - *x;
  }
  if (*x == 0) // problemy se zapornou nulou
  {
      return 0;
  }
  return sign;
}

/**
 * funkce tangens
 * prakticky vola pouze podil sinu a cosinu
 * @param x - tangens ceho se pocita
 * @param acc - s jakou presnosti
 * @return vraci hodnotu fce tan, mimo def. obor NAN
 */
double tangens(double x, double acc)
{
  if ( !(isfinite(x)) )
  {
    return 0.0/0.0;
  }
  acc=acc/100 ;
  int sign;
  sign = tan_converg(&x);
  if (x==IZP_PI_2)
  {
    return 0.0/0.0;
  }
/* - tahle heuristika se prokazala zalostne neucinou
  else if ( x>1.5 )
  { // tan(2x)=2*tan(x)/(1-(tan(x)^2))
    double help;
    help=tangens(x/2,acc) ;
    x = (2*help)/(1-help*help) ;
  } */
  else
  {
    x = sinus(x, acc)/cosin(x, acc);
  }
  return x*sign;
}

/**
 * Vytiskne kód chyby nebo text popisující stav programu.
 * @param code - kód chyby nebo stavu programu
 */
void printCode(int code)
{
  if (code < EOK || code >= ENDCODE)
  { code = EUNKNOWN; }

  FILE *outstream = stdout;
  if (code <= EUNKNOWN)
  {
    outstream = stderr;
  }

  fprintf(outstream, "%s", CODEMSG[code]);
}

/**
 * fce signum... puvodne jsem nechtel vubec pouzivat matematickou knihovnu
 * @param x - ktereho cisla chci signum
 * @return 1/0/-1 v zavislosti na vysledku
 */
int signum(double x)
{
  if ( x == 0 )
  {
    return 0;
  }
  else if ( x > 0 )
  {
    return 1;
  }
  else
  {
    return -1;
  }
}

/**
 * radar prevadejici pozici na zrychleni
 * vzorec:
 * nova draha= puvodni draha + puvodni rychlost*ubehly cas + 0.5*zrychleni
 *             *ubehly cas na druhou
 * nova rychlost= puvodni rychlost + zrychleni*ubehly cas
 * @param accel - hodnota zrychleni, kterou prevadim
 * @param *rad -pointer na strukturu, kde ukladam prubeznou polohu a rychlost
 * @return - aktualni pozice
 */
double accel2position(double accel, TRADAR *rad)
{
  if (!(isfinite(accel)))
  {
    return 0.0/0.0 ;
  }
  rad->position=(*rad).position+(*rad).speed/2 +TIME_SPAN_SQ_MUL_HALF * accel;
  rad->speed = TIME_SPAN * accel + (*rad).speed;
  return (*rad).position;
}

/**
 * radar prevadejici z pozice na zrychleni
 * vzorec, viz. predchozi radar
 * @param position - poloha, kterou prevadim an zrcyhleni
 * @param *rad -pointer na strukturu, kde ukladam prubeznou polohu a rychlost
 * @return aktualni zrychleni
 */
double position2accel(double position, TRADAR *rad)
{ // stx = s0 + v0t + 0.5 a t^2 && vtx = v0 + at
  double accel;
  if (!(isfinite(position)))
  {
    return 0.0/0.0 ;
  }
  accel = ((position-(*rad).position) - TIME_SPAN * (*rad).speed) / TIME_SPAN_SQ_MUL_HALF ;
  rad->speed = (*rad).speed + accel * TIME_SPAN ;
  rad->position = position;
  return accel ;
}

/**
 * cyklus nacitani a tisku hodnot pro radary
 * @param mode - ktery ze dvou radaru se pouzije
 */
void rad_cykl(int mode)
{
  int i = 0;
  int j = 1;
  TRADAR radar[DIMENSION];
  radar[0].position=0.0 , radar[0].speed=0.0 ;
  while (j < DIMENSION)
  { // nulovani hodnot pro vsechny radary (prvky pole radaru)
    radar[j]=radar[0];
    j++;
  }
  double (*rad_fce)(double,TRADAR *); // pointer na pozadovanou fci
  rad_fce = (mode==RAD1)? accel2position : position2accel ;
  char buffer[128];
  double out,in;
  while ( fscanf(stdin,"%127s", buffer) != EOF )
  {
    in = ext_strtod(buffer);
    out = rad_fce(in,&radar[i]);
    printf("%.10le\n", out);
    i++;
    i %= DIMENSION;
    //i = (++i) % DIMENSION; - nechapu proc tohle nefungovalo
  }
}

/**
 * cyklus nacitani a tisku hodnot pro matematicke fce
 * @param mode - ktera funkce se pouzije
 * @param acc - presnost s jakou se bude pocitat
 */
void math_cycle(int mode,double acc)
{
  double (*mat_fce)(double, double);
  switch (mode)
  { // ktera funkce se bude pouzivat
    case LN : mat_fce = nat_log; break;
    case TAN : mat_fce = tangens; break;
    case SIN : mat_fce = sinus; break;
    case COS : mat_fce = cosin; break;
  }
  double in;
  double out;
  char buffer[128];
  while ( fscanf(stdin,"%127s", buffer) != EOF )
  {
    in = ext_strtod(buffer);
    out = mat_fce(in,acc);
    printf("%.10le\n", out);
  }
}



/**
 * vypocet fce sinus, za pomoci taylorova rozvoje pro sin(x)
 * @param x - sinus ceho se ma vypocitat
 * @param acc - (accuracy) - s jakou presnosti
 * @return hodnota fce
 */
double sinus(double x, double acc)
{
  int sign;
  if ( !(isfinite(x)) )
  {
    return 0.0/0.0;
  }
  sign = sin_converg(& x);
  /** dal je to obdobne, jako vypocet cosinu na demonstracnim cviceni */
  double sum, x_sqd, help;
  x_sqd = x * x;
  help = x;
  sum = help;
  acc = (acc< MIN_ACC && acc>0)? acc : MIN_ACC;
  int j=3;
  do {
    help *= -x_sqd/ (j * (j-1));
    sum += help;
    j += 2;
  } while( absolute(help) > acc);
  return (sum * sign);
}


/**
 * nebyt toho, ze potrebuju pointer na tuto fci,
 * stacilo by mi makro, viz. fce sinus
 */
double cosin(double x, double acc)
{
  return sinus(x + IZP_PI_2, acc );
}


/**
 * tady predelavam logaritmus do intervallu nejlepsi konvergence
 * pro mala cisla pouzivam vzorec ln(x^y)=y*ln(x) (y=-1)
 * pro vetsi nez 1.5 zase ln(a*b)=ln(a)+ln(b) (b=e)
 * @param *x - pointer na cislo, ktere mam dostat do inter. nejlepsi. konverg.
 * @param *sign - urcuje, zda se bude navratova hodnota pozdeji k vysledku
 *                pricitat, nebo od nej odecitat, viz. vzorec pro mala cisla
 * @return kolikrat tam bude "b" z druheho vzroce
 */
unsigned int nat_log_converg(double *x, int *sign)
{
  if (*x<0.55)
  {
    *x = 1.0 / *x ;
    *sign = -1 ;
  }
  int shift = 0;
  while ( *x > 1.5 ) // 1.5/e ~= 0.552
  {
    *x = *x / IZP_E ;
    shift++ ;
  }
  *x -= 1 ;
  return shift;
}


/**
 * vypocet fce ln pocitam na tomto intervalu hodnotu fce s co nejnizsim
 * poctem operaci  za pomoci taylorova rozvoje pro ln(x+1)
 * @param x - ln ceho se ma vypocitat
 * @param acc - (accuracy) - s jakou presnosti
 * @return vypocitana hodnota
 */
double nat_log(double x, double acc)
{
  if (x<0 || isnan(x) )
  { // kdy vratit NaN
    return 0.0/0.0 ;
  }
  if ( isinf(x) )
  { // kdy vratit Inf
    return 1.0/0.0 ;
  }
  if (x==0)
  { // kdy vratit -Inf
    return -(1.0/0.0) ;
  }
  unsigned int shift;
  int sign = 1;
  shift = nat_log_converg(&x, &sign);
  acc = (acc< MIN_ACC && acc>0)? acc : MIN_ACC;
  double sum, help, diff;
  sum = x ;
  help = x;
  diff = x;
  int i = 1;
  while ( absolute(diff) > acc)
  {
    i ++ ;
    help *= -x ;
    diff = help/i ;
    sum += diff ;
  }
  return sign*(sum+shift) ;
}


/////////////////////////////////////////////////////////////////
/**
 * Hlavní program.
 * nejdriv vola fci na analyzu argument cmd, pak pripadne skonci
 * podle nich zavola bud funkci na radary nebo na mat. fce
 */
int main(int argc, char *argv[])
{

  TPARAM params = getParams(argc, argv);
  if (params.state != EOK)
  { // nìco nestandardního
    printCode(params.state);
    return (params.state == CHELP)? EXIT_SUCCESS : EXIT_FAILURE;
  }
  if (params.mode>RAD2)
  {
    math_cycle(params.mode,params.acc);
  }
  else
  {
    rad_cykl(params.mode);
  }

//    printf("%.13Le\n\n", (double)( ((double)0) * ((int)(-1)) ) ); // = -0.0
    return EXIT_SUCCESS;
}