/*************************************************************** __ST_LOG.C This file contains source code of functions for LOG 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 #endif #include /* errno(3) */ #include /* strcpy(3) */ #include /* bzero(3) */ #include #include #include #include #include #include #include #include #include #include #include /*************************************************************** Кодировка имен файлов: Трехзначное десятичное число, представляющее количество 128-и битных слов, из которых состоят вещественные числа размещенные в массивах: размер чисел в битах кодировка -------------------- --------- 128 001 256 002 512 004 1024 008 2048 016 4096 032 8192 064 16384 128 32768 256 65536 512 (это предел); ПРИМЕРЫ: ------- ei_log_001_emu32lsb.dfn - 128-бит, ei_log_512_emu32lsb.dfn - 65536-бит. ***************************************************************/ #if MPU_MATH_FN_LIMIT >= 128 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 128 */ #if MPU_MATH_FN_LIMIT >= 256 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 256 */ #if MPU_MATH_FN_LIMIT >= 512 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 512 */ #if MPU_MATH_FN_LIMIT >= 1024 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 1024 */ #if MPU_MATH_FN_LIMIT >= 2048 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 2048 */ #if MPU_MATH_FN_LIMIT >= 4096 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 4096 */ #if MPU_MATH_FN_LIMIT >= 8192 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 8192 */ #if MPU_MATH_FN_LIMIT >= 16384 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 16384 */ #if MPU_MATH_FN_LIMIT >= 32768 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 32768 */ #if MPU_MATH_FN_LIMIT >= 65536 #if MPU_WORD_ORDER_BIG_ENDIAN == 0 #include #else #include #endif #endif /* MPU_MATH_FN_LIMIT >= 65536 */ static int _get_n_log( 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_L128; break; #endif /* MPU_MATH_FN_LIMIT >= 128 */ #if MPU_MATH_FN_LIMIT >= 256 case NBR_256 : rc = N_LOG_L256; break; #endif /* MPU_MATH_FN_LIMIT >= 256 */ #if MPU_MATH_FN_LIMIT >= 512 case NBR_512 : rc = N_LOG_L512; break; #endif /* MPU_MATH_FN_LIMIT >= 512 */ #if MPU_MATH_FN_LIMIT >= 1024 case NBR_1024 : rc = N_LOG_L1024; break; #endif /* MPU_MATH_FN_LIMIT >= 1024 */ #if MPU_MATH_FN_LIMIT >= 2048 case NBR_2048 : rc = N_LOG_L2048; break; #endif /* MPU_MATH_FN_LIMIT >= 2048 */ #if MPU_MATH_FN_LIMIT >= 4096 case NBR_4096 : rc = N_LOG_L4096; break; #endif /* MPU_MATH_FN_LIMIT >= 4096 */ #if MPU_MATH_FN_LIMIT >= 8192 case NBR_8192 : rc = N_LOG_L8192; break; #endif /* MPU_MATH_FN_LIMIT >= 8192 */ #if MPU_MATH_FN_LIMIT >= 16384 case NBR_16384: rc = N_LOG_L16384; 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_log() */ static EMUSHORT *_get_log__L_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__L_128_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 128 */ #if MPU_MATH_FN_LIMIT >= 256 case NBR_256 : rc = (EMUSHORT *)&_ei_log__L_256_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 256 */ #if MPU_MATH_FN_LIMIT >= 512 case NBR_512 : rc = (EMUSHORT *)&_ei_log__L_512_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 512 */ #if MPU_MATH_FN_LIMIT >= 1024 case NBR_1024 : rc = (EMUSHORT *)&_ei_log__L_1024_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 1024 */ #if MPU_MATH_FN_LIMIT >= 2048 case NBR_2048 : rc = (EMUSHORT *)&_ei_log__L_2048_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 2048 */ #if MPU_MATH_FN_LIMIT >= 4096 case NBR_4096 : rc = (EMUSHORT *)&_ei_log__L_4096_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 4096 */ #if MPU_MATH_FN_LIMIT >= 8192 case NBR_8192 : rc = (EMUSHORT *)&_ei_log__L_8192_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 8192 */ #if MPU_MATH_FN_LIMIT >= 16384 case NBR_16384: rc = (EMUSHORT *)&_ei_log__L_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__L_ptr() */ void ei_log__L( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_log__L() Работает с internal e-type data struct. log( 1+x ) - 2S Concepts : RETURN in EIY = ----------------- S x where EIX = S*S, S = -------, 2 + x 0 <= EIX <= 0.0294372515228594143797... = = (sqrt(2)/(2+sqrt(2)))^4. METHOD : Берем S = x/(2+x) и основываясь на том, что log(1+x) = log(1+S) - log(1-S) = = 2*S + (2/3)*S^3 + (2/5)*S^5 + ... = L сможем рассчитать величину (log(1+x)-2*S)/S как L - 2*S --------- = S 2*S + (2/3)*S^3 + (2/5)*S^5 + ... - 2*S = ----------------------------------------- = S = (2/3)*S^2 + (2/5)*S^4 + (2/7)*S^6 + ... . Положив Z = S*S имеем log__L(Z) = = (2/3)*Z + (2/5)*Z^2 + (2/7)*Z^3 + ... = = L1*Z + L2*Z^2 + L3*Z^3 + L4*Z^4 + ... . Create a polinomial approximation to log__L(z) = = z*(L1 + z*(L2 + z*(L3 + ... )) ... ). 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_log( 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__L_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__L() */ /*************************************************************** Hide internal symbols: ***************************************************************/ __mpu_hidden_decl(ei_log__L); /* End of hide internal symbols. ***************************************************************/