sin cos tan arcsin arccos arctan (fonction infinie vers fonction convergente) C/C++

[sylvainmahe]

Bonjour

J'ai reproduis quelques fonctions de la classe math.h en langage C, pour essayer de comprendre comment cela fonctionnait, mais j'aimerais optimiser mes fonctions et surtout rendre ma fonction tan, arcsin, arccos et arctan plus précise sur certaines valeur, je m'explique:

Voila mon approximation de sin (précise sur toutes valeurs):

alpha = angle * 0.0174532925199432;

alpha2 = alpha * alpha;
alpha3 = alpha2 * alpha;
alpha5 = alpha2 * alpha3;

sinus = alpha - (0.1666666666666666 * alpha3) + (0.0083333333333333 * alpha5) - (0.0001984126984126 * (alpha5 * alpha2));

Mon approximation de cos (précise sur toutes valeurs):

alpha = angle * 0.0174532925199432;

alpha2 = alpha * alpha;
alpha4 = alpha2 * alpha2;

cosinus = 1.0 - ((0.5 * alpha2) - (0.0416666666666666 * alpha4) + (0.0013888888888888 * (alpha2 * alpha4)) - (0.0000248015873015 * (alpha4 * alpha4)));

Une fois dans des fonctions, ces calculs s’exécutent aussi rapidement que celles de math.h, donc c'est très bien

On serait tenté pour calculer tan, d'écrire:

sinus / cosinus

C'est bien mais ça prend un peu moins du double de temps à calculer par rapport à math.h, donc j'ai cherché sur wikipedia (ou autre), et j'ai trouvé ceci:

alpha = angle * 0.0174532925199432;

tangente = alpha + ((0.3333333333333333 * (alpha * alpha * alpha)) + (0.1333333333333333 * (alpha * alpha * alpha * alpha * alpha)) + (0.0539682539682539 * (alpha * alpha * alpha * alpha * alpha * alpha * alpha)));

Explication:


La c'est très très rapide à calculer, sauf que pour des valeurs d'angle supérieurs à 45 degrés, les résultats se dégradent très rapidement en précision car il faut de plus en plus d'itérations (fonction infinie) pour obtenir une précision acceptable.

Normalement on parle de séries de Taylor, pour sinus on voit qu'il faut très peut d'itérations pour avoir une précision assez grande (pour de la robotique), de même que cosinus, par contre ce n'est pas du tout le cas pour tangente.

Pour sinus et cosinus le calcul converge, il est constitué de - + - + - etc..., par contre pour calculer tangente ce sont des + + + +.


Je veux bien calculer tangente avec une division de sinus / cosinus, elle prend environ 200 micro seconde de plus que math.h, c'est pas trop grave. Par contre le soucis est le même pour arcsin, arccos, et arctan.

Jusqu'à aujourd'hui je calculais arcsin, arccos, et arctan en faisant des "tests" de différentes valeurs de sinus, cosinus, et tangente, ça marche bien, mais niveau temps de calcul on entre dans plusieurs milli (pas micro) secondes. Donc quand on veut faire un algorithme qui réagit dans la milliseconde on est foutu.

La solution est d'utiliser des séries:
arcsin:

return (SINUS + ((0.5 * ((SINUS * SINUS * SINUS) * 0.3333333333333333)) + (0.375 * ((SINUS * SINUS * SINUS * SINUS * SINUS) * 0.2)) + (0.3125 * ((SINUS * SINUS * SINUS * SINUS * SINUS * SINUS * SINUS) * 0.1428571428571428)))) * _RADIAN_TO_DEGREE;

arccos:

return (1.5707963267948966 - (COSINUS + ((0.5 * ((COSINUS * COSINUS * COSINUS) * 0.3333333333333333)) + (0.375 * ((COSINUS * COSINUS * COSINUS * COSINUS * COSINUS) * 0.2)) + (0.3125 * ((COSINUS * COSINUS * COSINUS * COSINUS * COSINUS * COSINUS * COSINUS) * 0.1428571428571428))))) * _RADIAN_TO_DEGREE;

arctan:

return (TANGENT - (((TANGENT * TANGENT * TANGENT) * 0.3333333333333333) + ((TANGENT * TANGENT * TANGENT * TANGENT * TANGENT) * 0.2) - ((TANGENT * TANGENT * TANGENT * TANGENT * TANGENT * TANGENT * TANGENT) * 0.1428571428571428))) * _RADIAN_TO_DEGREE;

Mais le phénomène étant bien qu'il faille itérer un grand nombre de fois pour certaines valeurs, pour obtenir une précision acceptable.

Auriez vous une idée de comment faire? (mise à part la blague: utilises math.h)
-----

[sylvainmahe]

Pour la fonction tan il faudrait faire un calcul de Cordic: http://www.trigofacile.com/maths/tri...dic/cordic.htm

Mais je ne comprend pas vraiment comment l'écrire en C... ?

[jiherve]

Bonsoir,
Il est trés facile de porter le CORDIC en C car cela ne demande qu'un tableau de constantes précalculées et un petit bout de code pour la gestion qui est une boucle et un test. L'avantage du CORDIC c'est qu'il est possible de calculer plusieurs fonction avec le même noyau.
On doit trouver le code tout fait sur le net.
https://people.sc.fsu.edu/~jburkardt...ic/cordic.html
J'ai une version en assembleur ATMEGA si cela t’intéresse
JR

[sylvainmahe]

Super merci

Je ne connaissais pas ce lien et il y a pleins d'exemples, je vais étudier ça demain.

