Math Processor Unit Library

diff --git a/mpu/st-ln.c b/mpu/st-ln.c
new file mode 100644
index 0000000..525784d
--- /dev/null
+++ b/mpu/st-ln.c
@@ -0,0 +1,420 @@
+
+/***************************************************************
+  __ST_LN.C
+
+       This file contains source code of functions for
+       LN 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_ln_001_emu32lsb.dfn -   128-бит,
+      ei_ln_512_emu32lsb.dfn - 65536-бит.
+
+ ***************************************************************/
+
+#if MPU_MATH_FN_LIMIT >= 128
+#if MPU_WORD_ORDER_BIG_ENDIAN == 0
+#include <math/ln/emu00128/ei_ln_001_emu32lsb.dfn>
+#else
+#include <math/ln/emu00128/ei_ln_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/ln/emu00256/ei_ln_002_emu32lsb.dfn>
+#else
+#include <math/ln/emu00256/ei_ln_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/ln/emu00512/ei_ln_004_emu32lsb.dfn>
+#else
+#include <math/ln/emu00512/ei_ln_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/ln/emu01024/ei_ln_008_emu32lsb.dfn>
+#else
+#include <math/ln/emu01024/ei_ln_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/ln/emu02048/ei_ln_016_emu32lsb.dfn>
+#else
+#include <math/ln/emu02048/ei_ln_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/ln/emu04096/ei_ln_032_emu32lsb.dfn>
+#else
+#include <math/ln/emu04096/ei_ln_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/ln/emu08192/ei_ln_064_emu32lsb.dfn>
+#else
+#include <math/ln/emu08192/ei_ln_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/ln/emu16384/ei_ln_128_emu32lsb.dfn>
+#else
+#include <math/ln/emu16384/ei_ln_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/ln/emu32768/ei_ln_256_emu32lsb.dfn>
+#else
+#include <math/ln/emu32768/ei_ln_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/ln/emu65536/ei_ln_512_emu32lsb.dfn>
+#else
+#include <math/ln/emu65536/ei_ln_512_emu32msb.dfn>
+#endif
+#endif /* MPU_MATH_FN_LIMIT >= 65536 */
+
+
+static int _get_n_ln( 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_LOG_N128;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 128 */
+#if MPU_MATH_FN_LIMIT >= 256
+    case NBR_256  :
+      rc = N_LOG_N256;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 256 */
+#if MPU_MATH_FN_LIMIT >= 512
+    case NBR_512  :
+      rc = N_LOG_N512;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 512 */
+#if MPU_MATH_FN_LIMIT >= 1024
+    case NBR_1024 :
+      rc = N_LOG_N1024;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 1024 */
+#if MPU_MATH_FN_LIMIT >= 2048
+    case NBR_2048 :
+      rc = N_LOG_N2048;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 2048 */
+#if MPU_MATH_FN_LIMIT >= 4096
+    case NBR_4096 :
+      rc = N_LOG_N4096;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 4096 */
+#if MPU_MATH_FN_LIMIT >= 8192
+    case NBR_8192 :
+      rc = N_LOG_N8192;
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 8192 */
+#if MPU_MATH_FN_LIMIT >= 16384
+    case NBR_16384:
+      rc = N_LOG_N16384;
+      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_ln() */
+
+
+static EMUSHORT *_get_log__N_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_log__N_128_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 128 */
+#if MPU_MATH_FN_LIMIT >= 256
+    case NBR_256  :
+      rc = (EMUSHORT *)&_ei_log__N_256_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 256 */
+#if MPU_MATH_FN_LIMIT >= 512
+    case NBR_512  :
+      rc = (EMUSHORT *)&_ei_log__N_512_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 512 */
+#if MPU_MATH_FN_LIMIT >= 1024
+    case NBR_1024 :
+      rc = (EMUSHORT *)&_ei_log__N_1024_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 1024 */
+#if MPU_MATH_FN_LIMIT >= 2048
+    case NBR_2048 :
+      rc = (EMUSHORT *)&_ei_log__N_2048_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 2048 */
+#if MPU_MATH_FN_LIMIT >= 4096
+    case NBR_4096 :
+      rc = (EMUSHORT *)&_ei_log__N_4096_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 4096 */
+#if MPU_MATH_FN_LIMIT >= 8192
+    case NBR_8192 :
+      rc = (EMUSHORT *)&_ei_log__N_8192_[0][0];
+      break;
+#endif /* MPU_MATH_FN_LIMIT >= 8192 */
+#if MPU_MATH_FN_LIMIT >= 16384
+    case NBR_16384:
+      rc = (EMUSHORT *)&_ei_log__N_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_log__N_ptr() */
+
+
+void ei_log__N( EMUSHORT *eiy, EMUSHORT *eix, int nb )
+/***************************************************************
+
+ Description        : ei_log__N() Работает с
+                                  internal e-type data struct.
+
+ Concepts           : SEE: ei_log() [polinom A].
+
+ ACCURACY           : In the absence of rounding error,
+                      the appriximation has absolute error
+                      less than EPSILON [see: FLOATP.H].
+
+ NOTE               : The coefficient L'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_ln( 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_log__N_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_log__N() */
+
+
+/***************************************************************
+  Hide internal symbols:
+ ***************************************************************/
+
+__mpu_hidden_decl(ei_log__N);
+
+
+/*
+  End of hide internal symbols.
+ ***************************************************************/