Math Processor Unit Library

libmpu – library of arithmetic functions for integer, real, and complex numbers of increased digit capacity

16 Commits   0 Branches   2 Tags

/***************************************************************
  __ST_ATAN2.C

       This file contains source code of functions for
       ATAN2 constants operations.

       PART OF : MPU - library .

       USAGE   : Internal only .

       NOTE    : NONE .

       Copyright (C) 2000 - 2024  by Andrew V.Kosteltsev.
       All Rights Reserved.
 ***************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <errno.h>   /* errno(3)  */
#include <string.h>  /* strcpy(3) */
#include <strings.h> /* bzero(3)  */
#include <stdlib.h>

#include <libmpu.h>
#include <mpu-context.h>

#include <mpu-emutype.h>
#include <mpu-integer.h>
#include <mpu-real.h>
#include <mpu-floatp.h>

#include <mpu-char.h>
#include <mpu-symbols.h>

#include <mpu-math-errno.h>
#include <mpu-mtherr.h>


/***************************************************************
  Кодировка имен файлов:

  Трехзначное десятичное число, представляющее количество
  128-и битных слов, из которых состоят вещественные числа
  размещенные в массивах:

    размер чисел в битах  кодировка
    --------------------  ---------
                     128  001
                     256  002
                     512  004
                    1024  008
                    2048  016
                    4096  032
                    8192  064
                   16384  128
                   32768  256
                   65536  512 (это предел);

    ПРИМЕРЫ:
    -------
      ei_atan2_001_emu32lsb.dfn -   128-бит,
      ei_atan2_512_emu32lsb.dfn - 65536-бит.

 ***************************************************************/

#if MPU_MATH_FN_LIMIT >= 128
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu00128/ei_atan2_001_emu32lsb.dfn>
#else
#include <math/atan2/emu00128/ei_atan2_001_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 128 */

#if MPU_MATH_FN_LIMIT >= 256
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu00256/ei_atan2_002_emu32lsb.dfn>
#else
#include <math/atan2/emu00256/ei_atan2_002_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 256 */

#if MPU_MATH_FN_LIMIT >= 512
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu00512/ei_atan2_004_emu32lsb.dfn>
#else
#include <math/atan2/emu00512/ei_atan2_004_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 512 */

#if MPU_MATH_FN_LIMIT >= 1024
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu01024/ei_atan2_008_emu32lsb.dfn>
#else
#include <math/atan2/emu01024/ei_atan2_008_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 1024 */

#if MPU_MATH_FN_LIMIT >= 2048
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu02048/ei_atan2_016_emu32lsb.dfn>
#else
#include <math/atan2/emu02048/ei_atan2_016_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 2048 */

#if MPU_MATH_FN_LIMIT >= 4096
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu04096/ei_atan2_032_emu32lsb.dfn>
#else
#include <math/atan2/emu04096/ei_atan2_032_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 4096 */

#if MPU_MATH_FN_LIMIT >= 8192
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu08192/ei_atan2_064_emu32lsb.dfn>
#else
#include <math/atan2/emu08192/ei_atan2_064_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 8192 */

#if MPU_MATH_FN_LIMIT >= 16384
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu16384/ei_atan2_128_emu32lsb.dfn>
#else
#include <math/atan2/emu16384/ei_atan2_128_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 16384 */

#if MPU_MATH_FN_LIMIT >= 32768
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu32768/ei_atan2_256_emu32lsb.dfn>
#else
#include <math/atan2/emu32768/ei_atan2_256_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 32768 */

#if MPU_MATH_FN_LIMIT >= 65536
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/atan2/emu65536/ei_atan2_512_emu32lsb.dfn>
#else
#include <math/atan2/emu65536/ei_atan2_512_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 65536 */


static int _get_n_atan2( int nb )
{
  int  rc = 0;

  if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT )
  {
    /* error: Invalid size of operand(s) */
    __real_error_no = __R_ESIZE__;
    __STIND; /* Set REAL ind-produsing operation Flag */
    return( rc );
  }

  switch( nb )
  {
#if MPU_MATH_FN_LIMIT >= 128
    case NBR_32   :
    case NBR_64   :
    case NBR_128  :
      rc = N_ATAN2_A128;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 128 */
#if MPU_MATH_FN_LIMIT >= 256
    case NBR_256  :
      rc = N_ATAN2_A256;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 256 */
#if MPU_MATH_FN_LIMIT >= 512
    case NBR_512  :
      rc = N_ATAN2_A512;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 512 */
#if MPU_MATH_FN_LIMIT >= 1024
    case NBR_1024 :
      rc = N_ATAN2_A1024;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 1024 */
#if MPU_MATH_FN_LIMIT >= 2048
    case NBR_2048 :
      rc = N_ATAN2_A2048;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 2048 */
#if MPU_MATH_FN_LIMIT >= 4096
    case NBR_4096 :
      rc = N_ATAN2_A4096;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 4096 */
#if MPU_MATH_FN_LIMIT >= 8192
    case NBR_8192 :
      rc = N_ATAN2_A8192;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 8192 */
#if MPU_MATH_FN_LIMIT >= 16384
    case NBR_16384:
      rc = N_ATAN2_A16384;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 16384 */

    default:
    {
      /* error: Invalid size of operand(s) */
      __real_error_no = __R_ESIZE__;
      __STIND; /* Set REAL ind-produsing operation Flag */
      break;
    }

  } /* End of switch( nb ) */

  return( rc );

} /* End of _get_n_atan2() */