Pour ton exemple en assembleur, franchement je veux bien, même si je ne m'en servirais pas. J'avais commencé à coder quelques trucs en assembleur pour le 328p, c'est vraiment instructif c'est le moins que l'on puisse dire

Après comme tu dis "très facile", je dirais qu'il y a très facile et très facile (pour un pro), moi quand je vois la facilité je vois qu'on arrive à résoudre sinus en une ligne. Par contre un arcsinus à l'air de demander déjà plus de calculs. Mais c'est à voir, peut être que quand j'aurais bien épluché ce que tu m'a donné ça fera comme d'habitude, les trucs en apparence assez complexes deviendrons faciles...

[A voir en vidéo sur Futura]

[sylvainmahe]

Ok, donc à la lecture du lien que tu m'as donné, et après avoir joué à ce que j'appelle maintenant "le jeu du chat et de la souris" en programmation, ou comment s'amuser à résoudre les dépendances des fonctions les unes avec les autres, voila ce que ça donne (la fonction arcsinus est tout en bas ):

int i4_mach ( int i )
{
  int value;

  if ( i < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "I4_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 16.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0;
    exit ( 1 );
  }
  else if ( i == 1 )
  {
    value = 5;
  }
  else if ( i == 2 )
  {
    value = 6;
  }
  else if ( i == 3 )
  {
    value = 7;
  }
  else if ( i == 4 )
  {
    value = 6;
  }
  else if ( i == 5 )
  {
    value = 32;
  }
  else if ( i == 6 )
  {
    value = 4;
  }
  else if ( i == 7 )
  {
    value = 2;
  }
  else if ( i == 8 )
  {
    value = 31;
  }
  else if ( i == 9 )
  {
    value = 2147483647;
  }
  else if ( i == 10 )
  {
    value = 2;
  }
  else if ( i == 11 )
  {
    value = 24;
  }
  else if ( i == 12 )
  {
    value = -125;
  }
  else if ( i == 13 )
  {
    value = 128;
  }
  else if ( i == 14 )
  {
    value = 53;
  }
  else if ( i == 15 )
  {
    value = -1021;
  }
  else if ( i == 16 )
  {
    value = 1024;
  }
  else if ( 16 < i )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "I4_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 16.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0;
    exit ( 1 );
  }

  return value;
}

float r4_pak ( float y, int n )
{
  static float aln210 = 3.321928094887362E+00;
  float aln2b;
  static int nmax = 0;
  static int nmin = 0;
  int nsum;
  int ny;
  float value;

  if ( nmin == 0 )
  {
    aln2b = 1.0E+00;
    if ( i4_mach ( 10 ) != 2 )
    {
      aln2b = r4_mach ( 5 ) * aln210;
    }
    nmin = aln2b * ( float ) ( i4_mach ( 12 ) );
    nmax = aln2b * ( float ) ( i4_mach ( 13 ) );
  }

  r4_upak ( y, &value, &ny );

  nsum = n + ny;

  if ( nsum < nmin )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_PAK - Warning!\n" );
    fprintf ( stderr, "  Packed number underflows.\n" );
    value = 0.0E+00;
    return value;
  }

  if ( nmax < nsum )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_PAK - Fatal error!\n" );
    fprintf ( stderr, "  Packed number overflows.\n" );
    exit ( 1 );
  }

  while ( nsum < 0 )
  {
    value = 0.5E+00 * value;
    nsum = nsum + 1;
  }

  while ( 0 < nsum )
  {
    value = 2.0E+00 * value;
    nsum = nsum - 1;
  }

  return value;
}

void r4_upak ( float x, float *y, int *n )
{
  float absx;

  absx = r4_abs ( x );
  *n = 0;
  *y = 0.0;

  if ( x == 0.0 )
  {
    return;
  }

  while ( absx < 0.5 )
  {
    *n = *n - 1;
    absx = absx * 2.0;
  }

  while ( 1.0 <= absx )
  {
    *n = *n + 1;
    absx = absx * 0.5;
  }

  if ( x < 0.0 )
  {
    *y = - absx;
  }
  else
  {
    *y = + absx;
  }

  return;
}

float r4_log ( float x )
{
  static float aln2 = 0.068147180559945309E+00;
  static float alncen[5] = {
    0.0E+00,
    +0.223143551314209755E+00,
    +0.405465108108164381E+00,
    +0.559615787935422686E+00,
    +0.693147180559945309E+00 };
  static float alncs[6] = {
    1.3347199877973882E+00,
    0.000693756283284112E+00,
    0.000000429340390204E+00,
    0.000000000289338477E+00,
    0.000000000000205125E+00,
    0.000000000000000150E+00 };
  static float center[4] = {
    1.0E+00,
    1.25E+00,
    1.50E+00,
    1.75E+00 };
  int n;
  static int nterms = 0;
  int ntrval;
  float t;
  float t2;
  float value;
  float xn;
  float y;

  if ( nterms == 0 )
  {
    nterms = r4_inits ( alncs, 6, 28.9E+00 * r4_mach ( 3 ) );
  }

  if ( x <= 0.0E+00 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_LOG - Fatal error!\n" );
    fprintf ( stderr, "  X <= 0.0\n" );
    exit ( 1 );
  }

  r4_upak ( x, &y, &n );

  xn = ( float ) ( n - 1 );
  y = 2.0E+00 * y;
  ntrval = ( int ) ( 4.0E+00 * y - 2.5E+00 );

  if ( ntrval == 5 )
  {
    t = ( ( y - 1.0E+00 ) - 1.0E+00 ) / ( y + 2.0E+00 );
  }
  else if ( ntrval < 5 )
  {
    t = ( y - center[ntrval-1] ) / ( y + center[ntrval-1] );
  }

  t2 = t * t;

  value = 0.625E+00 * xn + ( aln2 * xn + alncen[ntrval-1] 
    + 2.0E+00 * t 
    + t * t2 * r4_csevl ( 578.0E+00 * t2 - 1.0, alncs, nterms ) );

  return value;
}

