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

typedef int* reg;
typedef int* dbl;

int  nbdits  = 0;




void addreg( reg d, reg x, reg y )
{ int ret=0;
  int i;
  for( i = 0; i < nbdits; i++) {
    ret = x[i] + y[i] + ret;
    if ( ret > 9  ) {
      d[i] = ret - 10;
      ret  = 1;
    }  
    else { d[i] = ret; ret = 0;}
  }
}

void cdiv( reg x, int c)
{ int i, ret = 0;
  for( i = nbdits - 1; i >=0; i--){
   ret = ret * 10 + x[i];
   x[i] = ret / c;
   ret =  ret - c * x[i];
 } 
}

void carre( dbl y, reg x )
{ int i, j, ret;
  for( i = 0; i < nbdits * 2; i++) y[i] = 0;
  for( i = 0; i < nbdits; i++)
    for( j = 0; j < nbdits; j++)
      y[ i + j] += x[i] * x[j];
  ret = 0;
  for( i = 0; i < nbdits * 2; i++) {
    ret += y[i];
    y[i] = ret % 10;
    ret /= 10;
  }
}

int cmpdbl( dbl x, dbl y)
{ 
  int i = 2 * nbdits - 1;
  while ( ( i >=0 )  && x[i] == y[i] ) i--;
  if ( i < 0 )       return 0;
  if ( x[i] > y[i] ) return 1;
  return -1;
}

void cpreg( reg d, reg s )
{
  int i;
  for( i = 0; i < nbdits; i++)
    d[i] = s[i];
}

reg  initreg( void )
{ reg res; 
  res = ( reg ) calloc( nbdits, sizeof(int) );
  return res;
}

dbl  initdbl( void )
{ dbl res; 
  res = ( dbl ) calloc( nbdits * 2, sizeof(int) );
  return res;
}


dbl  atodbl( char *str )
{ dbl res;
  int i, dpl;
  res = initdbl();
  dpl = strlen( str ) - 1;
  for( i = 0; i <= dpl ; i++)
    res[i] = str[ dpl - i ] - '0';
  return res;
}

dbl  randdbl( void )
{ dbl res;
  int i;
  res = initdbl();
  for( i = 0; i < 2 * nbdits - 4; i++)
    res[i] = 1 + (random() % 8);
  return res;
}


void printreg( reg x )
{ int i = nbdits - 1; 
  while ( i > 0 && x[i] == 0 ) i--;
  while ( i >=0 ) {
    printf("%d", x[i]);
    i--;
  }
}

void printregln( reg x )
{ int i = nbdits - 1;
  printf("\n");
  while ( i > 0 && x[i] == 0 ) i--;
  while ( i >=0 ) {
    printf("%d", x[i]);
    i--;
  }
}

void printdbl( dbl x )
{ int i = 2 * nbdits - 1; 
  while ( i > 0 && x[i] == 0 ) i--;
  while ( i >=0 ) {
    printf("%d", x[i]);
    i--;
  }
}

void usage( void )
{
  printf("\nusage:");
  printf("\n   -t:taille");  
  printf("\n   -h:help"); 
  printf("\n");
  exit(1);
}
void option( char *str )
{ 
  if ( *str != '-') usage(); 
  str++;
  switch( *str) {
  case 't': str++; nbdits = atoi(str); break;
  case 'h':
  default :usage();
  }
}
int main( int argc, char* argv[] )
{ 
  reg a, b, x;
  dbl y, t;
  int i, cpt;
  for( i = 1; i < argc; i++)
    option( argv[i] );
  if ( ! nbdits ) {
          nbdits = 1 + (1 + strlen( argv[1] ) ) / 2;  
          y = atodbl( argv[1] );
  } else {
    y = randdbl( );
    printdbl(y);
  }
  a = initreg(); b = initreg(); x = initreg();
  t = initdbl();
  b[ nbdits - 2 ] = 9;

  cpt = nbdits * log( 10 ) / log(2);
  printdbl( y );
  printf("\nIterations :%d", cpt);
  while ( cpt--) {
    addreg( x, a, b );
    cdiv( x, 2 ); 
    carre( t, x );
    if ( cmpdbl( t, y) > 0 ) cpreg( b, x );
    else cpreg( a, x);
  }

  printf("\nRacine:");
  printregln( a );  
  printregln( x );  
  printregln( b );
  free(a); free(b); free(x);  
  free(y); free(t);
  printf("\n");
  return 0;
}