/*************************************************************** __ST_SIN.C This file contains source code of functions for SIN 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_sin_001_emu32lsb.dfn - 128-бит, ei_sin_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_sin( 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_SIN_S128; break; #endif /* MPU_MATH_FN_LIMIT >= 128 */ #if MPU_MATH_FN_LIMIT >= 256 case NBR_256 : rc = N_SIN_S256; break; #endif /* MPU_MATH_FN_LIMIT >= 256 */ #if MPU_MATH_FN_LIMIT >= 512 case NBR_512 : rc = N_SIN_S512; break; #endif /* MPU_MATH_FN_LIMIT >= 512 */ #if MPU_MATH_FN_LIMIT >= 1024 case NBR_1024 : rc = N_SIN_S1024; break; #endif /* MPU_MATH_FN_LIMIT >= 1024 */ #if MPU_MATH_FN_LIMIT >= 2048 case NBR_2048 : rc = N_SIN_S2048; break; #endif /* MPU_MATH_FN_LIMIT >= 2048 */ #if MPU_MATH_FN_LIMIT >= 4096 case NBR_4096 : rc = N_SIN_S4096; break; #endif /* MPU_MATH_FN_LIMIT >= 4096 */ #if MPU_MATH_FN_LIMIT >= 8192 case NBR_8192 : rc = N_SIN_S8192; break; #endif /* MPU_MATH_FN_LIMIT >= 8192 */ #if MPU_MATH_FN_LIMIT >= 16384 case NBR_16384: rc = N_SIN_S16384; 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_sin() */ static EMUSHORT *_get_sin__S_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_sin__S_128_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 128 */ #if MPU_MATH_FN_LIMIT >= 256 case NBR_256 : rc = (EMUSHORT *)&_ei_sin__S_256_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 256 */ #if MPU_MATH_FN_LIMIT >= 512 case NBR_512 : rc = (EMUSHORT *)&_ei_sin__S_512_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 512 */ #if MPU_MATH_FN_LIMIT >= 1024 case NBR_1024 : rc = (EMUSHORT *)&_ei_sin__S_1024_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 1024 */ #if MPU_MATH_FN_LIMIT >= 2048 case NBR_2048 : rc = (EMUSHORT *)&_ei_sin__S_2048_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 2048 */ #if MPU_MATH_FN_LIMIT >= 4096 case NBR_4096 : rc = (EMUSHORT *)&_ei_sin__S_4096_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 4096 */ #if MPU_MATH_FN_LIMIT >= 8192 case NBR_8192 : rc = (EMUSHORT *)&_ei_sin__S_8192_[0][0]; break; #endif /* MPU_MATH_FN_LIMIT >= 8192 */ #if MPU_MATH_FN_LIMIT >= 16384 case NBR_16384: rc = (EMUSHORT *)&_ei_sin__S_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_sin__S_ptr() */ void ei_sin__S( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_sin__S() Работает с internal e-type data struct. sin( EIX*k ) - EIX Concepts : RETURN in EIY = -------------------- EIX 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 (sin(k*x)-x)/x = z*(S0 + S1*z^1 + S2*z^2 + ... ). Then ei_sin__S(x*x) = z*(S0 + S1*z^1 + S2*z^2 + ... ). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. NOTE : The coefficient S'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_sin( 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_sin__S_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_sin__S() */ /*************************************************************** Hide internal symbols: ***************************************************************/ __mpu_hidden_decl(ei_sin__S); /* End of hide internal symbols. ***************************************************************/