float r4_mach ( int i )
{
  float value;

  if ( i < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 5.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0.0;
    exit ( 1 );
  }
  else if ( i == 1 )
  {
    value = 1.1754944E-38;
  }
  else if ( i == 2 )
  {
    value = 3.4028235E+38;
  }
  else if ( i == 3 )
  {
    value = 5.9604645E-08;
  }
  else if ( i == 4 )
  {
    value = 1.1920929E-07;
  }
  else if ( i == 5 )
  {
    value = 0.3010300E+00;
  }
  else if ( 5 < i )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 5.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0.0;
    exit ( 1 );
  }
 
  return value;
}

int r4_inits ( float dos[], int nos, float eta )
{
  float err;
  int i;
  int value;

  if ( nos < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_INITS - Fatal error!\n" );
    fprintf ( stderr, "  Number of coefficients < 1.\n" );
    exit ( 1 );
  }

  err = 0.0;

  for ( i = nos - 1; 0 <= i; i-- )
  {
    err = err + r4_abs ( dos[i] );
    if ( eta < err )
    {
      value = i + 1;
      return value;
    }
  }

  value = i;
  fprintf ( stderr, "\n" );
  fprintf ( stderr, "R4_INITS - Warning!\n" );
  fprintf ( stderr, "  ETA may be too small.\n" );

  return value;
}

float r4_sqrt ( float x )
{
  int irem;
  int iter;
  int ixpnt;
  int n;
  static int niter = 0;
  static float sqrt2[3] = {
     0.70710678118654752E+00,
     1.0E+00,
     1.41421356237309505E+00 };
  float value;
  float y;

  if ( niter == 0 )
  {
    niter = 1.443E+00 * r4_log ( - 0.104E+00 
      * r4_log ( 0.1E+00 * r4_mach ( 3 ) ) ) + 1.0E+00;
  }

  if ( x < 0.0E+00 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_SQRT - Fatal error!\n" );
    fprintf ( stderr, "  X < 0.0\n" );
    exit ( 1 );
  }
  else if ( x == 0.0E+00 )
  {
    value = 0.0E+00;
  }
  else
  {
    r4_upak ( x, &y, &n );
    ixpnt = n / 2;
    irem = n - 2 * ixpnt + 2;
    value = 0.261599E+00 + y * ( 1.114292E+00 
      + y * ( -0.516888E+00 + y * 0.141067E+00 ) );

    for ( iter = 1; iter <= niter; iter++ )
    {
      value = value + 0.5E+00 * ( y - value * value ) / value;
    }
    value = r4_pak ( sqrt2[irem-1] * value, ixpnt );
  }
  return value;
}

float r4_abs ( float x )
{
  float value;

  if ( 0.0 <= x )
  {
    value = x;
  } 
  else
  {
    value = -x;
  }
  return value;
}

 float r4_csevl ( float x, float a[], int n )
{
  float b0;
  float b1;
  float b2;
  int i;
  float twox;
  float value;

  if ( n < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms <= 0.\n" );
    exit ( 1 );
  }

  if ( 1000 < n )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms greater than 1000.\n" );
    exit ( 1 );
 }

  if ( x < -1.1 || 1.1 < x )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R4_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  X outside (-1,+1).\n" );
    exit ( 1 );
  }

  twox = 2.0 * x;
  b1 = 0.0;
  b0 = 0.0;

  for ( i = n - 1; 0 <= i; i-- )
  {
    b2 = b1;
    b1 = b0;
    b0 = twox * b1 - b2 + a[i];
  }

  value = 0.5 * ( b0 - b2 );

  return value;
}

float r4_asin ( float x )
{
	static float asincs[20] = {
		0.10246391753227159E+00,
		0.054946487221245833E+00,
		0.004080630392544969E+00,
		0.000407890068546044E+00,
		0.000046985367432203E+00,
		0.000005880975813970E+00,
		0.000000777323124627E+00,
		0.000000106774233400E+00,
		0.000000015092399536E+00,
		0.000000002180972408E+00,
		0.000000000320759842E+00,
		0.000000000047855369E+00,
		0.000000000007225128E+00,
		0.000000000001101833E+00,
		0.000000000000169476E+00,
		0.000000000000026261E+00,
		0.000000000000004095E+00,
		0.000000000000000642E+00,
		0.000000000000000101E+00,
		0.000000000000000016E+00 };
	static int nterms = 0;
	const float pi2 = 1.57079632679489661923E+00;
	static float sqeps = 0.0;
	float value;
	float y;
	float z;
	
	if ( nterms == 0 )
	{
		nterms = r4_inits ( asincs, 20, 0.1 * r4_mach ( 3 ) );
		sqeps = r4_sqrt ( 6.0 * r4_mach ( 3 ) );
	}

	y = r4_abs ( x );

	if ( x < - 1.0 - sqeps )
	{
		fprintf ( stderr, "\n" );
		fprintf ( stderr, "R4_ASIN - Fatal error!\n" );
		fprintf ( stderr, "	X < - 1.0\n" );
		exit ( 1 );
	}
	else if ( x < - 1.0 )
	{
		value = - pi2;
	}
	else if ( x < 1.0 )
	{
		z = 0.0;
		if ( sqeps < y )
		{
			z = y * y;
		}

		if ( z <= 0.5 )
		{
			value = x * ( 1.0 + r4_csevl ( 4.0 * z - 1.0, asincs, nterms ) );
		}
		else
		{
			value = pi2 - r4_sqrt ( 1.0 - z ) * ( 1.0 + 
				r4_csevl ( 3.0 - 4.0 * z, asincs, nterms ) );
		}

		if ( x < 0.0 )
		{
			value = - r4_abs ( value );
		}
		else if ( 0.0 < x )
		{
			value = + r4_abs ( value );
		}
	}
	else if ( x < 1.0 + sqeps )
	{
		value = pi2;
	}
	else
	{
		fprintf ( stderr, "\n" );
		fprintf ( stderr, "R4_ASIN - Fatal error!\n" );
		fprintf ( stderr, "	1.0 < X\n" );
		exit ( 1 );
	}
	return value;
}