static EMUSHORT *_get_atan2__A_ptr( int nb )
{
  EMUSHORT *rc = (EMUSHORT *)NULL;

  if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT )
  {
    /* error: Invalid size of operand(s) */
    __real_error_no = __R_ESIZE__;
    __STIND; /* Set REAL ind-produsing operation Flag */
    return( rc );
  }

  switch( nb )
  {
#if MPU_MATH_FN_LIMIT >= 128
    case NBR_32   :
    case NBR_64   :
    case NBR_128  :
      rc = (EMUSHORT *)&_ei_atan2__A_128_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 128 */
#if MPU_MATH_FN_LIMIT >= 256
    case NBR_256  :
      rc = (EMUSHORT *)&_ei_atan2__A_256_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 256 */
#if MPU_MATH_FN_LIMIT >= 512
    case NBR_512  :
      rc = (EMUSHORT *)&_ei_atan2__A_512_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 512 */
#if MPU_MATH_FN_LIMIT >= 1024
    case NBR_1024 :
      rc = (EMUSHORT *)&_ei_atan2__A_1024_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 1024 */
#if MPU_MATH_FN_LIMIT >= 2048
    case NBR_2048 :
      rc = (EMUSHORT *)&_ei_atan2__A_2048_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 2048 */
#if MPU_MATH_FN_LIMIT >= 4096
    case NBR_4096 :
      rc = (EMUSHORT *)&_ei_atan2__A_4096_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 4096 */
#if MPU_MATH_FN_LIMIT >= 8192
    case NBR_8192 :
      rc = (EMUSHORT *)&_ei_atan2__A_8192_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 8192 */
#if MPU_MATH_FN_LIMIT >= 16384
    case NBR_16384:
      rc = (EMUSHORT *)&_ei_atan2__A_16384_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 16384 */

    default:
    {
      /* error: Invalid size of operand(s) */
      __real_error_no = __R_ESIZE__;
      __STIND; /* Set REAL ind-produsing operation Flag */
      break;
    }

  } /* End of switch( nb ) */

  return( rc );

} /* End of _get_atan2__A_ptr() */


void ei_atan2__A( EMUSHORT *eiy, EMUSHORT *eix, int nb )
/***************************************************************

 Description        : ei_atan2__A() Работает с
                                    internal e-type data struct.

 Concepts           : SEE: ei_atan2() [polinom A].

 ACCURACY           : In the absence of rounding error,
                      the appriximation has absolute error
                      less than EPSILON [see: FLOATP.H].

 NOTE               : The coefficient C's are obtained by
                      a special Remez algorithm.

            ===============================================
            В INTERNET я нашел следующие алгоритмы:

            409  cacm  355  356 14  5  May 1971 e2 A60
            ------------------------------------------
               H. Schmitt;
                  Discrete {Chebychev} Curve Fit
                  approximation;Chebyshev approximation;
                  Chebyshev curve fitting;
                 +Chebyshev polynomial;
                  curve approximation;curve fitting;
                  exchange algorithm;
                 +polynomial approximation;Remez algorithm;

            501  toms   95   97  2  1  March 1976 e2 F K2
            ---------------------------------------------
               J. C. Simpson;
                  {FORTRAN} Translation of Algorithm 409
                  Discrete {Chebyshev} Curve Fit
                  approximation;polynomial approximation;
                  exchange algorithm;
                 +Chebyshev approximation;
                  polynomial approximation;
                  R,toms,95,4,1,March,1978,F. Futrell;

            последний из которых я перевел на "С", затем
            на язык операций повышенной разрядности.

            ===============================================

 Use Global Variable:

 Use Functions      : internal_np( nb );          | real.c
                      ei_copy( eiy, eix, nb );    | real.c
                      _gen_zero( eic, nb );       | real.c
                      ei_add( eic,eia,eib,nb );   | real.c
                      ei_mul( eic,eia,eib,nb );   | real.c

 Parameters         : EMUSHORT *eiy; - указатель на
                                       internal e-type
                                       data struct.
                                       TARGET;
                      EMUSHORT *eix; - указатель на
                                       internal e-type
                                       data struct.
                                       SOURCE;
                      int nb;        - количество бит в
                                       external e-type
                                       data struct.

 Return             : [void]

 ***************************************************************/
{
   EMUSHORT *x = NULL, *y = NULL;
   EMUSHORT *p;
   int       np;
   int       i, n = _get_n_atan2( nb );

   np = internal_np( nb );

   /*** Allocate memory for x, y . *****************************/
   x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) );
   if( !x )
   {
     /* fatal error */
     return;
   }

   y = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) );
   if( !y )
   {
     /* fatal error */

     /* FREE x *****************/
     __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) );
     /**************************/

     return;
   }
   /************************************************************/

   _gen_zero( y, nb );
   ei_copy( x, eix, nb );
   p = _get_atan2__A_ptr( nb );
   for( i = 0; i < n; i++ )
   {
     ei_add( y, y, p, nb );
     ei_mul( y, y, x, nb );
     p += np;
   }

   ei_copy( eiy, y, nb );

   /* FREE x *****************/
   /* FREE y *****************/
   __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) );
   /**************************/

} /* End of ei_atan2__A() */


/***************************************************************
  Hide internal symbols:
 ***************************************************************/

__mpu_hidden_decl(ei_atan2__A);


/*
  End of hide internal symbols.
 ***************************************************************/