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
Index
Trunk
Branches
Tags
Trunk
Branches
Tags
Home page
Home page
libs/libmpu.git/trunk/mpu/st-cos.c
log
diff
blame

/***************************************************************
  __ST_COS.C

       This file contains source code of functions for
       COS 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_cos_001_emu32lsb.dfn -   128-бит,
      ei_cos_512_emu32lsb.dfn - 65536-бит.

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

#if MPU_MATH_FN_LIMIT >= 128
#if MPU_WORD_ORDER_BIG_ENDIAN == 0
#include <math/cos/emu00128/ei_cos_001_emu32lsb.dfn>
#else
#include <math/cos/emu00128/ei_cos_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/cos/emu00256/ei_cos_002_emu32lsb.dfn>
#else
#include <math/cos/emu00256/ei_cos_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/cos/emu00512/ei_cos_004_emu32lsb.dfn>
#else
#include <math/cos/emu00512/ei_cos_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/cos/emu01024/ei_cos_008_emu32lsb.dfn>
#else
#include <math/cos/emu01024/ei_cos_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/cos/emu02048/ei_cos_016_emu32lsb.dfn>
#else
#include <math/cos/emu02048/ei_cos_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/cos/emu04096/ei_cos_032_emu32lsb.dfn>
#else
#include <math/cos/emu04096/ei_cos_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/cos/emu08192/ei_cos_064_emu32lsb.dfn>
#else
#include <math/cos/emu08192/ei_cos_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/cos/emu16384/ei_cos_128_emu32lsb.dfn>
#else
#include <math/cos/emu16384/ei_cos_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/cos/emu32768/ei_cos_256_emu32lsb.dfn>
#else
#include <math/cos/emu32768/ei_cos_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/cos/emu65536/ei_cos_512_emu32lsb.dfn>
#else
#include <math/cos/emu65536/ei_cos_512_emu32msb.dfn>
#endif
#endif /* MPU_MATH_FN_LIMIT >= 65536 */


static int _get_n_cos( 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_COS_C128;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 128 */
#if MPU_MATH_FN_LIMIT >= 256
    case NBR_256  :
      rc = N_COS_C256;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 256 */
#if MPU_MATH_FN_LIMIT >= 512
    case NBR_512  :
      rc = N_COS_C512;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 512 */
#if MPU_MATH_FN_LIMIT >= 1024
    case NBR_1024 :
      rc = N_COS_C1024;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 1024 */
#if MPU_MATH_FN_LIMIT >= 2048
    case NBR_2048 :
      rc = N_COS_C2048;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 2048 */
#if MPU_MATH_FN_LIMIT >= 4096
    case NBR_4096 :
      rc = N_COS_C4096;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 4096 */
#if MPU_MATH_FN_LIMIT >= 8192
    case NBR_8192 :
      rc = N_COS_C8192;
      break;
#endif /* MPU_MATH_FN_LIMIT >= 8192 */
#if MPU_MATH_FN_LIMIT >= 16384
    case NBR_16384:
      rc = N_COS_C16384;
      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_cos() */


static EMUSHORT *_get_cos__C_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_cos__C_128_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 128 */
#if MPU_MATH_FN_LIMIT >= 256
    case NBR_256  :
      rc = (EMUSHORT *)&_ei_cos__C_256_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 256 */
#if MPU_MATH_FN_LIMIT >= 512
    case NBR_512  :
      rc = (EMUSHORT *)&_ei_cos__C_512_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 512 */
#if MPU_MATH_FN_LIMIT >= 1024
    case NBR_1024 :
      rc = (EMUSHORT *)&_ei_cos__C_1024_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 1024 */
#if MPU_MATH_FN_LIMIT >= 2048
    case NBR_2048 :
      rc = (EMUSHORT *)&_ei_cos__C_2048_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 2048 */
#if MPU_MATH_FN_LIMIT >= 4096
    case NBR_4096 :
      rc = (EMUSHORT *)&_ei_cos__C_4096_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 4096 */
#if MPU_MATH_FN_LIMIT >= 8192
    case NBR_8192 :
      rc = (EMUSHORT *)&_ei_cos__C_8192_[0][0];
      break;
#endif /* MPU_MATH_FN_LIMIT >= 8192 */
#if MPU_MATH_FN_LIMIT >= 16384
    case NBR_16384:
      rc = (EMUSHORT *)&_ei_cos__C_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_cos__C_ptr() */


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

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

                                                          EIX*EIX
 Concepts           : RETURN in EIY = cos( EIX*k ) - 1 + --------- 
                                                             2
                      on [-PI/4, PI/4], where k = pi/PI,
                      PI is the rounted value of pi in machine
                      precision.

 METHOD             : Let z = x*x. Create a polinomial
                      approximation to
                         cos(k*x)-1+z/2 =
                              z*z*(C0 + C1*z^1 + C2*z^2 + ... ).
                      Then
                         ei_cos__C(x*x) =
                              z*z*(C0 + C1*z^1 + C2*z^2 + ... ).

 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_cos( 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_cos__C_ptr( nb );
   for( i = 0; i < n; i++ )
   {
     ei_add( y, y, p, nb );
     ei_mul( y, y, x, nb );
     p += np;
   }

   ei_mul( y, y, x, nb );

   ei_copy( eiy, y, nb );

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

} /* End of ei_cos__C() */


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

__mpu_hidden_decl(ei_cos__C);


/*
  End of hide internal symbols.
 ***************************************************************/
log
diff
blame
cGit-ui™ 0.1.7
Git™ 2.40.0
Nginx™ 1.22.1

Where any material of this site is being reproduced, published or issued to others the reference to the source is obligatory.

© Andrey V. Kosteltsev, 2019 – 2025.