Je ne veux pas encore rentrer dans des polémiques, mais, le web est truffé de cette programmation la, et beaucoup de gens sont fiers de ce genre de chose. Pour ma part je trouve cela juste affligeant cette programmation...

...Bref, je vais me débrouiller avec ça, ou bien tu as une méthode avec "un tableau de constantes précalculées, un petit bout de code pour la gestion qui est une boucle et un test" ??

[Antoane]

Bonjour,

Est-ce que cette discussion ne serait pas davantage à sa place en informatique ? Oui si il est purement question d'algo, non si la discussion porte davantage sur comment mettre cet algo dans un µC http://forums.futura-sciences.com/pr...algorithmique/

PS : switch i / case, c'est plus propre que if /elseif/elseif/elseif... Non ?

[sylvainmahe]

Bonjour Antoane,

Oui tu as raison, à l'origine j'avais posé une question similaire dans le forum de mathématiques, mais les mathématiciens on en général du mal à me répondre parce que pour eux donner un exemple avec un calcul qui créé une approximation (et non pas quelque chose d'éventuellement juste sur l'infini ?? mais la on va loin..), c'est m'expliquer quelque chose de forcément faux, donc ils n'ont pas envie de m'expliquer quelque chose de faux, même si en programmation de systèmes embarqués l’approximation est incontournable... (c'est surtout impossible de faire autrement).

Alors sinon pour ta remarque du if else if else, ce code n'est pas de moi J'ai simplement essayé de retrouver toutes les fonctions dépendantes de la fonction arcsinus.


Bon j'ai un petit peu avancé Jiherve, j'avais pas vu que dans ton lien ce fichier était bien plus clair:
https://people.sc.fsu.edu/~jburkardt...dic/cordic.cpp

Du coup voila la fonction arccosinus qu'on y trouve:

double arccos_cordic ( double t, int n )

//****************************************************************************80
//
//  Purpose:
//
//    ARCCOS_CORDIC returns the arccosine of an angle using the CORDIC method.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Jean-Michel Muller,
//    Elementary Functions: Algorithms and Implementation,
//    Second Edition,
//    Birkhaeuser, 2006,
//    ISBN13: 978-0-8176-4372-0,
//    LC: QA331.M866.
//
//  Parameters:
//
//    Input, double T, the cosine of an angle.  -1 <= T <= 1.
//
//    Input, int N, the number of iterations to take.
//    A value of 10 is low.  Good accuracy is achieved with 20 or more
//    iterations.
//
//    Output, double ARCCOS_CORDIC, an angle whose cosine is T.
//
//  Local Parameters:
//
//    Local, double ANGLES(60) = arctan ( (1/2)^(0:59) );
//
{

  double angle;
  double angles[60] = {
    0.78539816339744830962, 
    0.46364760900080611621, 
    0.24497866312686415417, 
    0.12435499454676143503, 
    0.062418809995957348474, 
    0.031239833430268276254, 
    0.015623728620476830803, 
    0.0078123410601011112965, 
    0.0039062301319669718276, 
    0.0019531225164788186851, 
    0.00097656218955931943040, 
    0.00048828121119489827547, 
    0.00024414062014936176402, 
    0.00012207031189367020424, 
    0.000061035156174208775022, 
    0.000030517578115526096862, 
    0.000015258789061315762107, 
    0.0000076293945311019702634, 
    0.0000038146972656064962829, 
    0.0000019073486328101870354, 
    0.00000095367431640596087942, 
    0.00000047683715820308885993, 
    0.00000023841857910155798249, 
    0.00000011920928955078068531, 
    0.000000059604644775390554414, 
    0.000000029802322387695303677, 
    0.000000014901161193847655147, 
    0.0000000074505805969238279871, 
    0.0000000037252902984619140453, 
    0.0000000018626451492309570291, 
    0.00000000093132257461547851536, 
    0.00000000046566128730773925778, 
    0.00000000023283064365386962890, 
    0.00000000011641532182693481445, 
    0.000000000058207660913467407226, 
    0.000000000029103830456733703613, 
    0.000000000014551915228366851807, 
    0.0000000000072759576141834259033, 
    0.0000000000036379788070917129517, 
    0.0000000000018189894035458564758, 
    0.00000000000090949470177292823792, 
    0.00000000000045474735088646411896, 
    0.00000000000022737367544323205948, 
    0.00000000000011368683772161602974, 
    0.000000000000056843418860808014870, 
    0.000000000000028421709430404007435, 
    0.000000000000014210854715202003717, 
    0.0000000000000071054273576010018587, 
    0.0000000000000035527136788005009294, 
    0.0000000000000017763568394002504647, 
    0.00000000000000088817841970012523234, 
    0.00000000000000044408920985006261617, 
    0.00000000000000022204460492503130808, 
    0.00000000000000011102230246251565404, 
    0.000000000000000055511151231257827021, 
    0.000000000000000027755575615628913511, 
    0.000000000000000003877787807814456755, 
    0.0000000000000000069388939039072283776, 
    0.0000000000000000034694469519536141888, 
    0.0000000000000000017347234759768070944};
  int i;
  int j;
  double poweroftwo;
  double sigma;
  double sign_z2;
  double theta;
  double x1;
  double x2;
  double y1;
  double y2;


  theta = 0.0;
  x1 = 1.0;
  y1 = 0.0;
  poweroftwo = 1.0;

  for ( j = 1; j <= n; j++ )
  {
    if ( y1 < 0.0 )
    {
      sign_z2 = -1.0;
    }
    else
    {
      sign_z2 = +1.0;
    }

    if ( t <= x1 )
    {
      sigma = + sign_z2;
    }
    else
    {
      sigma = - sign_z2;
    }

    if ( j <= 60 )
    {
      angle = angles[j-1];
    }
    else
    {
      angle = angle * 0.5;
    }

    for ( i = 1; i <= 2; i++ )
    {
      x2 =                      x1 - sigma * poweroftwo * y1;
      y2 = sigma * poweroftwo * x1 +                      y1;

      x1 = x2;
      y1 = y2;
    }

    theta  = theta + 2.0 * sigma * angle;

    t = t + t * poweroftwo * poweroftwo;

    poweroftwo = poweroftwo * 0.5;
  }

  return theta;
}

Il est dit d'utiliser une itération au moins supérieur à 10, 20 pour être bien, j'ai fait quelques essais, par exemple pour arccosinus de 0.5 j'obtiens bien 60 degrés à partir de 12 itérations, par contre c'est niveau temps de calcul que ça ne va pas:

arccosinus de 0.5 = 1893 micro secondes

Alors que ma fonction qui fait des tests de valeurs de cosinus fait la même chose en 1675 micro seconde (ce qui déjà est très long).

?

[CM63]

Bonjour,

Pour améliorer la précision pour le calcul des tangentes des angles supérieurs à 45°, il suffit de calculer (pour )

[Paraboloide_Hyperbolique]

Je ne veux pas encore rentrer dans des polémiques, mais, le web est truffé de cette programmation la, et beaucoup de gens sont fiers de ce genre de chose. Pour ma part je trouve cela juste affligeant cette programmation...

Bonjour,

C'est malheureusement le prix a payer si l'on veut quelque chose qui soit à la fois performant et robuste. Le code que vous montrez tient compte des problèmes d'arrondis numériques, d'overflow et d'underflow; ce qui n'est pas le cas quand on code simplement une expansion de Taylor.

[jiherve]

Bonjour
quel processeur utilises tu ?
as tu besoin de calculer en 64 bits ?
Le nombre d’itérations est proportionnel à la taille des variables donc en float(32bits) il y en aura deux fois moins et si ton processeur n'est pas un 64 bits cela ira mieux aussi.
JR

[sylvainmahe]

Oui je comprend bien ce que tu dis Paraboloide_Hyperbolique mais c'est un peu fou tout ce code juste pour asin...


jiherve euh comment dire, j'utilise l'atmega328p ...


CM63 j'avais pensé à ce truc la, mais j'avais demandé à mon ancien prof de math, sans doute il ne m'avais pas compris, parce qu'il m'avais dit que ce n'était pas possible. (enfin on parlait des fonctions inverse je crois ce qui peut être n'est pas la même chose...)

Du coup j'ai testé ton truc et ca fonctionne vraiment nickel ! La précision est la et la vitesse aussi:
tan:

const float _DEGREE_TO_RADIAN = 0.0174532925199432;

//si angle inférieur ou égale à 45°:
alpha = ANGLE * _DEGREE_TO_RADIAN;
return alpha + ((0.3333333333333333 * (alpha * alpha * alpha)) + (0.1333333333333333 * (alpha * alpha * alpha * alpha * alpha)) + (0.0539682539682539 * (alpha * alpha * alpha * alpha * alpha * alpha * alpha)));

//si angle supérieur à 45°:
alpha = (90.0 - ANGLE) * _DEGREE_TO_RADIAN;
return 1.0 / (alpha + ((0.3333333333333333 * (alpha * alpha * alpha)) + (0.1333333333333333 * (alpha * alpha * alpha * alpha * alpha)) + (0.0539682539682539 * (alpha * alpha * alpha * alpha * alpha * alpha * alpha))));

Du coup d'après toi, c'est possible aussi pour arcsin supérieur à 0.5 (ou à 0.707 -> 45°)? :

const float _RADIAN_TO_DEGREE = 57.2957795130823208;

return (SINUS + ((0.5 * ((SINUS * SINUS * SINUS) * 0.3333333333333333)) + (0.375 * ((SINUS * SINUS * SINUS * SINUS * SINUS) * 0.2)) + (0.3125 * ((SINUS * SINUS * SINUS * SINUS * SINUS * SINUS * SINUS) * 0.1428571428571428)))) * _RADIAN_TO_DEGREE;

Merci pour votre aide

[invite82078308]

Bonjour,

Est-ce que cette discussion ne serait pas davantage à sa place en informatique ? Oui si il est purement question d'algo, non si la discussion porte davantage sur comment mettre cet algo dans un µC http://forums.futura-sciences.com/pr...algorithmique/
(...)

Le mathématicien peut avoir son mot à dire, après tout les fonctions trigonométriques, c'est lui qui les a inventées.
L'instinct mathématique le pousse à résoudre un problème par lui même plutôt qu'à aller chercher ailleurs.
La question ici est d'avoir à la fois une bonne efficacité et une bonne précision, et pour les solutions présentées, sous forme de programme, le résultat dépend de l'implémentation du langage.
Voyons plutôt ce qu'on peut utiliser comme algorithme.
Il me semble intéressant d'utiliser les propriétés de l'exponentielle complexe.
On a: e^{i x}= cos(x) + i sin(x), qu'on peut approcher par la formule:
e^z= lim_{n-> +∞} (1+ z/n)ⁿ
Comme cela, c'est moyen mais on peut aller bien plus vite en prenant n de la forme 2^m, diviser par une puissance de deux est rudimentaire et élever à la puissance 2^m se fait en élevant le nombre m fois au carré.
Évaluer l'erreur commise en fonction de x et m ...
Il y la question des arrondis.
n étant choisi, me semble intéressant pour cela de calculer, plutôt que u_h= (1+ z/n)^h, de calculer la suite v_h = u_h -1 en réels flottants par récurrence, voire utiliser des méthodes plus méthodes plus puissantes pour éviter de perdre trop de chiffres après la virgule.

Tout cela dépend en fait des spécifications, en particulier sur la précision.

Pour les fonctions trigonométrique inverses, je ne sais pas trop ...
Utiliser la méthode de newton pour résoudre numériquent par exemple en t: a= cos t ???
en prenant la valeur donné par la fonction arc cos standard en première approximation ...

On notera que les relations trigonométriques habituelles permettent de tester la précision des calculs (prévoir aussi une bonne approximation de )

[CM63]

Bonjour,

Du coup d'après toi, c'est possible aussi pour arcsin supérieur à 0.5 (ou à 0.707 -> 45°)? :

Non mais pour l'arctg oui : si le nombre est supérieur à 1, on prend l'inverse, on calcule l'arctg en utilisant la série, et ensuite on prend

[jiherve]

Bonjour
sur une machine utilisant un ATMEGA 328 il faut utiliser des calculs en entier ou en complément à deux, le tout en assembleur, si l'on veut que cela "pousse".
JR

[sylvainmahe]

Schrodies-cat, c'est pas que c'est du chinois ce que tu écris, c'est qu'il me manque des choses dans ma culture pour tout comprendre

Par exemple la lettre m, je ne sais pas ce que c'est, donc je vais voir sur wikipedia et on me dit:
" Pour tout entier naturel n, il existe un autre entier naturel m tel que m soit supérieur ou égal à n.
Tout entier naturel est inférieur ou égal à au moins un autre entier naturel. "

Bon faut que je cogite...


CM63, ok. Donc j'ai fait un essai pour arctan, bizarrement le "seuil" ou la précision semble se dégrader de façon égale pour les 2 façons de calculer ce situe à une tangente égale à environ 0.66:

float Core::arctan (const float TANGENT)
{
	float tangent2 = 0;
	float tangent3 = 0;
	float tangent5 = 0;
	const float BND = 0.66;
	float angle = 0;
	
	tangent2 = TANGENT * TANGENT;
	tangent3 = tangent2 * TANGENT;
	tangent5 = tangent2 * tangent3;
	
	if (TANGENT < BND)
	{
		angle = (TANGENT - ((tangent3 * 0.3333333333333333) + (tangent5 * 0.2) - ((tangent2 * tangent5) * 0.1428571428571428) + ((tangent2 * tangent2 * tangent5) * 0.11111111111111111111111111111111) - ((TANGENT * tangent5 * tangent5) * 0.09090909090909090909091) + ((tangent3 * tangent5 * tangent5) * 0.07692307692307692307692) - ((tangent5 * tangent5 * tangent5) * 0.06666666666666666666667))) * _RADIAN_TO_DEGREE;
	}
	else
	{
		angle = (1.5707963267948966 - (1.0 / (TANGENT - ((tangent3 * 0.3333333333333333) + (tangent5 * 0.2) - ((tangent2 * tangent5) * 0.1428571428571428)) + ((tangent2 * tangent2 * tangent5) * 0.11111111111111111111111111111111) - ((TANGENT * tangent5 * tangent5) * 0.09090909090909090909091) + ((tangent3 * tangent5 * tangent5) * 0.07692307692307692307692) - ((tangent5 * tangent5 * tangent5) * 0.06666666666666666666667)))) * _RADIAN_TO_DEGREE + 45;
	}
	
	return angle;
}

Malgré ça dans certaines valeurs l'erreur max est d'environ 2 degrés. Mais bon je ne vais pas vous embêter plus longtemps avec mes histoires de reproduire quelques fonctions de math.h, vous m'avez déjà bien aidé, merci


jiherve je sais bien qu'on devrait tout coder en assembleur (blague ou pas), mais pour prendre un exemple, mes fonctions sinus, cosinus et tangente s’exécutent de 130 à 180 micro seconde sur l'atmega328p, c'est pas si pourri que ça !

Si tu veux recoder 5100 lignes de c en assembleur (d'un projet open source), ta participation sera appréciée

Bon après effectivement rien ne m’empêche de coder CERTAINES chose de ma lib en assembleur, quand comme tu le dis on a vraiment besoin de plus de vitesse... Pour l'instant le problème ne s'est posé que pour les fonctions de math évoquées ici (en plus d'un sqrt que j'ai réussi à descendre dans les 200 micro seconde en regardant quelques exemples sur internet).

[inviteb9f49292]

jiherve je sais bien qu'on devrait tout coder en assembleur (blague ou pas), mais pour prendre un exemple, mes fonctions sinus, cosinus et tangente s’exécutent de 130 à 180 micro seconde sur l'atmega328p, c'est pas si pourri que ça !

J'irai pas jusqu'à l'assembleur (sauf comme tu dis éventuellement pour des "coeurs d'algo") par contre je plussoie l'utilisation de représentation virgule fixe plutôt que flottante. Si tu veux être un peu "pro" dans le domaine du micro-contrôleur c'est une notion fondamentale.

[jiherve]

Bonjour,

Si tu veux recoder 5100 lignes de c en assembleur (d'un projet open source), ta participation sera appréciée

non merci j'ai du coder en langages divers et variés durant la quarantaine d'année de ma carrière quelques centaines de milliers de lignes j'ai donc bien donné.
Ceci dit pour les besoin d'un calculateur destiné à la balistique j'ai redéveloppé toute une librairie mathématique pour ATMEGA avec un format flottant "propriétaire" inspiré de l'IEEE 754.
JR

[sylvainmahe]

Merci pour votre aide en tout cas vous m'avez bien fait comprendre certaines choses

[sylvainmahe]

Pour lou_ibmix_xi, avant je faisais encore autrement que la virgule fixe (mais peut être ça s'appelle la virgule fixe), je multipliais par exemple par 1000 les calculs concernés, et je faisais donc tout avec des entiers.

[CM63]

Bonjour,

Je te propose une autre méthode:
- ramener les calculs entre 0 et pi/2 (ou pi/4 pour la tg) grâce aux formules d'équivalence,
- pour le cos utiliser un produit infini (recherche dans Wikipedia) plutôt que le développement en polynôme : la convergence est beaucoup plus rapide,
- pour le sin utiliser
- pour la tg : sin/cos
etc.

[jiherve]

Re
Non les methodes usuelles de calcul utilisent soit des approximations polynomiales soit le CORDIC ,cas de nombreuse calculettes car il y a de nombreux pièges dans les calculs utilisant des formats à la résolution limités ce qui est le cas des "flottants", ces deux méthodes n'ont pas vraiment de soucis de convergence(au sens durée) car les algorithmes sont déterministes, pour la précision (donc convergence vers la valeur "vraie") c'est autre chose les CORDIC sont très sensibles aux valeurs de leurs tables de constantes il faut faire très attention aux arrondis effectués lors de leurs définition,idem pour les coeff des polynomes.
JR

[inviteb9f49292]

Pour lou_ibmix_xi, avant je faisais encore autrement que la virgule fixe (mais peut être ça s'appelle la virgule fixe), je multipliais par exemple par 1000 les calculs concernés, et je faisais donc tout avec des entiers.

Je ne pense pas que ça revienne au même. La virgule fixe utilise les entiers de la CPU, c'est toi qui l'interprète comme un nombre à virgule, si tu es en 32 bits, tu peux par exemple décider de prendre 15 bits après la virgule (le plus petit représentable est 1/2^15), tu auras donc 16 bits entiers (le plus gros sera 2^16). C'est vraiment une technique très utile

[Ludwig1]

Je ne pense pas que ça revienne au même. La virgule fixe utilise les entiers de la CPU, c'est toi qui l'interprète comme un nombre à virgule, si tu es en 32 bits, tu peux par exemple décider de prendre 15 bits après la virgule (le plus petit représentable est 1/2^15), tu auras donc 16 bits entiers (le plus gros sera 2^16). C'est vraiment une technique très utile

Salut,

Comme je me suis fait traiter de vieux C.. il n'y a pas longtemps, je propose tout de même la méthode que les vieux C.. utilisaient il y a fort longtemps.
On fait une table 0 à pi/2 par exemple les valeurs de sin avec 4 chiffres après la virgule, ( nombres entiers, on incrémente la table ° par ° )
pour trouver le sin on prend les valeurs min max dans la table et on fait un calcul d'interpolation, c'est très rapide.
et c'est très précis.
Dans la pratique, 3 chiffres après la virgule sont suffisants.
Pour les calculs, à la fin, ne pas oublier de diviser par 1000 ou 10000 selon la position de la virgule.
On se fabriquera dans le langage avec lequel on travaille, selon le processeur cible, les routines Necessaires.

Cordialement

Ludwig

[sylvainmahe]

Ludwig1 c'est justement ce que j'étais en train d'écrire dans ma programmation

J'ai fait pas mal de tests depuis sur les autres possibles, ça n'est pas assez rapide, donc je pense que ceci est une solution, j'ai pas encore testé la rapidité de ce truc la mais ça ne saurait tarder...

Par contre moi j'écris ça en vieux C++.., pas en vieux C.. (blague )

[inviteb9f49292]

Ludwig1 c'est justement ce que j'étais en train d'écrire dans ma programmation

J'ai fait pas mal de tests depuis sur les autres possibles, ça n'est pas assez rapide, donc je pense que ceci est une solution, j'ai pas encore testé la rapidité de ce truc la mais ça ne saurait tarder...

Utiliser une "look-up table" ou cordic ou des polynomes c'est orthogonal avec la représentation que tu utilises pour tes nombres... Et la précision nécessaire dépends franchement de l'application visée, donc non, dans la pratique 3 chiffres après la virgule ne sont pas toujours suffisant... Surtout que le problème se situe essentiellement dans la diffusion de l'erreur, les erreurs d'arrondis des calculs intermédiaires font boule de neige et tu te retrouves avec un résultat complètement foireux... Bref, selon le type de calculs, pour avoir 3 chiffres significatifs au résultat, il faut que tu en gardes 5 ou 6 durant tes calculs intermédiaires... Pour ce qui est de la représentation, on peut faire ce qu'on veut mais ce n'est pas une raison... Soit tu as un coeur flottant, et a part exception tu travailles en flottant, sinon je te conseille la virgule fixe qui est la "solution d'école", il se peut que pour certains cas particulier il y ait des représentations plus adaptés à la problématique, mais amha c'est un choix d'expert, et un expert doit connaître ses outils de bases pour savoir s'ils sont adaptés ou pas...

Par contre moi j'écris ça en vieux C++.., pas en vieux C.. (blague )

Oui je sens la "private joke" suite à une autre discussion... Mais j'en profite pour rebondir là-dessus et faire mon vieux c** (non pas un langage, mais bien l'adjectif): en embarqué nous n'aimons pas beaucoup le C++, je n'entrerai pas dans le débat ici c'est hors-sujet, mais c'est du C (le langage cette fois ci), parfois de l'ADA, et il y a pas mal de raisons objectives à cet état de fait.

[sylvainmahe]

Ce que tu dis ne changera rien, et je pense que tu le sais
Il y a 2 ans si j'avais écouté ce genre de chose je ne serais arrivé à rien, c'est d’ailleurs ce qu'ILS (les gens de la haute sphère de l'électronique embarqué) prédisaient, mais ils s'étaient trompés. Mais ça je n'y peu rien.

[inviteb9f49292]

Ce que tu dis ne changera rien, et je pense que tu le sais
Il y a 2 ans si j'avais écouté ce genre de chose je ne serais arrivé à rien, c'est d’ailleurs ce qu'ILS (les gens de la haute sphère de l'électronique embarqué) prédisaient, mais ils s'étaient trompés. Mais ça je n'y peu rien.

C'est hors sujet je n'insisterai pas (surtout que je n'ai rien à y gagner), mais "les gens de la haute sphère de l'électronique embarqué" c'est leur gagne pain, ce n'est pas parce que ton exigence de qualité/fiabilité n'est pas au même niveau que le leurs qu'ils racontent des conneries. Après c'est sûr, pour bricoler un bidule où la fiabilité la maintenance etc... ne sont pas des contraintes majeures on fait ce que l'on veut.
Et donc tu crois vraiment que ce que tu as écris en C++ ne pouvait pas l'être en C ? Et tu penses que les gens qui expliquent que le C++ n'est pas adapté à l'embarqué n'ont en tête que l'aspect "rédactionnel" du développement logiciel ? C'est que à priori tu n'as pas compris quels sont les problèmes que posent le C++ dans le monde de l'embarqué. Mais je le répète, c'est toi qui y gagnerai à écouter les gens du métiers...

[sylvainmahe]

Moi j'aime bien partir d'un projet en manipulant des objets, je trouve ça très pratique, voila tout. J'ai à peine été regardé les inconvénients ou les avantages de cette programmation objet sur des microcontrôleurs, puisque je n'ai jamais eu de bug ou de problème à cause de ça. Quand j'en aurais un j'irais voir.

J'impose rien à personne, chacun est libre de programmer comme il veut, et encore heureusement.

[Ludwig1]

C'est hors sujet je n'insisterai pas (surtout que je n'ai rien à y gagner), mais "les gens de la haute sphère de l'électronique embarqué" c'est leur gagne pain, ce n'est pas parce que ton exigence de qualité/fiabilité n'est pas au même niveau que le leurs qu'ils racontent des conneries. Après c'est sûr, pour bricoler un bidule où la fiabilité la maintenance etc... ne sont pas des contraintes majeures on fait ce que l'on veut.
Et donc tu crois vraiment que ce que tu as écris en C++ ne pouvait pas l'être en C ? Et tu penses que les gens qui expliquent que le C++ n'est pas adapté à l'embarqué n'ont en tête que l'aspect "rédactionnel" du développement logiciel ? C'est que à priori tu n'as pas compris quels sont les problèmes que posent le C++ dans le monde de l'embarqué. Mais je le répète, c'est toi qui y gagnerai à écouter les gens du métiers...

Salut,
je dirai que chacun voit midi à sa porte, simplement je dirai aussi qu'il y a bien des façons de résoudre ce type de PB. Dans le monde de l'embarqué on utilise bien des façons de faire. A.m.h.a, il n'y a pas de méthode standard, tous dépend du PB à traiter.
Connaissant l'histoire, comment Ritchie et ses deux compères on pondus C je ne suis ni convaincu par C ni C++ mais chacun fait comme il pense.
Mes PB se situant dans le domaine de l'automatique, j'utilise des outils mis au point pour ce genre d'applications.
Cordialement
Ludwig

[Ludwig1]

Moi j'aime bien partir d'un projet en manipulant des objets, je trouve ça très pratique, voila tout. J'ai à peine été regardé les inconvénients ou les avantages de cette programmation objet sur des microcontrôleurs, puisque je n'ai jamais eu de bug ou de problème à cause de ça. Quand j'en aurais un j'irais voir.

J'impose rien à personne, chacun est libre de programmer comme il veut, et encore heureusement.

Salut,
Pour ce qui est des objets je partage ton avis, je fais de même, simplement pas avec les méthodes usuelles.

Cordialement

Ludwig