/*************************************************************** __MPU_MATH.C This file contains source code of functions for REAL MATH 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 #include #include #include #include #include #include #include #include #include #include #include #include #include static void ei_copysign( EMUSHORT *eiy, EMUSHORT *eix, EMUSHORT *eis, int nb ) /*************************************************************** Description : ei_copysign() Работает с internal e-type data struct. Concepts : Copy VALUE from EIX into EIY and SIGN from EIS into EIY. NOTE : Use Global Variable: Use Functions : internal_np( nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-real.c Parameters : EMUSHORT *eiy; - указатель на internal e-type data struct. TARGET; EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE VALUE; EMUSHORT *eis; - указатель на internal e-type data struct. SOURCE SIGN; int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT sign; int ps = 0; #if MPU_WORD_ORDER_BIG_ENDIAN == 0 ps = internal_np( nb ) - 1; #endif sign = eis[ps]; ei_copy( eiy, eix, nb ); eiy[ps] = sign; } /* End of ei_copysign() */ static void ei_drem( EMUSHORT *eiy, EMUSHORT *eix, EMUSHORT *eid, int nb ) /*************************************************************** Description : ei_drem() Работает с internal e-type data struct. Concepts : Reduce EIX into [-EID/2, EID/2] and return result in EIY. NOTE : Use Global Variable: Use Functions : internal_np( nb ); | mpu-real.c _gen_two( eic, nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-real.c ei_remain( r,q,y,x,nb ); | mpu-real.c Parameters : EMUSHORT *eiy; - указатель на internal e-type data struct. TARGET; EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE VALUE; EMUSHORT *eid; - указатель на internal e-type data struct. SOURCE + 2*DIAPASON; int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT *x = NULL, *b = NULL, *d = NULL, *two = NULL; EMUSHORT sign; int ps = 0; int np; np = internal_np( nb ); #if MPU_WORD_ORDER_BIG_ENDIAN == 1 ps = np - 1; #endif /*** Allocate memory for x, b, d, two . *********************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } b = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !b ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } d = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !d ) { /* fatal error */ /* FREE x *****************/ /* FREE b *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE b *****************/ /* FREE d *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_two( two, nb ); sign = eix[ps]; ei_copy( x, eix, nb ); x[ps] = (EMUSHORT)0; /* ei_abs( x, nb ); */ ei_copy( d, eid, nb ); d[ps] = (EMUSHORT)0; /* ei_abs( d, nb ); */ ei_div( b, d, two, nb ); if( ei_cmp( x, d, nb ) > 0 ) ei_remain( x, (EMUSHORT *)0, x, d, nb ); if( ei_cmp( x, b, nb ) >= 0 ) { ei_sub( x, x, d, nb ); } x[ps] ^= sign; ei_copy( eiy, x, nb ); /* return value */ /* FREE x *****************/ /* FREE b *****************/ /* FREE d *****************/ /* FREE two ***************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_drem() */ #if BITS_PER_EMUSHORT == 16 static EMUSHORT bzmask[] = { 0xffff, 0xfffe, 0xfffc, 0xfff8, 0xfff0, 0xffe0, 0xffc0, 0xff80, 0xff00, 0xfe00, 0xfc00, 0xf800, 0xf000, 0xe000, 0xc000, 0x8000, 0x0000, }; #else /* not (BITS_PER_EMUSHORT == 16) */ #if BITS_PER_EMUSHORT == 32 static EMUSHORT bzmask[] = { 0xffffffff, 0xfffffffe, 0xfffffffc, 0xfffffff8, 0xfffffff0, 0xffffffe0, 0xffffffc0, 0xffffff80, 0xffffff00, 0xfffffe00, 0xfffffc00, 0xfffff800, 0xfffff000, 0xffffe000, 0xffffc000, 0xffff8000, 0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000, 0xfff00000, 0xffe00000, 0xffc00000, 0xff800000, 0xff000000, 0xfe000000, 0xfc000000, 0xf8000000, 0xf0000000, 0xe0000000, 0xc0000000, 0x80000000, 0x00000000, }; #else /* not (BITS_PER_EMUSHORT == 32) */ #if BITS_PER_EMUSHORT == 64 static EMUSHORT bzmask[] = { 0xffffffffffffffff, 0xfffffffffffffffe, 0xfffffffffffffffc, 0xfffffffffffffff8, 0xfffffffffffffff0, 0xffffffffffffffe0, 0xffffffffffffffc0, 0xffffffffffffff80, 0xffffffffffffff00, 0xfffffffffffffe00, 0xfffffffffffffc00, 0xfffffffffffff800, 0xfffffffffffff000, 0xffffffffffffe000, 0xffffffffffffc000, 0xffffffffffff8000, 0xffffffffffff0000, 0xfffffffffffe0000, 0xfffffffffffc0000, 0xfffffffffff80000, 0xfffffffffff00000, 0xffffffffffe00000, 0xffffffffffc00000, 0xffffffffff800000, 0xffffffffff000000, 0xfffffffffe000000, 0xfffffffffc000000, 0xfffffffff8000000, 0xfffffffff0000000, 0xffffffffe0000000, 0xffffffffc0000000, 0xffffffff80000000, 0xffffffff00000000, 0xfffffffe00000000, 0xfffffffc00000000, 0xfffffff800000000, 0xfffffff000000000, 0xffffffe000000000, 0xffffffc000000000, 0xffffff8000000000, 0xffffff0000000000, 0xfffffe0000000000, 0xfffffc0000000000, 0xfffff80000000000, 0xfffff00000000000, 0xffffe00000000000, 0xffffc00000000000, 0xffff800000000000, 0xffff000000000000, 0xfffe000000000000, 0xfffc000000000000, 0xfff8000000000000, 0xfff0000000000000, 0xffe0000000000000, 0xffc0000000000000, 0xff80000000000000, 0xff00000000000000, 0xfe00000000000000, 0xfc00000000000000, 0xf800000000000000, 0xf000000000000000, 0xe000000000000000, 0xc000000000000000, 0x8000000000000000, 0x0000000000000000, }; #else /* not (BITS_PER_EMUSHORT == 64) */ #error mpu-math.c: Cannot use that size of EMUSHORT type #endif /* BITS_PER_EMUSHORT == 64 */ #endif /* BITS_PER_EMUSHORT == 32 */ #endif /* BITS_PER_EMUSHORT == 16 */ void ei_trunc( EMUSHORT *eix, unsigned int bz, int nb ) /*************************************************************** Description : ei_trunc() Работает с internal e-type data struct. Concepts : ZERO to BZ last bits of Significand of EIX. NOTE : Use Global Variable: Use Functions : internal_ne( nb ); | mpu-real.c internal_ns( nb ); | mpu-real.c ei_issignull( eix, nb ); | mpu-real.c Parameters : EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE, TARGET; unsigned int bz; - количество обнуляемых младших бит мантиссы в external e-type data struct. int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT *p; int ne, ns; unsigned int k = (unsigned)NSBITS(nb); if( nb < NBR_32 ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"trunc" ); return; } if( bz > k ) { /* error: Invalid number of TRUNC bits(BZ) */ /* количество обнуляемых бит больше чем количество бит мантиссы */ __real_truncate_error(); return; } /* Если на входе знаковый ноль. */ if( ei_issignull( eix, nb ) ) { return; } ne = internal_ne( nb ); ns = internal_ns( nb ); #if MPU_WORD_ORDER_BIG_ENDIAN == 1 p = &eix[ne+ns+2]; /* lgw */ *p-- = (EMUSHORT)0; /* p -> low part of Significand */ #else p = &eix[0]; /* lgw */ *p++ = (EMUSHORT)0; /* p -> low part of Significand */ #endif while( bz >= BITS_PER_EMUSHORT ) { #if MPU_WORD_ORDER_BIG_ENDIAN == 1 *p-- = (EMUSHORT)0; #else *p++ = (EMUSHORT)0; #endif bz -= BITS_PER_EMUSHORT; } /* Clear the remaining bits. */ *p &= bzmask[bz]; } /* End of ei_TRUNC() */ static void __ei_TRUNC( EMUSHORT *eix, int nb ) { unsigned int z = 0; switch( nb ) { case NBR_32: case NBR_64: case NBR_128: z = 44; break; case NBR_256: z = 95; break; case NBR_512: z = 184; break; case NBR_1024: z = 388; break; case NBR_2048: z = 770; break; case NBR_4096: z = 1585; break; case NBR_8192: z = 3964; break; case NBR_16384: z = 6424; break; case NBR_32768: z = 12842; break; case NBR_65536: z = 15882; break; case NBR_131072: z = 51758; break; default: /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"TRUNC" ); break; } /* End of switch( nb ) */ ei_trunc( eix, z, nb ); } /* End of __ei_TRUNC() */ void ei_sin( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_sin() Работает с internal e-type data struct. Concepts : Return in EIY = SIN( EIX ). NOTE : See: ei_sin__S() & ei_cos__C(). Use Global Variable: Use Functions : ei_copy( eiy, eix, nb ); | mpu-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 *a = NULL, *c = NULL, *z = NULL, *x = NULL, *half = NULL, *one = NULL, *small = NULL; EMUSHORT *pi, *pi2, *pi_2, *pi_4, *pi3_4; int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"sin" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* SIN(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* SIN(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* SIN(Infinity) must by NaN */ if( ei_isinfin( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL Infinity - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } pi = _get_m_pi_ptr ( nb ); pi2 = _get_m_2pi_ptr ( nb ); pi_2 = _get_m_pi_2_ptr ( nb ); pi_4 = _get_m_pi_4_ptr ( nb ); pi3_4 = _get_m_3pi_4_ptr( nb ); /*** Allocate memory for a, c, z, half, one, small . ******/ a = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !a ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE half **************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); /********************************************* small: 1.0 + small^2 == 1.0; big : big = 1/(small^2); *********************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); ei_drem( x, eix, pi2, nb ); /* reduse EIX into [-PI, PI]. */ ei_copysign( a, x, one, nb ); if( ei_cmp( a, pi_4, nb ) >= 0 ) { if( ei_cmp( a, pi3_4, nb ) >= 0 ) /* ... in [3PI/4, PI]. */ { ei_sub( a, pi, a, nb ); ei_copysign( x, a, x, nb ); } else /* ... in [PI/4, 3PI/4]. */ { /* return: sign(x)*COS(PI/2-|x|).*/ ei_sub( a, pi_2, a, nb ); ei_mul( z, a, a, nb ); ei_cos__C( c, z, nb ); ei_mul( z, z, half, nb ); ei_sub( a, z, c, nb ); ei_sub( a, one, a, nb ); ei_copysign( eiy, a, x, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* return: SIN(x) */ if( ei_cmp( a, small, nb ) < 0 ) { ei_copy( eiy, x, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_mul( z, x, x, nb ); ei_sin__S( z, z, nb ); ei_mul( z, z, x, nb ); ei_add( eiy, x, z, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_sin() */ void ei_cos( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_cos() Работает с internal e-type data struct. Concepts : Return in EIY = COS( EIX ). NOTE : See: ei_sin__S() & ei_cos__C(). Use Global Variable: Use Functions : ei_copy( eiy, eix, nb ); | mpu-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 *a = NULL, *c = NULL, *z = NULL, *x = NULL, *s = NULL, *half = NULL, *one = NULL, *small = NULL; EMUSHORT *pi, *pi2, *pi_2, *pi_4, *pi3_4; int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"cos" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* COS(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"cos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* COS(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"cos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* COS(Infinity) must by NaN */ if( ei_isinfin( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"cos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL Infinity - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } pi = _get_m_pi_ptr ( nb ); pi2 = _get_m_2pi_ptr ( nb ); pi_2 = _get_m_pi_2_ptr ( nb ); pi_4 = _get_m_pi_4_ptr ( nb ); pi3_4 = _get_m_3pi_4_ptr( nb ); /*** Allocate memory for a, c, z, s, half, one, small . *****/ a = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !a ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } s = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE half **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); _gen_one ( s, nb ); /********************************************* small: 1.0 + small^2 == 1.0; big : big = 1/(small^2); *********************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); ei_drem( x, eix, pi2, nb ); /* reduse EIX into [-PI, PI]. */ ei_copysign( a, x, one, nb ); if( ei_cmp( a, pi_4, nb ) >= 0 ) { if( ei_cmp( a, pi3_4, nb ) >= 0 ) /* ... in [3PI/4, PI]. */ { ei_sub( a, pi, a, nb ); ei_neg( s, nb ); } else /* ... in [PI/4, 3PI/4]. */ { /* return: SIN(PI/2-|x|).*/ ei_sub( a, pi_2, a, nb ); ei_mul( z, a, a, nb ); ei_sin__S( z, z, nb ); ei_mul( z, z, a, nb ); ei_add( eiy, a, z, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* return: s*COS(a) */ if( ei_cmp( a, small, nb ) < 0 ) { ei_copy( eiy, s, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_mul( z, a, a, nb ); ei_cos__C( c, z, nb ); ei_mul( z, z, half, nb ); ei_sub( a, z, c, nb ); ei_sub( a, one, a, nb ); ei_copysign( eiy, a, s, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_cos() */ void ei_tan( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_tan() Работает с internal e-type data struct. Concepts : Return in EIY = TAN( EIX ). NOTE : See: ei_sin__S() & ei_cos__C(). Use Global Variable: Use Functions : ei_copy( eiy, eix, nb ); | mpu-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 *a = NULL, *c = NULL, *z = NULL, *ss = NULL, *cc = NULL, *x = NULL, *half = NULL, *one = NULL, *small = NULL; EMUSHORT *pi, *pi_2, *pi_4; int np; int k; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"tan" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* TAN(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"tan", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* TAN(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"tan", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* TAN(Infinity) must by NaN */ if( ei_isinfin( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"tan", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL Infinity - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } pi = _get_m_pi_ptr ( nb ); pi_2 = _get_m_pi_2_ptr ( nb ); pi_4 = _get_m_pi_4_ptr ( nb ); /*** Allocate memory for a, c, z, ss, cc, half, one, small. */ a = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !a ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ss = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ss ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } cc = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !cc ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ /* FREE half **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); /********************************************* small: 1.0 + small^2 == 1.0; big : big = 1/(small^2); *********************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); ei_drem( x, eix, pi, nb ); /* reduse EIX into [-PI/2, PI/2]. */ ei_copysign( a, x, one, nb ); /* ... = abs(x) */ if( ei_cmp( a, pi_4, nb ) >= 0 ) { k = 1; ei_sub( a, pi_2, a, nb ); ei_copysign( x, a, x, nb ); } else { k = 0; if( ei_cmp( a, small, nb ) < 0 ) { ei_copy( eiy, x, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } ei_mul( z, x, x, nb ); ei_cos__C( cc, z, nb ); ei_sin__S( ss, z, nb ); ei_mul( z, z, half, nb ); /* Next get c = cos(x) accurately */ ei_sub( c, z, cc, nb ); ei_sub( c, one, c, nb ); if( k == 0 ) { /* a use as temp variable */ /* return x+(x*(z-(cc-ss)))/c; ... sin/cos */ ei_sub( a, cc, ss, nb ); ei_sub( a, z, a, nb ); ei_mul( a, x, a, nb ); ei_div( a, a, c, nb ); ei_add( a, x, a, nb ); ei_copy( eiy, a, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* a use as temp variable */ /* return c/(x+x*ss); ... cos/sin */ ei_mul( a, x, ss, nb ); ei_add( a, x, a, nb ); ei_div( a, c, a, nb ); ei_copy( eiy, a, nb ); /* FREE x *****************/ /* FREE a *****************/ /* FREE c *****************/ /* FREE z *****************/ /* FREE ss ****************/ /* FREE cc ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( k == 0 ) */ } /* End of ei_tan() */ void ei_log1p( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_log1p() Работает с internal e-type data struct. Concepts : Return in EIY the LOGARITHM of ( 1 + EIX ). METHOD : 1. Argument Reduction: find k and f such that 1+x = 2^k * (1+f), whwre sqrt(2)/2 < 1+f < sqrt(2). 2. Let s = f/(2+f); Based on log(1+f) = log(1+s) - log(1-s) = = 2*S + (2/3)*S^3 + (2/5)*S^5 + ... = L, log(1+f) is computed by log(1+f) = 2*s + s*log__L(s*s), where log__L(z) = = z*(L1 + z*(L2 + z*(L3 + ... )) ... ). See: ei_log__L() for the values of the coefficients. 3. Finaly, log(1+x) = k*ln2 + log(1+f). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. NOTE : 1. В шаге 3 n*ln2 всегда сохраняется в двух числах n*ln2hi + n*ln2lo, где ln2hi ищется таким, чтобы последние биты были равны 0. Это обеспечивает точное представление произведения n*ln2hi в разрядной сетке машины. Количество младших обнуляемых бит можно увидеть в './service/emuXXXXX/SxxxCxxX.c'. 2. In step 1, f may not representable. A correction term c for f is computed. It follows that the correction term for f - t (the leading term of log(1+f) in step 2) is c-c*x. We add this correction term to n*ln2lo to attenuate the error. SPECIAL CASES : log1p(x) = -InD с выставлением DOMAIN flag [если( x < -1 )]; log1p(+NaN) = +NaN с выставлением DOMAIN flag; log1p(-NaN) = -NaN с выставлением DOMAIN flag; log1p(-inf) = -InD с выставлением DOMAIN flag; log1p(-1) = -inf с выставлением DOMAIN flag; log1p(+inf) = +inf [норма]; log1p(0) = 0 [норма]; Use Global Variable: Use Functions : ei_copy( eiy, eix, nb ); | mpu-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 *z = NULL, *s = NULL, *t = NULL, *c = NULL, *x = NULL, *tx = NULL, /* temp */ *tr = NULL, /* temp */ *half = NULL, *one = NULL, *negone = NULL, *two = NULL, *zero = NULL; static EMUSHORT *k = NULL, /* for Exponent */ *tk = NULL; /* temp */ EMUSHORT *ln2hi, *ln2lo, *sqrt2, *small; int np, ne; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"log1p" ); ei_ind( eiy, nb ); return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); ne = internal_ne( nb ) + 1; /*** Allocate memory for x, zero . **************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE zero **************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); /* Befor test +/-inf */ ei_copy( x, eix, nb ); /*************************** Test for EIX. ***************************/ /* LOG1P(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log1p", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /*Set REAL InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* LOG1P(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log1p", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL NaN - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* LOG1P(+Infinity) must by Infinity; LOG1P(-Infinity) must by InD; */ if( ei_isinfin( eix, nb ) ) { if( ei_cmp( eix, zero, nb ) < 0 ) /* ( EIX < 0 ) */ { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log1p", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL -Infinity - produsing Domain Flag */ } else ei_copy( eiy, eix, nb ); /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for z, s, t, c, tx, tr . ***************/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } s = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for half, one, negone, two . ***********/ half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } negone = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !negone ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for k, tk . ****************************/ k = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !k ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tk = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !tk ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); _gen_one ( negone, nb ); ei_neg ( negone, nb ); _gen_two ( two, nb ); ln2hi = _get_m_ln2hi_ptr( nb ); /* Service Constant */ ln2lo = _get_m_ln2lo_ptr( nb ); /* Service Constant */ sqrt2 = _get_m_sqrt2_ptr( nb ); /* Math Constant */ small = _get_epsilon_ptr( nb ); /* See: FLOATP.H */ if( ei_cmp( x, negone, nb ) > 0 ) /* ( x > -1.0 ) */ { /******************** Argument reduction ********************/ ei_copysign( tx, x, one, nb ); if( ei_cmp( tx, small, nb ) < 0 ) /* if( x < small ) */ { ei_copy( eiy, x, nb ); /* return( x ); */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_add( tx, one, x, nb ); ei_logb( k, tx, ne, nb ); /* k=logb(1+x); */ ei_nege( tk, k, ne ); ei_ldexp( z, tk, x, ne, nb ); /* z=ldexp( x,-k); */ ei_ldexp( t, tk, one, ne, nb ); /* t=ldexp(one,-k); */ ei_add( tx, z, t, nb ); if( ei_cmp( tx, sqrt2, nb ) >= 0 ) /* if( z+t >= sqrt2 ) */ { ei_ince( k, k, ne ); /* k += 1; */ ei_mul( z, z, half, nb ); /* z *= half; */ ei_mul( t, t, half, nb ); /* t *= half; */ } ei_add( t, t, negone, nb ); /* t += negone; */ ei_add( x, z, t, nb ); /* x = z + t; */ ei_sub( tx, t, x, nb ); /* Correction term for x */ ei_add( c, tx, z, nb ); /* c = (t-x) + z; */ /****************** Compute log(1+x) ******************/ ei_add( tx, two, x, nb ); ei_div( s, x, tx, nb ); /* s = x/(2+x); */ ei_mul( tx, x, x, nb ); ei_mul( t, tx, half, nb ); /* t = x*x*half; */ ei_ltor( tr, k, nb, ne ); ei_mul( tr, tr, ln2lo, nb ); ei_mul( tx, c, x, nb ); ei_sub( tr, tr, tx, nb ); ei_add( c, c, tr, nb ); /* c += (k*ln2lo-c*x); */ ei_mul( tr, s, s, nb ); ei_log__L( tx, tr, nb ); ei_add( tx, t, tx, nb ); ei_mul( tr, s, tx, nb ); ei_add( z, c, tr, nb ); /* z = c+s*(t+log_L(s*s)); */ ei_sub( tx, z, t, nb ); ei_add( x, x, tx, nb ); /* x += (z - t); */ ei_ltor( tr, k, nb, ne ); ei_mul( tr, tr, ln2hi, nb ); ei_add( eiy, tr, x, nb ); /* return( k*ln2hi+x ); */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else /* т.е. ( x <= -1.0 ) */ { if( ei_cmp( x, negone, nb ) == 0 ) /* ( x == -1.0 ) */ { /* x == -1, return -Infinity with signal */ /* "argument domain error" */ /* return: -Infinity */ ei_infin( eiy, (unsigned)1, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log1p", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* ( x == -1 ) - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* x < -1, return -InD with signal */ /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log1p", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* ( x < -1 ) - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ /* FREE z *****************/ /* FREE s *****************/ /* FREE t *****************/ /* FREE c *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x > -1.0 ) */ } /* End of ei_log1p() */ void ei_log( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_log() Работает с internal e-type data struct. Concepts : Return in EIY the LOGARITHM of ( EIX ). METHOD : * This code was derived, with minor * modification, from: * Teter Tang, "Table-Driven Implementation * of the Logarithm in IEEE Floating-Point * arithmetic." ACM Trans. Math Software, * vol 16. no 4, pp 378-400, Dec 1990). Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, where F = j/128 for j an integer in [0, 128]. log(2^m) = log2_hi*m + log2_tail*m since m is an integer, the dominant term is exact. m has at most nE*HOST_BITS_PER_EMUSHORT digits (for subnormal numbers), and log2_hi has nE*HOST_BITS_PER_EMUSHORT + 1 trailing zero bits. log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table[] logF_hi[] + 512 is exact. log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... the leading term is calculated to extra precision in two parts, the larger of which adds exactly to the dominant m and F terms. There are two cases: 1. when m, j are non-zero (m | j), use absolute precision for the leading term. 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). In this case, use a relative precision of [after __ei_TRUNC()] bits. (This is done differently in the original paper). Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. Used for generation of extend precision logarithms. Values for log(F) were generated using ei_log1p() function. EI_LOG__N() : Построение ряда: ===================================================================== log(1+f/F) = 2*f/(2*F+f) + 1/12 * (2*f/(2*F+f))**3 + ... q = f/F; s = q/(2+q) = (f/F)/(2+f/F) = f/(2F+f); тогда log(1+f/F) = 2s + (2/3)s^3 + (2/5)s^5 + ... = / f \ 2 / f \ ^3 2 / f \ ^5 = 2*| ---- | + ---*| ---- | + ---*| ---- | + ... = \2F+f/ 3 \2F+f/ 5 \2F+f/ / f \ 2*2^2 / f \ ^3 2*2^4 / f \ ^5 = 2*| ---- | + -----*| ---- | + -----*| ---- | + ... = \2F+f/ 3*2^2 \2F+f/ 5*2^4 \2F+f/ / 2f \ 1 / 2f \ ^3 1 / 2f \ ^5 = | ---- | + -----*| ---- | + -----*| ---- | + ... . \2F+f/ 3*2^2 \2F+f/ 5*2^4 \2F+f/ 1 1 ei_log__N(x) = -----*x^3 + -----*x^5 + ... . 3*2^2 5*2^4 Общая формула получения коэффициентов для APPROX(): k = 1/(s*2^(s-1)), где s = 3, 5, 7, 9, 11, ... . Интервал для APPROX(): [0.0, 1/256]. ===================================================================== ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. NOTE : log2_hi = logF_head[n_log_table]; log2_lo = logF_tail[n_log_table]; Not use the service constants ln2hi, ln2lo. SPECIAL CASES : log(x) = -InD с выставлением DOMAIN flag [если( x < 0 )]; log(+NaN) = +NaN с выставлением DOMAIN flag; log(-NaN) = -NaN с выставлением DOMAIN flag; log(-inf) = -InD с выставлением DOMAIN flag; log( 0 ) = -inf с выставлением DOMAIN flag; log(+inf) = +inf [норма]; log(1) = 0 [норма]; Use Global Variable: Use Functions : ei_log__N( eiy, eix, nb ); | st-ln.c ei_copy( eiy, eix, nb ); | mpu-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, *F = NULL, *f = NULL, *g = NULL, *q = NULL, *u = NULL, *u1 = NULL, *u2 = NULL, *v = NULL, *zero = NULL, *one = NULL, *half = NULL, *two = NULL, *tx = NULL, /* temp */ *tm = NULL, /* temp */ *tr = NULL; /* temp */ EMUSHORT *log2_hi = NULL, *log2_lo = NULL; EMUSHORT *m = NULL, /* for Exponent */ *je = NULL, *min_bin_exp = NULL, *mt = NULL; /* temp */ EMUSHORT *p, *logF_head, *logF_tail; __mpu_int32_t j; int np, ne, i; int n_log_table; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"log" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); ne = internal_ne( nb ) + 1; /*** Allocate memory for x, y . *****************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); /* Befor test +/-inf */ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* LOG(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* LOG(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL NaN - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* LOG(+Infinity) must by Infinity; LOG(-Infinity) must by InD; */ if( ei_isinfin( eix, nb ) ) { if( ei_cmp( eix, zero, nb ) < 0 ) /* ( EIX < 0 ) */ { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* Set REAL -Infinity - produsing Domain Flag */ } else ei_copy( eiy, eix, nb ); /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for F, f, g, q, u, u1, u2, v . *********/ F = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !F ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } f = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !f ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } g = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !g ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } q = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !q ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u1 = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u1 ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u2 = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u2 ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } v = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !v ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for one, half, two, tx, tm, tr . *******/ one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tm = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tm ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(14*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for log2_hi, log2_lo . *****************/ log2_hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !log2_hi ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(16*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } log2_lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !log2_lo ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ __mpu_sbrk( -(int)(17*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for m, je, min_bin_exp, mt . ***********/ m = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !m ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } je = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !je ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } min_bin_exp = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !min_bin_exp ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } mt = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !mt ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ __mpu_sbrk( -(int)(3*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ n_log_table = _get_n_log_table( nb ); logF_head = _get_logF_head_ptr( nb ); logF_tail = _get_logF_tail_ptr( nb ); ei_cpye( min_bin_exp, _get_min_2_exp_ptr( nb ), ne, ne ); _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* Service Constant: hight entries of log_table[] */ ei_copy( log2_hi, logF_head + np*n_log_table, nb ); ei_copy( log2_lo, logF_tail + np*n_log_table, nb ); /* ei_copy( x, eix, nb ); */ if( ei_cmp( x, zero, nb ) > 0 ) /* ( x > 0.0 ) */ { /******************** Argument reduction ********************/ ei_logb( m, x, ne, nb ); /* m = logb( x ); */ ei_nege( mt, m, ne ); ei_ldexp( g, mt, x, ne, nb ); /* g = ldexp( x, -m ); */ /* if( m == min_bin_exp ) */ if( ei_cmpe( m, min_bin_exp, ne ) == 0 ) { ei_logb( je, g, ne, nb ); /* je = logb( g ); */ ei_adde( m, m, je, ne ); /* m += je; */ ei_nege( mt, je, ne ); ei_ldexp( g, mt, g, ne, nb ); /* g = ldexp( g, -je ); */ } /* begin: ************************ j = n_log_table*(g-1) + .5; F = (1.0/n_log_table) * j + 1; f = g - F; *********************************/ j = (__mpu_int32_t)n_log_table; ei_ltor( tr, (EMUSHORT *)&j, nb, 1 ); ei_sub( tx, g, one, nb ); ei_mul( tr, tr, tx, nb ); ei_add( tr, tr, half, nb ); /* tr=j; */ ei_ltor( tx, (EMUSHORT *)&j, nb, 1 ); ei_div( tx, one, tx, nb ); /* begin : ********************************* Skip digits after decimal point. *******************************************/ ei_rtol_frac( (EMUSHORT *)&j, /* j = tr; */ (EMUSHORT *)0, tr, 1, nb ); ei_ltor( tr, (EMUSHORT *)&j, nb, 1 ); /* end : *********************************** Skip digits after decimal point. *******************************************/ ei_mul( tx, tx, tr, nb ); ei_add( F, tx, one, nb ); /* F*128 is an integer in [128,512] */ ei_sub( f, g, F, nb ); /* For following Code. */ ei_rtol_frac( je, /* je = tr; */ (EMUSHORT *)0, tr, ne, nb ); /* end: ************************** j = n_log_table*(g-1) + .5; F = (1.0/n_log_table) * j + 1; f = g - F; *********************************/ /************************************************** Approximate expansion for log( 1+f/F ) ~= u + q. **************************************************/ ei_mul( tr, two, F, nb ); ei_add( tr, tr, f, nb ); ei_div( g, one, tr, nb ); /* g = 1/(2*F+f); */ ei_mul( u, two, f, nb ); ei_mul( u, u, g, nb ); /* u = 2*f*g; */ ei_mul( v, u, u, nb ); /* v = u*u; */ ei_log__N( tx, v, nb ); ei_mul( q, u, tx, nb ); /* q = u*v*(A1+v*(A2+v*(A3+...))); */ if( ei_cmp0e( m, ne ) || ei_cmp0e( je, ne ) ) { /* Формирование константы 0x2000.....0001; */ #if MPU_WORD_ORDER_BIG_ENDIAN == 1 /* hight part */ p = mt; *p++ = (EMUSHORT)0; *p++ = HIGHT_EXTWO >> 1; /* 0x2000... */ for( i = 0; i < ne - 2; i++ ) *p++ = (EMUSHORT)0; p--; *p |= 1; #else /* hight part */ p = mt + ne - 1; *p-- = (EMUSHORT)0; *p-- = HIGHT_EXTWO >> 1; for( i = 0; i < ne - 2; i++ ) *p-- = (EMUSHORT)0; p++; *p |= 1; #endif ei_ltor( tr, mt, nb, ne ); ei_add( u1, u, tr, nb ); /* u1 = u + 513; */ ei_sub( u1, u1, tr, nb ); /* u1 -= 513; */ } else { ei_copy( u1, u, nb ); __ei_TRUNC( u1, nb ); } ei_mul( tr, F, u1, nb ); ei_sub( tr, f, tr, nb ); ei_mul( tr, two, tr, nb ); ei_mul( tx, u1, f, nb ); ei_sub( tr, tr, tx, nb ); ei_mul( u2, tr, g, nb ); /* u2 = (2.0*(f - F*u1) - u1*f) * g; */ /****************************************** u1 + u2 = 2f/(2F+f); to extra precision. ******************************************/ /*********************************************************** log(x) = log(2^m*F*(1+f/F)) = (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); (exact) + (tiny). ***********************************************************/ ei_ltor( tm, m, nb, ne ); ei_mul( tr, tm, log2_hi, nb ); ei_add( tr, tr, logF_head+np*j, nb ); ei_add( u1, u1, tr, nb ); /* u1 += m*logF_head[N] + logF_head[j]; */ ei_add( tr, u2, logF_tail+np*j, nb ); ei_add( u2, tr, q, nb ); /* u2 = (u2 + logF_tail[j]) + q; */ ei_mul( tr, log2_lo, tm, nb ); ei_add( u2, u2, tr, nb ); /* u2 += logF_tail[N]*m; */ ei_add( eiy, u1, u2, nb ); /* return( u1 + u2 ); */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ /* FREE mt ****************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else /* т.е. ( x <= 0.0 ) */ { if( ei_cmp( x, zero, nb ) == 0 ) /* ( x == 0.0 ) */ { /* x == 0.0, return -Infinity with signal */ /* "argument domain error" */ /* return: -Infinity */ ei_infin( eiy, (unsigned)1, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* ( x == 0 ) - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ /* FREE mt ****************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* x < 0.0, return -InD with signal */ /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"log", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* ( x < 0 ) - produsing Domain Flag */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ /* FREE mt ****************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x > 0.0 ) */ } /* End of ei_log() */ typedef struct Double Double; struct Double { EMUSHORT *a; EMUSHORT *b; }; static void __ei_log__D( Double diy, EMUSHORT *eix, int nb ) { EMUSHORT *x = NULL, *F = NULL, *f = NULL, *g = NULL, *q = NULL, *u = NULL, *u1 = NULL, *u2 = NULL, *v = NULL, *zero = NULL, *one = NULL, *half = NULL, *two = NULL, *tx = NULL, /* temp */ *tm = NULL, /* temp */ *tr = NULL; /* temp */ EMUSHORT *log2_hi = NULL, *log2_lo = NULL; EMUSHORT *m = NULL, /* for Exponent */ *je = NULL, *min_bin_exp = NULL, *mt = NULL; /* temp */ EMUSHORT *p, *logF_head, *logF_tail; __mpu_int32_t j; int np, ne, i; int n_log_table; 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 */ /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); ne = internal_ne( nb ) + 1; /*** Allocate memory for x, y . *****************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for F, f, g, q, u, u1, u2, v . *********/ F = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !F ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } f = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !f ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } g = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !g ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } q = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !q ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u1 = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u1 ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } u2 = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !u2 ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } v = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !v ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for one, half, two, tx, tm, tr . *******/ one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tm = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tm ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(14*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for log2_hi, log2_lo . *****************/ log2_hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !log2_hi ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(16*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } log2_lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !log2_lo ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ __mpu_sbrk( -(int)(17*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for m, je, min_bin_exp, mt . ***********/ m = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !m ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } je = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !je ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } min_bin_exp = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !min_bin_exp ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } mt = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !mt ) { /* fatal error */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ __mpu_sbrk( -(int)(3*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ n_log_table = _get_n_log_table( nb ); logF_head = _get_logF_head_ptr( nb ); logF_tail = _get_logF_tail_ptr( nb ); ei_cpye( min_bin_exp, _get_min_2_exp_ptr( nb ), ne, ne ); _gen_zero( zero, nb ); _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* Service Constant: hight entries of log_table[] */ ei_copy( log2_hi, logF_head + np*n_log_table, nb ); ei_copy( log2_lo, logF_tail + np*n_log_table, nb ); ei_copy( x, eix, nb ); /******************** Argument reduction ********************/ ei_logb( m, x, ne, nb ); /* m = logb( x ); */ ei_nege( mt, m, ne ); ei_ldexp( g, mt, x, ne, nb ); /* g = ldexp( x, -m ); */ /* if( m == min_bin_exp ) */ if( ei_cmpe( m, min_bin_exp, ne ) == 0 ) { ei_logb( je, g, ne, nb ); /* je = logb( g ); */ ei_adde( m, m, je, ne ); /* m += je; */ ei_nege( mt, je, ne ); ei_ldexp( g, mt, g, ne, nb ); /* g = ldexp( g, -je ); */ } /* begin: ************************ j = n_log_table*(g-1) + .5; F = (1.0/n_log_table) * j + 1; f = g - F; *********************************/ j = (__mpu_int32_t)n_log_table; ei_ltor( tr, (EMUSHORT *)&j, nb, 1 ); ei_sub( tx, g, one, nb ); ei_mul( tr, tr, tx, nb ); ei_add( tr, tr, half, nb ); /* tr=j; */ ei_ltor( tx, (EMUSHORT *)&j, nb, 1 ); ei_div( tx, one, tx, nb ); /* begin : ********************************* Skip digits after decimal point. *******************************************/ ei_rtol_frac( (EMUSHORT *)&j, /* j = tr; */ (EMUSHORT *)0, tr, 1, nb ); ei_ltor( tr, (EMUSHORT *)&j, nb, 1 ); /* end : *********************************** Skip digits after decimal point. *******************************************/ ei_mul( tx, tx, tr, nb ); ei_add( F, tx, one, nb ); /* F*128 is an integer in [128,512] */ ei_sub( f, g, F, nb ); /* For following Code. */ ei_rtol_frac( je, /* je = tr; */ (EMUSHORT *)0, tr, ne, nb ); /* end: ************************** j = n_log_table*(g-1) + .5; F = (1.0/n_log_table) * j + 1; f = g - F; *********************************/ /************************************************** Approximate expansion for log( 1+f/F ) ~= u + q. **************************************************/ ei_mul( tr, two, F, nb ); ei_add( tr, tr, f, nb ); ei_div( g, one, tr, nb ); /* g = 1/(2*F+f); */ ei_mul( u, two, f, nb ); ei_mul( u, u, g, nb ); /* u = 2*f*g; */ ei_mul( v, u, u, nb ); /* v = u*u; */ ei_log__N( tx, v, nb ); ei_mul( q, u, tx, nb ); /* q = u*v*(A1+v*(A2+v*(A3+...))); */ if( ei_cmp0e( m, ne ) || ei_cmp0e( je, ne ) ) { /* Формирование константы 0x2000.....0001; */ #if MPU_WORD_ORDER_BIG_ENDIAN == 1 /* hight part */ p = mt; *p++ = (EMUSHORT)0; *p++ = HIGHT_EXTWO >> 1; /* 0x2000... */ for( i = 0; i < ne - 2; i++ ) *p++ = (EMUSHORT)0; p--; *p |= 1; #else /* hight part */ p = mt + ne - 1; *p-- = (EMUSHORT)0; *p-- = HIGHT_EXTWO >> 1; for( i = 0; i < ne - 2; i++ ) *p-- = (EMUSHORT)0; p++; *p |= 1; #endif ei_ltor( tr, mt, nb, ne ); ei_add( u1, u, tr, nb ); /* u1 = u + 513; */ ei_sub( u1, u1, tr, nb ); /* u1 -= 513; */ } else { ei_copy( u1, u, nb ); __ei_TRUNC( u1, nb ); } ei_mul( tr, F, u1, nb ); ei_sub( tr, f, tr, nb ); ei_mul( tr, two, tr, nb ); ei_mul( tx, u1, f, nb ); ei_sub( tr, tr, tx, nb ); ei_mul( u2, tr, g, nb ); /* u2 = (2.0*(f - F*u1) - u1*f) * g; */ /****************************************** u1 + u2 = 2f/(2F+f); to extra precision. ******************************************/ /*********************************************************** log(x) = log(2^m*F*(1+f/F)) = (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); (exact) + (tiny). ***********************************************************/ ei_ltor( tm, m, nb, ne ); ei_mul( tr, tm, log2_hi, nb ); ei_add( tr, tr, logF_head+np*j, nb ); ei_add( u1, u1, tr, nb ); /* u1 += m*logF_head[N] + logF_head[j]; */ ei_add( tr, u2, logF_tail+np*j, nb ); ei_add( u2, tr, q, nb ); /* u2 = (u2 + logF_tail[j]) + q; */ ei_mul( tr, log2_lo, tm, nb ); ei_add( u2, u2, tr, nb ); /* u2 += logF_tail[N]*m; */ /* Only difference is here */ ei_add( diy.a, u1, u2, nb ); /* diy.a = u1 + u2; */ __ei_TRUNC( diy.a, nb ); ei_sub( tr, u1, diy.a, nb ); ei_add( diy.b, tr, u2, nb ); /* diy.b = (u1 - diy.a) + u2; */ /* FREE x *****************/ /* FREE zero **************/ /* FREE F *****************/ /* FREE f *****************/ /* FREE g *****************/ /* FREE q *****************/ /* FREE u *****************/ /* FREE u1 ****************/ /* FREE u2 ****************/ /* FREE v *****************/ /* FREE one ***************/ /* FREE half **************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tr ****************/ /* FREE log2_hi ***********/ /* FREE log2_lo ***********/ __mpu_sbrk( -(int)(18*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE m *****************/ /* FREE je ****************/ /* FREE min_bin_exp *******/ /* FREE mt ****************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; /* return( diy ); */ } /* End of __ei_log__D() */ void ei_log10( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** log10( e ) = 1/ln( 10. ); Переход от натурального логарифма к десятичному: log10( N ) = log10( e )*ln( N ). ***************************************************************/ /******************************* mtherr("ei_log10"); ПОТОМ !!! *******************************/ { EMUSHORT *m_1_ln10; m_1_ln10 = _get_m_1_ln10_ptr( nb ); ei_log( eiy, eix, nb ); ei_mul( eiy, eiy, m_1_ln10, nb ); } void ei_log2( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Переход от натурального логарифма к логарифму по основанию 2: log2( N ) = log2( e )*ln( N ). ***************************************************************/ /****************************** mtherr("ei_log2"); ПОТОМ !!! ******************************/ { EMUSHORT *m_1_ln2; m_1_ln2 = _get_m_1_ln2_ptr( nb ); ei_log( eiy, eix, nb ); ei_mul( eiy, eiy, m_1_ln2, nb ); } static void __ei_exp__E( EMUSHORT *eiy, EMUSHORT *eix, EMUSHORT *eic, int nb ) /*************************************************************** Description : __ei_exp__E() Работает с internal e-type data struct. Concepts : ASSUMPTION: EIC << EIX SO THAT __ei_TRUNC( EIX+EIC) = EIX. (EIC is the correction term for EIX) ei_exp__E(x,c) RETURNS in EIY: / exp(x+c)-1-x , | EPS < |x| < .3465736... ei_exp__E(x,c) = | \ 0, |x| < EPS; where the constant .3465736... = ln(2)/2. METHOD : 1. Rational approximation. Let r=x+c. Based on 2 * sinh(r/2) exp(r) - 1 = -----------------------, cosh(r/2) - sinh(r/2) __ei_exp__E(r) is computed using x*x (x/2)*W-(Q-(2*P+x*P)) --- + (c + x*[----------------------- + c ]) 2 1 - W where P := p1*x^2 + p2*x^4 + ..., Q := q1*x^2 + q2*x^4 + ..., W := x/2-(Q-x*P), See the listing below for the values of p1,p2,q1,q2,q3. The polynomials P and Q may be regarded as the approximations to ei_sinh__S() and ei_cosh__C() : sinh(r/2) = r/2 + r * P, cosh(r/2) = 1 + Q.) EI_SINH__S() : Построение ряда: ============================================================== Известно, что sinh(x) = x + (1/3!)*x^3 + (1/5!)*x^5 + ... . Мы имеем следующий способ sinh(r/2) = r/2 + r*P; и sinh(r/2) = r/2 + (1/3!)*(r/2)^3 + (1/5!)*(r/2)^5 + ... . следовательно P = (1/3!)*(r^2/2^3) + (1/5!)*(r^4/2^5) + ... = 1 1 1 = ------*r^2 + ------*r^4 + ------*r^6 + ... . 2^3*3! 2^5*5! 2^7*7! z = r*r; P = p1*z + p2*z^2 + p3*z^3 + ... . Вот все, что нужно для APPROX(), чтобы получить p1, p2, p3, ... на интервале [0, ln(2)/2] по флгоритму Ремеза. ============================================================== EI_COSH__C() : Построение ряда: ============================================================== Известно, что cosh(x) = 1 + (1/2!)*x^2 + (1/4!)*x^4 + ... . Мы имеем следующий способ cosh(r/2) = 1 + Q; и cosh(r/2) = 1 + (1/2!)*(r/2)^2 + (1/4!)*(r/2)^4 + ... . следовательно 1 1 1 Q = ------*r^2 + ------*r^4 + ------*r^6 + ... . 2^2*2! 2^4*4! 2^6*6! z = r*r; Q = p1*z + p2*z^2 + p3*z^3 + ... . Вот все, что нужно для APPROX(), чтобы получить q1, q2, q3, ... на интервале [0, ln(2)/2] по флгоритму Ремеза. ============================================================== ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. NOTE : The coefficient hS's & hC's [ in __ei_sinh__S(), __ei_cosh__C()] 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 : ei_sinh__S( eiy, eix, nb ); | st-sinh.c ei_cosh__C( eiy, eix, nb ); | st-cosh.c ei_copy( eiy, eix, nb ); | mpu-real.c . . . Parameters : EMUSHORT *eiy; - указатель на internal e-type data struct. TARGET; EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE; EMUSHORT *eic; - указатель на internal e-type data struct. (EIC is the correction term for EIX) SOURCE & TARGET; int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT *z = NULL, *c = NULL, *p = NULL, *q = NULL, *xp = NULL, *xh = NULL, *w = NULL, *x = NULL, *tx = NULL, /* temp */ *tr = NULL, /* temp */ *half = NULL, *one = NULL, *zero = NULL; EMUSHORT *small; int np; 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 */ /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for z, c, p, q, xp, xh, w, x, tx, tr . */ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE z *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } p = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !p ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } q = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !q ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } xp = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !xp ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } xh = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !xh ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } w = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !w ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for half, one, zero . ******************/ half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); _gen_half( half, nb ); _gen_one ( one, nb ); small = _get_epsilon_ptr( nb ); /* See: FLOATP.H */ ei_copy( x, eix, nb ); ei_copy( c, eic, nb ); ei_copysign( tx, x, one, nb ); if( ei_cmp( tx, small, nb ) > 0 ) /* if( x > small ) */ { ei_mul( z, x, x, nb ); ei_sinh__S( p, z, nb ); ei_cosh__C( q, z, nb ); ei_mul( xp, x, p, nb ); /* xp = x*p; */ ei_mul( xh, x, half, nb ); /* xh = x*half; */ ei_sub( tx, q, xp, nb ); ei_sub( w, xh, tx, nb ); /* w = xh-(q-xp); */ ei_add( p, p, p, nb ); /* p = p + p; */ ei_add( tx, p, xp, nb ); ei_sub( tx, q, tx, nb ); ei_mul( tr, xh, w, nb ); ei_sub( tx, tr, tx, nb ); ei_sub( tr, one, w, nb ); ei_div( tx, tx, tr, nb ); ei_add( tx, tx, c, nb ); ei_mul( tx, x, tx, nb ); ei_add( c, c, tx, nb ); /* c+=x*[(xh*w-(q-(p+xp)))/(one-w)+c]; */ ei_mul( tr, z, half, nb ); ei_add( tr, tr, c, nb ); /* return( z*half+c ); */ ei_copy( eiy, tr, nb ); /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE zero **************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else /* т.е. ( |x| < small ) */ { if( ei_cmp( x, zero, nb ) != 0 ) { /* raise the inexact flag */ } /* return( copysign(zero,x) ); */ ei_copysign( eiy, zero, x, nb ); /* FREE z *****************/ /* FREE c *****************/ /* FREE p *****************/ /* FREE q *****************/ /* FREE xp ****************/ /* FREE xh ****************/ /* FREE w *****************/ /* FREE x *****************/ /* FREE tx ****************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE zero **************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End of __ei_exp__E() */ void ei_expm1( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_expm1() Работает с internal e-type data struct. Concepts : EI_EXPM1(EIX) RETURN in EIY THE EXPONENTIAL OF EIX MINUS ONE. METHOD : 1. Argument Reduction: given the input EIX, find r and integer k such that EIX = k*ln2 + r, |r| <= 0.5*ln2 . r will be represented as r = z+c for better accuracy. 2. Compute EI_EXPM1(r) = exp(r)-1 by EI_EXPM1(r=z+c) := z + __ei_exp__E(z,c) 3. EI_EXPM1(x) = 2^k*(EI_EXPM1(r) + 1-2^-k). Remarks: 1. When k=1 and z < -0.25, we use the following formula for better accuracy: EI_EXPM1(x) = 2*((z+0.5)+__ei_exp__E(z,c)). 2. To avoid rounding error in 1-2^-k where k is large, we use EI_EXPM1(x) = = 2^k*{[z+(__ei_exp__E(z,c)-2^-k )] + 1} when k > NSBITS(nb). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-float.h]. SPECIAL CASES : expm1(+NaN) = +NaN с выставлением DOMAIN flag; expm1(-NaN) = -NaN с выставлением DOMAIN flag; expm1(-inf) = -1 [норма]; expm1(+inf) = +inf [норма]; expm1(0) = 0 [норма]; Use Global Variable: Use Functions : __ei_exp__E( eiy, eix, nb ); | this file ei_copy( eiy, eix, nb ); | mpu-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, *z = NULL, *hi = NULL, *lo = NULL, *c = NULL, *zero = NULL, *half = NULL, *one = NULL, *negone = NULL, *quart = NULL, *tx = NULL, /* temp */ *tm = NULL, /* temp */ *tk = NULL; /* temp */ EMUSHORT *k = NULL, /* for Exponent */ *mk = NULL, *prec1 = NULL, *prec2 = NULL, *kone = NULL; /* temp 1 */ EMUSHORT *ln2hi, *ln2lo, *invln2, *lnhuge, *lntiny; EMUSHORT prec; int np, ne, i; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"expm1" ); /* ei_ind( eiy, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* EXPM1(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"expm1", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* EXPM1(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"expm1", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for z, hi, lo, c . *******************/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !hi ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !lo ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for zero, half, one, negone, quart . ***/ zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } negone = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !negone ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } quart = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !quart ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for tx, tm, tk . ***********************/ tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tm = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tm ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tk = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tk ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for k, mk, prec1, prec2, kone . ********/ k = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !k ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } mk = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !mk ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } prec1 = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !prec1 ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } prec2 = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !prec2 ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ __mpu_sbrk( -(int)(3*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } kone = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !kone ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /* _gen( kone = 1 ); */ for( i = 0; i < ne; i++ ) kone[i] = (EMUSHORT)0; #if MPU_WORD_ORDER_BIG_ENDIAN == 1 kone[ne-1] = (EMUSHORT)1; #else kone[0] = (EMUSHORT)1; #endif prec = NSBITS(nb); ei_cpye_unpack( prec1, (EMUSHORT *)&prec, ne, 1 ); ei_shln( prec2, prec1, ne, ne ); ln2hi = _get_m_ln2hi_ptr( nb ); /* Service Constant */ ln2lo = _get_m_ln2lo_ptr( nb ); /* Service Constant */ lnhuge = _get_m_ln_huge_ptr( nb ); /* Service Constant */ lntiny = _get_m_ln_tiny_ptr( nb ); /* Service Constant */ invln2 = _get_m_1_ln2_ptr( nb ); /* Math Constant */ _gen_zero( zero, nb ); _gen_half( half, nb ); _gen_one ( one, nb ); _gen_one ( negone, nb ); /* Befor test +/-inf */ ei_neg ( negone, nb ); _gen_two ( quart, nb ); ei_div( quart, half, quart, nb ); ei_neg ( quart, nb ); /* quart = -0.25; */ /* ei_copy( x, eix, nb ); */ if( ei_cmp( x, lnhuge, nb ) <= 0 ) { if( ei_cmp( x, lntiny, nb ) >= 0 ) { /******************** Argument reduction ********************/ ei_copysign( tm, half, x, nb ); ei_mul( tx, invln2, x, nb ); ei_add( tx, tx, tm, nb ); /* begin : ********************************* k = _REALtoINT(x/ln(2)). *******************************************/ ei_rtol_frac( k, /* k = tx; */ (EMUSHORT *)0, tx, ne, nb ); /* k = _RtoINT(x/ln2); */ ei_ltor( tk, k, nb, ne ); /* tk = k; */ /* end : *********************************** k = _REALtoINT(x/ln(2)). *******************************************/ ei_mul( tx, tk, ln2hi, nb ); ei_sub( hi, x, tx, nb ); /* hi = x-k*ln2hi; */ ei_mul( lo, tk, ln2lo, nb ); ei_sub( z, hi, lo, nb ); /* z = hi-(lo=k*ln2lo); */ ei_sub( tm, hi, z, nb ); ei_sub( c, tm, lo, nb ); /* c = (hi-z)-lo; */ if( ei_cmp0e( k, ne ) == 0 ) /* if( k == 0 ) */ { __ei_exp__E( tx, z, c, nb ); ei_add( eiy, z, tx, nb ); /* return(z+__ei_exp__E(z,c)); */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ /* FREE kone **************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( k == 0 ) */ if( ei_cmpe( k, kone, ne ) == 0 ) /* if( k == 1 ) */ { if( ei_cmp( z, quart, nb ) < 0 ) /* if( z < -0.25 ) */ { ei_add( x, z, half, nb ); /* x = z+half; */ __ei_exp__E( tx, z, c, nb ); ei_add( x, x, tx, nb ); /* x += __ei_exp__E(z,c); */ ei_add( eiy, x, x, nb ); /* return( x + x ); */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ /* FREE kone **************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { __ei_exp__E( tx, z, c, nb ); ei_add( z, z, tx, nb ); /* z += __ei_exp__E(z,c); */ ei_add( x, half, z, nb ); /* x = half+z; */ ei_add( eiy, x, x, nb ); /* return( x + x ); */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ /* FREE kone **************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( k == 1 ) */ else { ei_nege( mk, k, ne ); /* mk = -k; */ if( ei_cmpe( k, prec1, ne ) <= 0 ) /* if( k <= prec1 ) */ { ei_ldexp( tx, mk, one, ne, nb ); ei_sub( x, one, tx, nb ); /* x = one-ldexp(one,-k); */ __ei_exp__E( tx, z, c, nb ); ei_add( z, z, tx, nb ); /* z += __ei_exp__E(z,c); */ } else if( ei_cmpe( k, prec2, ne ) < 0 ) /* if( k < prec2 ) */ { /* x = __ei_exp__E(z,c)-ldexp(one,-k); */ __ei_exp__E( tx, z, c, nb ); ei_ldexp( tm, mk, one, ne, nb ); ei_sub( x, tx, tm, nb ); ei_add( x, x, z, nb ); /* x += z; */ ei_copy( z, one, nb ); /* z = one; */ } else { __ei_exp__E( tx, z, c, nb ); ei_add( x, tx, z, nb ); /* x = __ei_exp__E(z,c)+z; */ ei_copy( z, one, nb ); /* z = one; */ } ei_add( tx, x, z, nb ); ei_ldexp( eiy, k, tx, ne, nb ); /* return( ldexp(x+z,k) ); */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ /* FREE kone **************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x >= lntiny ) */ else { /* ( x <= lntiny ) */ /**************************************** EXPM1(-Infinity) must by -1.0 [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* -big# - rounded to -1.0 (inexact) */ /* return: -1.0 ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_copy( eiy, negone, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"expm1", __INEXACT__, eiy, x, (EMUSHORT *)0, nb ); __STINX; /* -big# - produsing Inexact Flag */ } ei_copy( eiy, negone, nb ); } /* End of ( x <= lntiny ) */ } /* End if( x <= lnhuge ) */ else { /* ( x >= lnhuge ) */ /**************************************** EXPM1(+Infinity) must by +Inf [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* +big# - overflow to +Inf */ /* return: +Infinity ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_infin( eiy, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную __error_no следующим образом. __error_no = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"expm1", __INEXACT__, eiy, x, (EMUSHORT *)0, nb ); __STINX; /* +big# - produsing Inexact Flag */ } ei_infin( eiy, (unsigned)0, nb ); } /* End of ( x >= lnhuge ) */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE quart *************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE prec2 *************/ /* FREE kone **************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_expm1() */ void ei_exp( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_exp() Работает с internal e-type data struct. Concepts : EI_EXP(EIX) RETURN in EIY THE EXPONENTIAL OF EIX. METHOD : 1. Argument Reduction: given the input EIX, find r and integer k such that EIX = k*ln2 + r, |r| <= 0.5*ln2 . r will be represented as r = z+c for better accuracy. 2. Compute EXP(r) by EXP(r) = 1 + r + r*R1/(2-R1), where R1 = x - x^2*(p1+x^2*(p2+x^2*( ... ))). 3. EI_EXP(x) = 2^k * EXP(r). EI_EXP__E() : Построение ряда: ============================================================== Известно, что exp(x) = 1 + x + (1/2!)*x^2 + (1/3!)*x^3 + ... . Мы имеем: exp(x) = 1 + x + x*R1/(2-R1); следовательно x*R1/(2-R1) = (1/2!)*x^2 + (1/3!)*x^3 + ... = E. тогда (2-R1)/x*R1 = 1/E, 2/(x*R1) - 1/x = 1/E далее 2/(x*R1) = (x+E)/(E*x), (x*R1)/2 = (E*x)/(x+E), R1 = (2*E*x)/(x*(x+E)). Это значит, что (2/2!)*x^3 + (2/3!)*x^4 + (2/4!)*x^5 + ... R1 = --------------------------------------------------. x^2 + (1/2!)*x^3 + (1/3!)*x^4 + (1/4!)*x^5 + ... Разделив полином числителя на полином знаменателя можно получить формулы вычисления коэффициентов p1, p2, p3, ... полинома R1 = p1*z + p2*z^2 + p3*z^3 + ... , где z = x*x. Вот все, что нужно для APPROX(), чтобы получить p1, p2, p3, ... на интервале [0, ln(2)/2] по флгоритму Ремеза. ============================================================== ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. SPECIAL CASES : exp(+NaN) = +NaN с выставлением DOMAIN flag; exp(-NaN) = -NaN с выставлением DOMAIN flag; exp(-inf) = 0 [норма]; exp(+inf) = +inf [норма]; exp(0) = 1 [норма]; Use Global Variable: Use Functions : ei_exp__E( eiy, eix, nb ); | st-exp.c ei_copy( eiy, eix, nb ); | mpu-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, *z = NULL, *hi = NULL, *lo = NULL, *c = NULL, *zero = NULL, *half = NULL, *one = NULL, *two = NULL, *tx = NULL, /* temp */ *tm = NULL, /* temp */ *tk = NULL; /* temp */ EMUSHORT *k = NULL; /* for Exponent */ EMUSHORT *ln2hi, *ln2lo, *invln2, *lnhuge, *lntiny; int np, ne; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"exp" ); /* ei_ind( eiy, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* EXP(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"exp", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* EXP(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"exp", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for z, hi, lo, c . *******************/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !hi ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !lo ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for zero, half, one, two . *************/ zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for tx, tm, tk . ***********************/ tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tm = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tm ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tk = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tk ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for k . ********************************/ k = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !k ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ ln2hi = _get_m_ln2hi_ptr( nb ); /* Service Constant */ ln2lo = _get_m_ln2lo_ptr( nb ); /* Service Constant */ lnhuge = _get_m_ln_huge_ptr( nb ); /* Service Constant */ lntiny = _get_m_ln_tiny_ptr( nb ); /* Service Constant */ invln2 = _get_m_1_ln2_ptr( nb ); /* Math Constant */ _gen_zero( zero, nb ); /* Befor test +/-inf */ _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* ei_copy( x, eix, nb ); */ if( ei_cmp( x, lnhuge, nb ) <= 0 ) { if( ei_cmp( x, lntiny, nb ) >= 0 ) { /************************************* Argument reduction: x --> x - k*ln2 *************************************/ ei_copysign( tm, half, x, nb ); ei_mul( tx, invln2, x, nb ); ei_add( tx, tx, tm, nb ); /* begin : ********************************* k = _REALtoINT(x/ln(2)). *******************************************/ ei_rtol_frac( k, /* k = tx; */ (EMUSHORT *)0, tx, ne, nb ); /* k = _RtoINT(x/ln2); */ ei_ltor( tk, k, nb, ne ); /* tk = k; */ /* end : *********************************** k = _REALtoINT(x/ln(2)). *******************************************/ /***************************************************** Express x-k*ln2 as hi-lo and let x = hi-lo rounded. *****************************************************/ ei_mul( tx, tk, ln2hi, nb ); ei_sub( hi, x, tx, nb ); /* hi = x-k*ln2hi; */ ei_mul( lo, tk, ln2lo, nb ); ei_sub( x, hi, lo, nb ); /* x = hi-(lo=k*ln2lo); */ /******************************** return( 2^k*[1+x+x*c/(2+c)] ); ********************************/ ei_mul( z, x, x, nb ); ei_exp__E( tx, z, nb ); ei_sub( c, x, tx, nb ); /* c = x - ei_exp__E(z); */ ei_mul( tx, x, c, nb ); ei_sub( tm, two, c, nb ); ei_div( tx, tx, tm, nb ); ei_sub( tx, lo, tx, nb ); ei_sub( tx, hi, tx, nb ); ei_add( tx, one, tx, nb ); /* return( ldexp(1.0+(hi-(lo-(x*c)/(2.0-c))), k) ); */ ei_ldexp( eiy, k, tx, ne, nb ); /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( x >= lntiny ) */ else { /* ( x <= lntiny ) */ /**************************************** EXP(-Infinity) must by 0.0 [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* -big# - underflows to ZERO */ /* return: 0.0 ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_copy( eiy, zero, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"exp", __INEXACT__, eiy, x, (EMUSHORT *)0, nb ); __STINX; /* -big# - produsing Inexact Flag */ } ei_copy( eiy, zero, nb ); } /* End of ( x <= lntiny ) */ } /* End if( x <= lnhuge ) */ else { /* ( x >= lnhuge ) */ /**************************************** EXP(+Infinity) must by +Inf [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* +big# - overflows to +Inf */ /* return: +Infinity ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_infin( eiy, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"exp", __INEXACT__, eiy, x, (EMUSHORT *)0, nb ); __STINX; /* +big# - produsing Inexact Flag */ } ei_infin( eiy, (unsigned)0, nb ); } /* End of ( x >= lnhuge ) */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_exp() */ static void __ei_exp__D( EMUSHORT *eiy, EMUSHORT *eix, EMUSHORT *eic, int nb ) /*************************************************************** Return: EXP( r=x+c ) for |c| < |x| with no overlap (без перекрытия). ***************************************************************/ { EMUSHORT *x = NULL, *z = NULL, *hi = NULL, *lo = NULL, *c = NULL, *zero = NULL, *half = NULL, *one = NULL, *two = NULL, *tx = NULL, /* temp */ *tm = NULL, /* temp */ *tk = NULL; /* temp */ EMUSHORT *k = NULL; /* for Exponent */ EMUSHORT *ln2hi, *ln2lo, *invln2, *lnhuge, *lntiny; int np, ne; 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 */ /* ei_ind( eiy, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* EXP__D(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"__ei_exp__D", __DOMAIN__, eiy, x, eic, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* EXP__D(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"__ei_exp__D", __DOMAIN__, eiy, x, eic, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for z, hi, lo, c . *******************/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !hi ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !lo ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } c = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !c ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for zero, half, one, two . *************/ zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for tx, tm, tk . ***********************/ tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tm = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tm ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tk = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tk ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for k . ********************************/ k = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !k ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ ln2hi = _get_m_ln2hi_ptr( nb ); /* Service Constant */ ln2lo = _get_m_ln2lo_ptr( nb ); /* Service Constant */ lnhuge = _get_m_ln_huge_ptr( nb ); /* Service Constant */ lntiny = _get_m_ln_tiny_ptr( nb ); /* Service Constant */ invln2 = _get_m_1_ln2_ptr( nb ); /* Math Constant */ _gen_zero( zero, nb ); /* Befor test +/-inf */ _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* ei_copy( x, eix, nb ); */ ei_copy( c, eic, nb ); if( ei_cmp( x, lnhuge, nb ) <= 0 ) { if( ei_cmp( x, lntiny, nb ) >= 0 ) { /************************************* Argument reduction: x --> x - k*ln2 *************************************/ ei_copysign( tm, half, x, nb ); ei_mul( tx, invln2, x, nb ); ei_add( tx, tx, tm, nb ); /* begin : ********************************* k = _REALtoINT(x/ln(2)). *******************************************/ ei_rtol_frac( k, /* k = tx; */ (EMUSHORT *)0, tx, ne, nb ); /* k = _RtoINT(x/ln2); */ ei_ltor( tk, k, nb, ne ); /* tk = k; */ /* end : *********************************** k = _REALtoINT(x/ln(2)). *******************************************/ /********************************************************* Express (x+c)-k*ln2 as hi-lo and let x = hi-lo rounded. *********************************************************/ ei_mul( tx, tk, ln2hi, nb ); ei_sub( hi, x, tx, nb ); /* hi = x-k*ln2hi; */ ei_mul( lo, tk, ln2lo, nb ); ei_sub( lo, lo, c, nb ); ei_sub( x, hi, lo, nb ); /* x = hi-(lo=k*ln2lo-c); */ /******************************** return( 2^k*[1+x+x*c/(2+c)] ); ********************************/ ei_mul( z, x, x, nb ); ei_exp__E( tx, z, nb ); ei_sub( c, x, tx, nb ); /* c = x - ei_exp__E(z); */ ei_mul( tx, x, c, nb ); ei_sub( tm, two, c, nb ); ei_div( c, tx, tm, nb ); /* c = (x*c)/(2.0-c); */ ei_sub( tx, lo, c, nb ); ei_sub( tx, hi, tx, nb ); ei_add( tx, one, tx, nb ); /* return( ldexp(1.0 + (hi - (lo - c)), k) ); */ ei_ldexp( eiy, k, tx, ne, nb ); /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( x >= lntiny ) */ else { /* ( x <= lntiny ) */ /**************************************** EXP__D(-Infinity) must by 0.0 [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* -big# - underflows to ZERO */ /* return: 0.0 ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_copy( eiy, zero, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"__ei_exp__D", __INEXACT__, eiy, x, eic, nb ); __STINX; /* -big# - produsing Inexact Flag */ } ei_copy( eiy, zero, nb ); } /* End of ( x <= lntiny ) */ } /* End if( x <= lnhuge ) */ else { /* ( x >= lnhuge ) */ /**************************************** EXP__D(+Infinity) must by +Inf [exact]; ****************************************/ if( !ei_isinfin( eix, nb ) ) { /* raise the inexact flag */ /* +big# - overflows to +Inf */ /* return: +Infinity ( after if(){} body ) */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_infin( eiy, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* Non change */ (__mpu_char8_t *)"__ei_exp__D", __INEXACT__, eiy, x, eic, nb ); __STINX; /* +big# - produsing Inexact Flag */ } ei_infin( eiy, (unsigned)0, nb ); } /* End of ( x >= lnhuge ) */ /* FREE x *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE c *****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE tm ****************/ /* FREE tk ****************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of __ei_exp__D() */ void ei_atan2( EMUSHORT *eic, EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_atan2() Работает с internal e-type data struct. Concepts : ATAN2( EIY,EIX ) RETURN ARG( EIX+iEIY). METHOD : 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). 2. Reduce x to positive by (if x and y are unexceptional): ARG(x+iy) = arctan(y/x) ... if x>0, ARG(x+iy) = pi-arctan[y/(-x)] ... if x<0, 3. According to the integer k=4t+0.25 truncated, t=y/x, the argument is further reduced to one of the following intervals and the arctangent of y/x is evaluated by the corresponding formula: [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) [11/16,19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) EI_ATAN2__A() : Построение ряда: ============================================================== Известно, что atan(x) = x - (1/3)*x^3 + (1/5)*x^5 - (1/7)*x^7 + ... . Мы имеем следующий способ atan(t) = t - t*(z*a1+z^2*a2+z^3*a3...) = t - t*A, где z = t*t. Следовательно a1 = +1/3; a2 = -1/5; a3 = +1/7; a4 = -1/9; ... . Вот все, что нужно для APPROX(), чтобы получить a1, a2, a3, ... на интервале [-7/16, +7/16] по флгоритму Ремеза. ============================================================== ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. Atan2(y,x) returns (PI/pi) * the exact ARG(x+iy) nearly rounded, where pi is rounded value of exact PI. NOTE : The coefficient A's [ in ei_atan2__A()] 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; последний из которых я перевел на "С", затем на язык операций повышенной разрядности. ========================================== SPECIAL CASES : Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). ARG( NAN, (anything) ) is NaN; ARG( (anything), NaN ) is NaN; ARG(+(anything but NaN), +-0) is +-0 ; ARG(-(anything but NaN), +-0) is +-PI; ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; ARG(+INF,+-(anything but INF and NaN) ) is +-0 ; ARG(-INF,+-(anything but INF and NaN) ) is +-PI; ARG(+INF,+-INF ) is +- PI/4; ARG(-INF,+-INF ) is +-3PI/4; ARG( (anything but,0,NaN,and INF),+-INF ) is +-PI/2; Use Global Variable: Use Functions : ei_atan2__A( eiy, eix, nb ); | st-atan2.c ei_copy( eiy, eix, nb ); | mpu-real.c . . . Parameters : EMUSHORT *eic; - указатель на internal e-type data struct. TARGET; EMUSHORT *eiy; - указатель на internal e-type data struct. SOURCE; EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE; int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT *x = NULL, *y = NULL, *z = NULL, *hi = NULL, *lo = NULL, *t = NULL, *signy = NULL, *signx = NULL, *zero = NULL, *one = NULL, *two = NULL, *small = NULL, *big = NULL, *tx = NULL, /* temp */ *ty = NULL; /* temp */ EMUSHORT *k = NULL, /* for Exponent */ *m = NULL, /* for Exponent */ *kt = NULL, *mk = NULL, *prec1 = NULL, *mprec2 = NULL; EMUSHORT *pi, *pi_2, *pi_4, *pi3_4; EMUSHORT *athfhi, *athflo, *at1fhi, *at1flo; EMUSHORT *c39_16, *c1_16; EMUSHORT prec, j; int np, ne; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"atan2" ); /* ei_ind( eiy, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; 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; } /************************************************************/ ei_copy( y, eiy, nb ); /* temp for _mtherr() */ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX, EIY. ***************************/ /* if x or y is InD ATAN2(x,y) must by InD. ******************************************/ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eic, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"atan2", __DOMAIN__, eic, y, x, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_isind( eiy, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eic, eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"atan2", __DOMAIN__, eic, y, x, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* if x or y is NaN ATAN2(x,y) must by NaN. ******************************************/ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eic, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"atan2", __DOMAIN__, eic, y, x, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_isnans( eiy, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eic, eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"atan2", __DOMAIN__, eic, y, x, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for z, hi, lo, t, signy, signx . *******/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } hi = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !hi ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } lo = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !lo ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } signy = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !signy ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } signx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !signx ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for zero, one, two, small, big . *******/ zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } big = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !big ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ __mpu_sbrk( -(int)(12*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for tx, ty . ***************************/ tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ __mpu_sbrk( -(int)(13*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ty = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ty ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(14*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for k, m, kt, mk, prec1, mprec2 . ******/ k = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !k ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } m = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !m ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } kt = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !kt ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ __mpu_sbrk( -(int)(2*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } mk = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !mk ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ __mpu_sbrk( -(int)(3*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } prec1 = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !prec1 ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ __mpu_sbrk( -(int)(4*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } mprec2 = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !mprec2 ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ __mpu_sbrk( -(int)(5*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /********************************************* small: 1.0 + small^2 == 1.0; big : big = 1/(small^2); *********************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); ei_div( big, one, _get_epsilon_ptr( nb ), nb ); pi = _get_m_pi_ptr ( nb ); pi_2 = _get_m_pi_2_ptr ( nb ); pi_4 = _get_m_pi_4_ptr ( nb ); pi3_4 = _get_m_3pi_4_ptr( nb ); /* Service Constants */ athfhi = _get_m_athfhi_ptr( nb ); athflo = _get_m_athflo_ptr( nb ); at1fhi = _get_m_at1fhi_ptr( nb ); at1flo = _get_m_at1flo_ptr( nb ); c39_16 = _get_m_39_16_ptr ( nb ); c1_16 = _get_m_1_16_ptr ( nb ); /* ei_copy( y, eiy, nb ); */ /* ei_copy( x, eix, nb ); */ /************************************ Copy down the sign of EIY and EIX. ************************************/ ei_copysign( signy, one, y, nb ); ei_copysign( signx, one, x, nb ); /************************** If x is 1.0, goto begin. **************************/ if( ei_cmp( x, one, nb ) == 0 ) { ei_copysign( y, y, one, nb ); ei_copy( t, y, nb ); if( !ei_isinfin( t, nb ) ) goto begin; } /************* If y is 0.0 *************/ if( ei_cmp( y, zero, nb ) == 0 ) { if( ei_cmp( signx, one, nb ) == 0 ) ei_copy( eic, y, nb ); else ei_copysign( eic, pi, signy, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************* If x is 0.0 *************/ if( ei_cmp( x, zero, nb ) == 0 ) { ei_copysign( eic, pi_2, signy, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /******************** When x is Infinity ********************/ if( ei_isinfin( x, nb ) ) { if( ei_isinfin( y, nb ) ) { if( ei_cmp( signx, one, nb ) == 0 ) ei_copysign( eic, pi_4, signy, nb ); else ei_copysign( eic, pi3_4, signy, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { if( ei_cmp( signx, one, nb ) == 0 ) ei_copysign( eic, zero, signy, nb ); else ei_copysign( eic, pi, signy, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x == Infinity ) */ /******************** When y is Infinity ********************/ if( ei_isinfin( y, nb ) ) { ei_copysign( eic, pi_2, signy, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( y == Infinity ) */ /************* Compute y/x *************/ /* prec1 = NSBITS(nb); mprec2 = -(3/2)*NSBITS(nb); */ prec = NSBITS(nb); ei_cpye_unpack( prec1, (EMUSHORT *)&prec, ne, 1 ); prec /= 2; ei_cpye_unpack( mprec2, (EMUSHORT *)&prec, ne, 1 ); ei_adde( mprec2, mprec2, prec1, ne ); ei_nege( mprec2, mprec2, ne ); ei_copysign( y, y, one, nb ); ei_copysign( x, x, one, nb ); ei_logb( k, y, ne, nb ); /* k = logb( y ); */ ei_logb( kt, x, ne, nb ); /* kt = logb( x ); */ ei_sube( m, k, kt, ne ); /* m = (k=logb(y)) - logb(x); */ ei_nege( mk, k, ne ); /* mk = -k; */ if( ei_cmpe( m, prec1, ne ) > 0 ) { ei_add( t, big, big, nb ); /* t = big+big; */ } else if( ei_cmpe( m, mprec2, ne ) < 0 ) { ei_div( t, y, x, nb ); /* t = y/x; */ } else { ei_div( t, y, x, nb ); /* t = y/x; */ ei_ldexp( y, mk, y, ne, nb ); /* y = ldexp( y, -k ); */ ei_ldexp( x, mk, x, ne, nb ); /* x = ldexp( x, -k ); */ } /************************************************************* BEGIN ARGUMENT REDUCTION: *************************************************************/ begin: if( ei_cmp( t, c39_16, nb ) < 0 ) /* if( t < 39/16 ) */ { /* truncate 4(t+1/16) to integer for branching */ ei_add( tx, t, c1_16, nb ); ei_mul( tx, tx, two, nb ); ei_mul( tx, tx, two, nb ); ei_rtol_frac( (EMUSHORT *)&j, /* j = 4*(t+0.0625); */ (EMUSHORT *)0, tx, 1, nb ); switch( j ) { /* t is in [0, 7/16] */ case 0: case 1: if( ei_cmp( t, small, nb ) < 0 ) { if( ei_cmp( signx, zero, nb ) > 0 ) { ei_copysign( eic, t, signy, nb ); } else { ei_sub( t, pi, t, nb ); ei_copysign( eic, t, signy, nb ); } /* raise the inexact flag */ /* big+small */ ei_copy( y, eiy, nb ); /* for &eiy == &eic */ ei_copy( x, eix, nb ); /* for &eix == &eic */ /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( (EMUSHORT *)0, /* Non change */ (__mpu_char8_t *)"atan2", __INEXACT__, eic, y, x, nb ); __STINX; /* produsing Inexact Flag */ /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_copy( hi, zero, nb ); ei_copy( lo, zero, nb ); break; /* t is in [7/16, 11/16] */ case 2: ei_copy( hi, athfhi, nb ); ei_copy( lo, athflo, nb ); ei_add( z, x, x, nb ); /* z = x+x; */ ei_add( ty, y, y, nb ); ei_sub( ty, ty, x, nb ); ei_add( tx, z, y, nb ); ei_div( t, ty, tx, nb ); /* t = ((y+y)-x)/(z+y); */ break; /* t is in [11/16, 19/16] */ case 3: case 4: ei_copy( hi, pi_4, nb ); ei_copy( lo, zero, nb ); ei_sub( ty, y, x, nb ); ei_add( tx, x, y, nb ); ei_div( t, ty, tx, nb ); /* t = (y-x)/(x+y); */ break; /* t is in [19/16, 39/16] */ default: ei_copy( hi, at1fhi, nb ); ei_copy( lo, at1flo, nb ); ei_sub( z, y, x, nb ); /* z = y-x; */ ei_add( ty, y, y, nb ); ei_add( y, ty, y, nb ); /* y = y+y+y; */ ei_add( t, x, x, nb ); /* t = x+x; */ ei_add( ty, z, z, nb ); ei_sub( ty, ty, x, nb ); ei_add( tx, t, y, nb ); ei_div( t, ty, tx, nb ); /* t = ((z+z)-x)/(t+y); */ break; } /* End of switch( j ) */ } /* End if( t < 39/16 ) */ else { ei_copy( hi, pi_2, nb ); ei_copy( lo, zero, nb ); /* t is in [39/16, big] */ if( ei_cmp( t, big, nb ) <= 0 ) { ei_copy( tx, x, nb ); ei_neg( tx, nb ); ei_div( t, tx, y, nb ); /* t = -x/y; */ } /* t is in [big, Infinity] */ else { /* raise the inexact flag */ /* big+small */ /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( (EMUSHORT *)0, /* Non change */ (__mpu_char8_t *)"atan2", __INEXACT__, (EMUSHORT *)0, eiy, eix, nb ); __STINX; /* produsing Inexact Flag */ ei_copy( t, zero, nb ); } } /************************************************************* :END OF ARGUMENT REDUCTION. *************************************************************/ /************************************************************* COMPUTE ATAN(t) FOR t in [-7/16, 7/16]. *************************************************************/ ei_mul( z, t, t, nb ); ei_atan2__A( tx, z, nb ); ei_mul( z, t, tx, nb ); ei_sub( z, lo, z, nb ); ei_add( z, z, t, nb ); ei_add( z, z, hi, nb ); if( ei_cmp( signx, zero, nb ) > 0 ) { ei_copysign( eic, z, signy, nb ); } else { ei_sub( z, pi, z, nb ); ei_copysign( eic, z, signy, nb ); } /* FREE x *****************/ /* FREE y *****************/ /* FREE z *****************/ /* FREE hi ****************/ /* FREE lo ****************/ /* FREE t *****************/ /* FREE signy *************/ /* FREE signx *************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE two ***************/ /* FREE small *************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(15*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE k *****************/ /* FREE m *****************/ /* FREE kt ****************/ /* FREE mk ****************/ /* FREE prec1 *************/ /* FREE mprec2 ************/ __mpu_sbrk( -(int)(6*ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End of ei_atan2() */ void ei_atan( EMUSHORT *eiy, EMUSHORT *eix, int nb ) { EMUSHORT *one = NULL; int np; np = internal_np( nb ); /*** Allocate memory for one . ******************************/ one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ return; } /************************************************************/ _gen_one( one, nb ); ei_atan2( eiy, eix, one, nb ); /* FREE one ***************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ } /* End of ei_atan() */ void ei_sinh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_sinh() Работает с internal e-type data struct. Concepts : SINH(EIX) RETURN THE HYPERBOLIC SINE OF EIX. METHOD : 1. Reduce x to non-negative by sinh(-x) = - sinh(x). 2. expm1(x) + expm1(x)/(expm1(x)+1) 0 <= x <= lnovfl: sinh(x) = --------------------------------; 2 lnovfl <=x<= lnovfl+ln2: sinh(x) = expm1(x)/2 (avoid overflow); lnovfl+ln2 < x < INF: overflow to INF; ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: FLOATP.H]. NOTE : SPECIAL CASES : sinh(-InD) = -InD с выставлением DOMAIN flag; sinh(+NaN) = +NaN с выставлением DOMAIN flag; sinh(-NaN) = -NaN с выставлением DOMAIN flag; sinh(-inf) = -Inf [норма]; sinh(+inf) = +inf [норма]; sinh(+big#)= +Inf с выставлением OVERFLOW flag; sinh(-big#)= -Inf с выставлением OVERFLOW flag; sinh(0) = 0 [норма]; Use Global Variable: Use Functions : ei_expm1( eiy, eix, nb ); | this file ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *sign = NULL, *ln2p = NULL, *one = NULL, *half = NULL, *tr = NULL; /* temp */ EMUSHORT *max = NULL; /* for MAX_BIN_Exponent */ EMUSHORT *ln2, *lnovfl; int np, ne; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"sinh" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* SINH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sinh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* SINH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sinh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* SINH(+Infinity) must by +Infinity; SINH(-Infinity) must by -Infinity; */ if( ei_isinfin( eix, nb ) ) { /* Normal Exit */ ei_copy( eiy, eix, nb ); /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for t, sign, ln2p, one, half, tr . *****/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } sign = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !sign ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ln2p = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ln2p ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ /* FREE half **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for max . ******************************/ max = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !max ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ /* FREE half **************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ ln2 = _get_m_ln2_ptr( nb ); ei_cpye( max, _get_max_2_exp_ptr( nb ), ne, ne ); _gen_half( half, nb ); _gen_one ( one, nb ); /* Service Constant: */ lnovfl = _get_m_ln_huge_ptr( nb ); ei_add( ln2p, lnovfl, ln2, nb ); /* ei_copy( x, eix, nb ); */ ei_copysign( sign, one, x, nb ); ei_copysign( x, x, one, nb ); if( ei_cmp( x, lnovfl, nb ) < 0 ) { ei_expm1( t, x, nb ); /* return( copysign((t+t/(one+t))*half, sign) ); */ ei_add( tr, one, t, nb ); ei_div( tr, t, tr, nb ); ei_add( tr, t, tr, nb ); ei_mul( tr, tr, half, nb ); ei_copysign( eiy, tr, sign, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ /* FREE half **************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { if( ei_cmp( x, ln2p, nb ) <= 0 ) /* if( x <= lnovfl+ln2 ) */ { /* Substract x by ln(2^(max+1)) and return 2^max*exp(x) to avoid unnecessary overflow. ******************************************************/ ei_sub( tr, x, lnovfl, nb ); ei_expm1( tr, tr, nb ); ei_add( tr, one, tr, nb ); ei_ldexp( tr, max, tr, ne, nb ); ei_copysign( eiy, tr, sign, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ /* FREE half **************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* SINH(+-big#) overflow to +-Infinity */ /* return: +-Infinity */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_expm1( tr, x, nb ); ei_mul( eiy, tr, sign, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную __error_no следующим образом. __error_no = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"sinh", __OVERFLOW__, eiy, x, (EMUSHORT *)0, nb ); __STOVF; /* +-big# - produsing Overflow Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE ln2p **************/ /* FREE one ***************/ /* FREE half **************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if/else( x <= lnovfl+ln2 ) */ } /* Enf if/else( x < lnovfl ) */ } /* End of ei_sinh() */ void ei_cosh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_cosh() Работает с internal e-type data struct. Concepts : COSH(EIX) RETURN THE HYPERBOLIC COSINE OF EIX. METHOD : 1. Replace x by |x|. 2. 0 <= x <= ln2/2 : [ exp(x) - 1 ]^2 cosh(x) = 1 + -------------------; 2*exp(x) ln2/2 <= x <= thovfl : exp(x) + 1/exp(x) cosh(x) = -------------------; 2 thovfl <= x <= lnovfl: cosh(x) = exp(x)/2; lnovfl<=x<=lnovfl+ln2: cosh(x) = exp(x)/2 (avoid overflow); ln2+lnovfl < x < INF : overflow to INF; ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : thovfl = log1p(2.0/EPSILON-2.0) / 2.0; [small == EPSILON]. SPECIAL CASES : cosh(-InD) = -InD с выставлением DOMAIN flag; cosh(+NaN) = +NaN с выставлением DOMAIN flag; cosh(-NaN) = -NaN с выставлением DOMAIN flag; cosh(-inf) = +Inf [норма]; cosh(+inf) = +Inf [норма]; cosh(+big#)= +Inf с выставлением OVERFLOW flag; cosh(-big#)= +Inf с выставлением OVERFLOW flag; cosh(0) = 1.0 [норма]; Use Global Variable: Use Functions : ei_exp( eiy, eix, nb ); | this file __ei_exp__E( eiy, eix, eic, nb ); | this file ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *ln2p = NULL, *ln2_2 = NULL, *zero = NULL, *half = NULL, *one = NULL, *two = NULL, *tr = NULL; /* temp */ EMUSHORT *max = NULL; /* for MAX_BIN_Exponent */ EMUSHORT *ln2, *lnovfl, *thovfl, *small; int np, ne; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"cosh" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* COSH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"cosh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* COSH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"cosh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* COSH(+Infinity) must by +Infinity; COSH(-Infinity) must by -Infinity; */ if( ei_isinfin( eix, nb ) ) { /* Normal Exit */ ei_copy( eiy, eix, nb ); /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ne = internal_ne( nb ) + 1; /*** Allocate memory for t, ln2p, ln2_2, tr . ***************/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ln2p = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ln2p ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ln2_2 = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ln2_2 ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for zero, half, one, two . *************/ zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for max . ******************************/ max = (EMUSHORT *)__mpu_sbrk( (int)(ne*SIZE_OF_EMUSHORT) ); if( !max ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ ln2 = _get_m_ln2_ptr( nb ); small = _get_epsilon_ptr( nb ); ei_cpye( max, _get_max_2_exp_ptr( nb ), ne, ne ); _gen_zero( zero, nb ); _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* Service Constant: */ lnovfl = _get_m_ln_huge_ptr( nb ); thovfl = _get_m_thovfl_ptr( nb ); ei_add( ln2p, lnovfl, ln2, nb ); ei_div( ln2_2, ln2, two, nb ); ei_copysign( x, x, one, nb ); if( ei_cmp( x, thovfl, nb ) <= 0 ) { if( ei_cmp( x, ln2_2, nb ) < 0 ) { if( ei_cmp( x, small, nb ) < 0 ) { ei_add( eiy, one, x, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { __ei_exp__E( tr, x, zero, nb ); ei_add( t, x, tr, nb ); /* t = x+__ei_exp__E(x,0.0); */ ei_add( x, t, t, nb ); /* x = t+t; */ ei_add( tr, two, x, nb ); ei_div( tr, t, tr, nb ); ei_mul( tr, t, tr, nb ); ei_add( eiy, one, tr, nb ); /* return(one+t*t/(2.0+x)); */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } else { /******************************* for x lies in [ln2/2, thovfl] *******************************/ ei_exp( t, x, nb ); /* return( (t+one/t)*half ); */ ei_div( tr, one, t, nb ); ei_add( tr, t, tr, nb ); ei_mul( eiy, tr, half, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x <= thovfl ) */ /* if( lnovfl <= x && x <= lnovfl+ln2 ) */ if( (ei_cmp( lnovfl, x, nb ) <= 0) && (ei_cmp( x, ln2p, nb ) <= 0) ) { /* Substract x by ln(2^(max+1)) and return 2^max*exp(x) to avoid unnecessary overflow. ******************************************************/ ei_sub( tr, x, lnovfl, nb ); ei_exp( tr, tr, nb ); ei_ldexp( eiy, max, tr, ne, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { ei_copy( t, eix, nb ); /* for &eiy == &eix */ ei_exp( tr, x, nb ); /* new version: */ ei_mul( eiy, tr, half, nb ); if( ei_cmp( x, ln2p, nb ) > 0 ) { /* for large (x != +-Infinity), cosh(x) = exp(x)/2; */ /* COSH(+-big#) overflow to +-Infinity */ /* return: +-Infinity ( after if(){} body ) */ /* previus version: ei_mul( eiy, tr, half, nb ); for _mtherr() */ /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, /* previus version: No change */ (__mpu_char8_t *)"cosh", __OVERFLOW__, eiy, t, (EMUSHORT *)0, nb ); __STOVF; /* +-big# - produsing Overflow Flag */ } /* previus version: ei_mul( eiy, tr, half, nb ); */ /* FREE x *****************/ /* FREE t *****************/ /* FREE ln2p **************/ /* FREE ln2_2 *************/ /* FREE tr ****************/ /* FREE zero **************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ /* FREE max ***************/ __mpu_sbrk( -(int)(ne*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if/else( lnovfl <= x && x <= lnovfl+ln2 ) */ } /* End of ei_cosh() */ void ei_tanh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_tanh() Работает с internal e-type data struct. Concepts : TANH(EIX) RETURN THE HYPERBOLIC TANGENT OF EIX. METHOD : 1. Reduce x to non-negative by tanh(-x) = -tanh(x). 2. 0 MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"tanh" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* TANH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"tanh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* TANH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"tanh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, sign, small, tr . ***************/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } sign = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !sign ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for half, one, two . *******************/ half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); _gen_half( half, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); /* Service Constant: */ thovfl = _get_m_thovfl_ptr( nb ); ei_copysign( sign, one, x, nb ); ei_copysign( x, x, one, nb ); /* TANH(+Infinity) must by +1.0; TANH(-Infinity) must by -1.0; */ if( ei_isinfin( eix, nb ) ) { /* Normal Exit */ ei_copy( eiy, sign, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_cmp( x, thovfl, nb ) < 0 ) { if( ei_cmp( x, one, nb ) > 0 ) { /* return( copysign( one-two/(expm1(x+x)+two), sign ) ); */ ei_add( tr, x, x, nb ); ei_expm1( tr, tr, nb ); ei_add( tr, tr, two, nb ); ei_div( tr, two, tr, nb ); ei_sub( tr, one, tr, nb ); ei_copysign( eiy, tr, sign, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { if( ei_cmp( x, small, nb ) > 0 ) { /* t = -expm1(-(x+x)); */ ei_add( tr, x, x, nb ); ei_neg( tr, nb ); ei_expm1( tr, tr, nb ); ei_neg( tr, nb ); ei_copy( t, tr, nb ); /* return( copysign( t/(two-t), sign ) ); */ ei_sub( tr, two, t, nb ); ei_div( tr, t, tr, nb ); ei_copysign( eiy, tr, sign, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { ei_copy( t, eix, nb ); /* for &eiy == &eix */ /* raise the inexact flag for non-zero x*/ /* return: +-x */ ei_copysign( eiy, x, sign, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( (EMUSHORT *)0, /* Non change */ (__mpu_char8_t *)"tanh", __INEXACT__, eiy, t, (EMUSHORT *)0, nb ); __STINX; /* produsing Inexact Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if/else( x > small ) */ } /* End if/else( x > one ) */ } else { ei_copy( t, eix, nb ); /* for &eiy == &eix */ /* TANH(+-big#) round to +-1.0 */ /* raise the inexact flag */ /* return: +-1.0 */ ei_copy( eiy, sign, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( (EMUSHORT *)0, /* Non change */ (__mpu_char8_t *)"tanh", __INEXACT__, eiy, t, (EMUSHORT *)0, nb ); __STINX; /* produsing Inexact Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE sign **************/ /* FREE small *************/ /* FREE tr ****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE two ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* Enf if/else( x < thovfl ) */ } /* End of ei_tanh() */ void ei_asinh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_asinh() Работает с internal e-type data struct. Concepts : ASINH(X) RETURN THE INVERSE HYPERBOLIC SINE OF EIX. METHOD : Based on asinh(x) = sign(x)*log[ |x|+sqrt(x*x+1) ] we have asinh(x) = x if(1+x*x=1), = sign(x)*(log1p(x)+ln2) if(sqrt(1+x*x)=x), else = sign(x)* log1p(|x|+ |x|/ (1/|x|+sqrt(1+(1/|x|)^2)) ). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : SPECIAL CASES : asinh(-InD) = -InD с выставлением DOMAIN flag; asinh(+NaN) = +NaN с выставлением DOMAIN flag; asinh(-NaN) = -NaN с выставлением DOMAIN flag; Use Global Variable: Use Functions : ei_log1p( eiy, eix, nb ); | this file ei_sqrt( eiy, eix, nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *s = NULL, *small = NULL, *big = NULL, *one = NULL, *ten = NULL, *tr = NULL; /* temp */ EMUSHORT *ln2hi, *ln2lo; int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"asinh" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* ASINH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"asinh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ASINH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"asinh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, s, small, big, one, ten, tr . ***/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } s = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } small = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !small ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } big = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !big ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ten = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ten ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ /* FREE one ***************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_one ( one, nb ); _gen_ten ( ten, nb ); /********************************************* small: 1.0 + small^2 == 1.0; big : big = 1/(small^2); & (big+1) == big; *********************************************/ ei_sqrt( small, _get_epsilon_ptr( nb ), nb ); ei_div( big, one, _get_epsilon_ptr( nb ), nb ); ei_mul( big, big, ten, nb ); /* Service Constant: */ ln2hi = _get_m_ln2hi_ptr( nb ); ln2lo = _get_m_ln2lo_ptr( nb ); ei_copysign( t, x, one, nb ); if( ei_cmp( t, small, nb ) > 0 ) { if( ei_cmp( t, big, nb ) < 0 ) { ei_div( s, one, t, nb ); /* s = one/t; */ /* return( copysign( log1p( t+t/(s+sqrt(one+s*s)) ), x ) ); */ ei_mul( tr, s, s, nb ); ei_add( tr, one, tr, nb ); ei_sqrt( tr, tr, nb ); ei_add( tr, s, tr, nb ); ei_div( tr, t, tr, nb ); ei_add( tr, t, tr, nb ); ei_log1p( tr, tr, nb ); ei_copysign( eiy, tr, x, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* ( |x| > big ) */ ei_log1p( tr, t, nb ); ei_add( s, tr, ln2lo, nb ); /* s = log1p(t)+ln2lo; */ ei_add( tr, s, ln2hi, nb ); /* return( copysign(s+ln2hi,x) ); */ ei_copysign( eiy, tr, x, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } else { /* ( |x| < small ) */ ei_copy( eiy, x, nb ); /* return( x ); */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE small *************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End of ei_asinh() */ void ei_acosh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_acosh() Работает с internal e-type data struct. Concepts : ASINH(X) RETURN THE INVERSE HYPERBOLIC COSINE OF EIX. METHOD : Based on acosh(x) = log[ x + sqrt(x*x-1) ] we have acosh(x) = log1p(x)+ln2 if(x>big), else = log1p( sqrt(x-1) * ( sqrt(x-1) + sqrt(x+1) ) ). These formulae avoid the over/underflow complication. ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : SPECIAL CASES : acosh(-InD) = -InD с выставлением DOMAIN flag; acosh(+NaN) = +NaN с выставлением DOMAIN flag; acosh(-NaN) = -NaN с выставлением DOMAIN flag; acosh( x<1 ) = -Ind с выставлением DOMAIN flag; Use Global Variable: Use Functions : ei_log1p( eiy, eix, nb ); | this file ei_sqrt( eiy, eix, nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *big = NULL, *one = NULL, *ten = NULL, *tr = NULL; /* temp */ EMUSHORT *ln2hi, *ln2lo; int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"acosh" ); /* ei_ind( eiy, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* ACOSH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acosh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ACOSH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acosh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, big, one, ten, tr . *************/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } big = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !big ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ten = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ten ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ /* FREE one ***************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_one ( one, nb ); /* ACOSH( x < 1 ) must by InD */ if( ei_cmp( eix, one, nb ) < 0 ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acosh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* (x < 1) - produsing Domain Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } _gen_ten ( ten, nb ); /********************************************* big : big = 1/(small^2); & (big+1) == big; *********************************************/ ei_div( big, one, _get_epsilon_ptr( nb ), nb ); ei_mul( big, big, ten, nb ); /* Service Constant: */ ln2hi = _get_m_ln2hi_ptr( nb ); ln2lo = _get_m_ln2lo_ptr( nb ); /* ei_copy( x, eix, nb ); */ if( ei_cmp( x, big, nb ) > 0 ) { /* return( log1p(x) + ln(2) ) if x is large. */ ei_log1p( t, x, nb ); ei_add( t, t, ln2lo, nb ); /* t = log1p(x)+ln2lo; */ ei_add( eiy, t, ln2hi, nb ); /* return( t+ln2hi ); */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_sub( tr, x, one, nb ); ei_sqrt( t, tr, nb ); /* t = sqrt(x-1.0); */ ei_add( tr, x, one, nb ); ei_sqrt( tr, tr, nb ); ei_add( tr, t, tr, nb ); ei_mul( tr, t, tr, nb ); ei_log1p( eiy, tr, nb ); /* return( log1p( t*(t+sqrt(x+1.0)) ) ); */ /* FREE x *****************/ /* FREE t *****************/ /* FREE big ***************/ /* FREE one ***************/ /* FREE ten ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End of ei_acosh() */ void ei_atanh( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_atanh() Работает с internal e-type data struct. Concepts : ATANH(X) RETURN THE INVERSE HYPERBOLIC TANGENT OF EIX. METHOD : Return 1 2x x atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------); 2 1 - x 1 - x ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : Для r в интервале [0, thovfl], проделаем следующее: t = tanh ( r ); t = atanh( t ). При этом для r в интервале [1.0, thovfl] с увеличением r растет и величина разности (t-r) чего в принципе не должно быть. СЛЕДОВАТЕЛЬНО ПОКА МЫ МОЖЕМ РЕКОМЕНДОВАТЬ проводить такие опыты лишь на промежутке [-1, 1]. SPECIAL CASES : atanh(-InD) = -InD с выставлением DOMAIN flag; atanh(+NaN) = +NaN с выставлением DOMAIN flag; atanh(-NaN) = -NaN с выставлением DOMAIN flag; atanh(|x|>1) = -Ind с выставлением DOMAIN flag; atanh( +-1 ) = +-Inf с выставлением OVERFLOW flag; Use Global Variable: Use Functions : ei_log1p( eiy, eix, nb ); | this file ei_copy( eiy, eix, nb ); | mpu-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, *z = NULL, *half = NULL, *one = NULL, *tr = NULL; /* temp */ int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"atanh" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* ATANH(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"atanh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ATANH(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"atanh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for z, half, one, tr . *****************/ z = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !z ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE half **************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); /* ei_copy( x, eix, nb ); */ ei_copysign( z, half, x, nb ); ei_copysign( x, x, one, nb ); /* ATANH( |x| > 1 ) must by InD */ if( ei_cmp( x, one, nb ) > 0 ) { /* "argument domain error" */ /* return: InD */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"atanh", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* (x > 1) - produsing Domain Flag */ /* FREE x *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ATANH( |x| == 1 ) must by sign(x)*Infinity */ if( ei_cmp( x, one, nb ) == 0 ) { /* ATANH(|x| == 1) overflow to +-Infinity */ /* "overflow range error" */ /* return: sign(x)*Infinity */ ei_copy( x, eix, nb ); /* for &eiy == &eix */ ei_infin( tr, (unsigned)0, nb ); ei_copysign( eiy, tr, z, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"atanh", __OVERFLOW__, eiy, x, (EMUSHORT *)0, nb ); __STOVF; /* (|x| == 1) - produsing Overflow Flag */ /* FREE x *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ei_sub( tr, one, x, nb ); ei_div( x, x, tr, nb ); /* x = x/(1.0-x); */ ei_add( tr, x, x, nb ); ei_log1p( tr, tr, nb ); ei_mul( eiy, z, tr, nb ); /* return( z*log1p(x+x) ); */ /* FREE x *****************/ /* FREE z *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End of ei_atanh() */ void ei_asin( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_asin() Работает с internal e-type data struct. Concepts : ASIN(X) RETURN THE INVERSE SINE OF EIX. METHOD : Return asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is computed as follows if( x < 0.5) 1-x*x; if( x >= 0.5) 2*(1-|x|)-(1-|x|)*(1-|x|). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : SPECIAL CASES : asin(-InD) = -InD с выставлением DOMAIN flag; asin(+NaN) = +NaN с выставлением DOMAIN flag; asin(-NaN) = -NaN с выставлением DOMAIN flag; asin(|x|>1) = -Ind с выставлением DOMAIN flag; Use Global Variable: Use Functions : ei_atan2( eiy, eix, nb ); | this file ei_sqrt( eiy, eix, nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *s = NULL, *half = NULL, *one = NULL, *tr = NULL; /* temp */ int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"asin" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* ASIN(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"asin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ASIN(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"asin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, s, half, one, tr . **************/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } s = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } half = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !half ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE half **************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_half( half, nb ); _gen_one ( one, nb ); ei_copysign( s, x, one, nb ); /* ATANH( |x| > 1 ) must by InD */ if( ei_cmp( s, one, nb ) > 0 ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"asin", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* (x > 1) - produsing Domain Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_cmp( s, half, nb ) <= 0 ) { /* return( atan2( x, sqrt( one - x*x ) ) ); */ ei_mul( tr, x, x, nb ); ei_sub( tr, one, tr, nb ); ei_sqrt( tr, tr, nb ); ei_atan2( eiy, x, tr, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { ei_sub( t, one, s, nb ); /* t = one - s; */ ei_add( s, t, t, nb ); /* s = t + t; */ /* return( atan2( x, sqrt( s - t*t ) ) ); */ ei_mul( tr, t, t, nb ); ei_sub( tr, s, tr, nb ); ei_sqrt( tr, tr, nb ); ei_atan2( eiy, x, tr, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE half **************/ /* FREE one ***************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End of ei_asin() */ void ei_acos( EMUSHORT *eiy, EMUSHORT *eix, int nb ) /*************************************************************** Description : ei_acos() Работает с internal e-type data struct. Concepts : ACOS(X) RETURN THE INVERSE COSINE OF EIX. METHOD : Return ________ / 1 - x acos(x) = 2*atan2( / --------, 1 ). \/ 1 + x ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : SPECIAL CASES : acos(-InD) = -InD с выставлением DOMAIN flag; acos(+NaN) = +NaN с выставлением DOMAIN flag; acos(-NaN) = -NaN с выставлением DOMAIN flag; acos(|x|>1) = -Ind с выставлением DOMAIN flag; Use Global Variable: Use Functions : ei_atan2( eiy, eix, nb ); | this file ei_sqrt( eiy, eix, nb ); | mpu-real.c ei_copy( eiy, eix, nb ); | mpu-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, *t = NULL, *s = NULL, *zero = NULL, *one = NULL, *negone = NULL, *tx = NULL, /* temp */ *tr = NULL; /* temp */ int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"acos" ); /* ei_ind( eiy, nb ); */ /* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for x . ********************************/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ return; } /************************************************************/ ei_copy( x, eix, nb ); /* temp for _mtherr() */ /*************************** Test for EIX. ***************************/ /* ACOS(InD) must by InD */ if( ei_isind( eix, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* ACOS(NaN) must by NaN */ if( ei_isnans( eix, nb ) ) { /* "argument domain error" */ /* return: NaN */ ei_copy( eiy, eix, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* NaN - produsing Domain Flag */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, s, zero, one, negone, tx, tr . **/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } s = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } negone = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !negone ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tr = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tr ) { /* fatal error */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); _gen_one ( one, nb ); _gen_one ( negone, nb ); ei_neg ( negone, nb ); ei_copysign( s, x, one, nb ); /* ATANH( |x| > 1 ) must by InD */ if( ei_cmp( s, one, nb ) > 0 ) { /* "argument domain error" */ /* return: InD */ ei_ind( eiy, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eiy, (__mpu_char8_t *)"acos", __DOMAIN__, eiy, x, (EMUSHORT *)0, nb ); __STDOM; /* (x > 1) - produsing Domain Flag */ /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE tx ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_cmp( x, negone, nb ) != 0 ) { /* t = atan2( sqrt( (one-x)/(one+x) ), one ); */ ei_sub( tx, one, x, nb ); ei_add( tr, one, x, nb ); ei_div( tx, tx, tr, nb ); ei_sqrt( tx, tx, nb ); ei_atan2( t, tx, one, nb ); } else { ei_atan2( t, one, zero, nb ); /* t = PI/2; */ } ei_add( eiy, t, t, nb ); /* FREE x *****************/ /* FREE t *****************/ /* FREE s *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE tx ****************/ /* FREE tr ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End of ei_acos() */ static void __ei_pow__P( EMUSHORT *eic, EMUSHORT *eix, EMUSHORT *eiy, int nb ) { EMUSHORT *s_a = NULL, *s_b = NULL, *t_a = NULL, *t_b = NULL; Double s, t; EMUSHORT *x = NULL, *y = NULL, *zero = NULL, *one = NULL, *big = NULL, *tx = NULL, /* temp */ *ty = NULL; /* temp */ int np; 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 */ /* ei_ind( eic, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; np = internal_np( nb ); /*** Allocate memory for s_a, s_b, t_a, t_b . ***************/ s_a = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s_a ) { /* fatal error */ return; } s_b = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !s_b ) { /* fatal error */ /* FREE s_a ***************/ __mpu_sbrk( -(int)(np*SIZE_OF_EMUSHORT) ); /**************************/ return; } t_a = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t_a ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } t_b = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t_b ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for x, y, zero, one, big, tx, ty . *****/ x = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !x ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } y = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !y ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } big = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !big ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ty = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ty ) { /* fatal error */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(10*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ s.a = s_a; s.b = s_b; t.a = t_a; t.b = t_b; _gen_zero( zero, nb ); _gen_one ( one, nb ); /********************************************* big : big = 1/(small^2); *********************************************/ ei_div( big, one, _get_epsilon_ptr( nb ), nb ); ei_copy( x, eix, nb ); ei_copy( y, eiy, nb ); if( ei_cmp( x, zero, nb ) == 0 ) { if( ei_cmp( y, zero, nb ) > 0 ) { /* #(12) */ ei_copy( eic, zero, nb ); /* Not set the FLAGS */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* (14) DIV_BY_ZERO */ /* "function singularity" */ /* return: +Infinity */ ei_infin( eic, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __SING__, eic, x, y, nb ); __STSNG; /* +0^(-(y!=0,NaN)) - produsing Sing Flag */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x == zero ) */ if( ei_cmp( x, one, nb ) == 0 ) { ei_copy( eic, one, nb ); /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } if( ei_isinfin( x, nb ) ) { if( ei_cmp( y, zero, nb ) < 0 ) { /* #(18) */ ei_copy( eic, zero, nb ); /* Not set the FLAGS */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* #(17) */ /* "overflow range error" */ /* "Result too large" */ /* return: +Infinity */ ei_infin( eic, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __OVERFLOW__, eic, x, y, nb ); __STOVF; /* produsing Overflow Flag */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x == Infinity ) */ if( ei_cmp( y, big, nb ) >= 0 ) { if( ei_cmp( x, one, nb ) < 0 ) { /* "underflow range error" */ /* return: ZERO */ ei_copy( eic, zero, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __UNDERFLOW__, eic, x, y, nb ); __STUDF; /* produsing Underflow Flag */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* "overflow range error" */ /* "Result too large" */ /* return: +Infinity */ ei_infin( eic, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __OVERFLOW__, eic, x, y, nb ); __STOVF; /* produsing Overflow Flag */ /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( y >= big ) */ /********************************************************* Return exp( y*log( x ) ), using simulated extended precision for the log() and the multiply. *********************************************************/ __ei_log__D( s, x, nb ); /* s = log__D( x ); */ ei_copy( t.a, y, nb ); /* t.a = y; */ __ei_TRUNC( t.a, nb ); /* _TRUNC( t.a ); */ ei_sub( t.b, y, t.a, nb ); /* t.b = y - t.a; */ ei_mul( ty, s.b, y, nb ); ei_mul( tx, t.b, s.a, nb ); ei_add( t.b, ty, tx, nb ); /* t.b = s.b*y + t.b*s.a; */ ei_mul( t.a, t.a, s.a, nb ); /* t.a *= s.a; */ ei_add( s.a, t.a, t.b, nb ); /* s.a = t.a + t.b; */ ei_sub( tx, t.a, s.a, nb ); ei_add( s.b, tx, t.b, nb ); /* s.b = (t.a - s.a) + t.b; */ /********* return( __ei_exp__D( s.a, s.b ) ); **********/ __ei_exp__D( eic, s.a, s.b, nb ); /* FREE s_a ***************/ /* FREE s_b ***************/ /* FREE t_a ***************/ /* FREE t_b ***************/ /* FREE x *****************/ /* FREE y *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE big ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(11*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End of __ei_pow__P() */ void ei_pow( EMUSHORT *eic, EMUSHORT *eix, EMUSHORT *eiy, int nb ) /*************************************************************** Description : ei_pow() Работает с internal e-type data struct. Concepts : EI_POW( EIC, EIX, EIY) RETURN EIC = EIX^EIY. METHOD : 1. Compute and return log(x) in three pieces: log(x) = n*ln2 + hi + lo, where n is an integer. 2. Perform y*log(x) by simulating muti-precision arithmetic and return the answer in three pieces: y*log(x) = m*ln2 + hi + lo, where m is an integer. 3. Return x^y = exp(y*log(x)) = 2^m*(exp(hi+lo)). ACCURACY : In the absence of rounding error, the appriximation has absolute error less than EPSILON [see: mpu-floatp.h]. NOTE : SPECIAL CASES : ==================================================================== | | X | Y | RESULT | | ====================|====================|========================== #( 1) (anything) ^ 0 = 1 ; (anything) ^ 1 = itself; (anything) ^ NaN = NaN; --------------------|--------------------|-------------------------- #( 4) NaN ^ (anything except 0) = NaN; --------------------|--------------------|-------------------------- #( 5) +(anything > 1) ^ +INF = +INF; -(anything > 1) ^ +INF = NaN; +-(anything > 1) ^ -INF = +0; +-(anything < 1) ^ +INF = +0; +(anything < 1) ^ -INF = +INF; -(anything < 1) ^ -INF = NaN; +-1 ^ +-INF = -InD and signal DOMAIN; --------------------|--------------------|-------------------------- #(12) +0 ^ +(anything except 0, NaN) = +0; --------------------|--------------------|-------------------------- #(13) -0 ^ +(anything except 0, NaN, odd integer) = +0; --------------------|--------------------|-------------------------- #(14) +0 ^ -(anything except 0, NaN) = +INF and signal DIV-BY-ZERO; --------------------|--------------------|-------------------------- #(15) -0 ^ -(anything except 0, NaN, odd integer) = +INF with signal; --------------------|--------------------|-------------------------- #(16) -0 ^ (odd integer) = -( +0^(odd integer) ); --------------------|--------------------|-------------------------- #(17) +INF ^ +(anything except 0, NaN) = +INF; --------------------|--------------------|-------------------------- #(18) +INF ^ -(anything except 0, NaN) = +0; --------------------|--------------------|-------------------------- #(19) -INF ^ (odd integer) = -( +INF^(odd integer) ); -INF ^ (even integer) = ( +INF^(even integer) ); --------------------|--------------------|-------------------------- #(21) -INF ^ -(anything except NaN, integer) = NaN with signal; --------------------|--------------------|-------------------------- #(22) -(x=anything) ^ (k=integer) = (-1)^k * (x^k); -(anything except 0) ^ (non-integer) = NaN with signal; ==================================================================== Use Global Variable: Use Functions : ei_drem( eiy, eix, eid, nb ); | this file __ei_pow__P( eic, eix, eiy, nb ); | this file ei_copy( eiy, eix, nb ); | mpu-real.c . . . Parameters : EMUSHORT *eic; - указатель на internal e-type data struct. TARGET; EMUSHORT *eix; - указатель на internal e-type data struct. SOURCE; EMUSHORT *eiy; - указатель на internal e-type data struct. SOURCE; int nb; - количество бит в external e-type data struct. Return : [void] ***************************************************************/ { EMUSHORT *x = NULL, *y = NULL, *t = NULL, *zero = NULL, *one = NULL, *negone = NULL, *two = NULL, *tx = NULL, /* temp */ *ty = NULL; /* temp */ int np; if( nb < NBR_32 || nb > MPU_MATH_FN_LIMIT ) { /* error: Invalid size of operand(s) */ __real_invalid_size( (__mpu_char8_t *)"pow" ); /* ei_ind( eic, nb ); *//* Invalid NB */ return; } /******************************* Hight Precision for r32, r64. *******************************/ if( nb < NBR_128 ) nb = NBR_128; 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; } /************************************************************/ ei_copy( x, eix, nb ); ei_copy( y, eiy, nb ); if( ei_isind( eix, nb ) || ei_isind( eiy, nb ) ) { /* "argument domain error" */ /* return: InD */ ei_ind( eic, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __DOMAIN__, eic, x, y, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /*** Allocate memory for t, zero, one, negone, two . ********/ t = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !t ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ __mpu_sbrk( -(int)(2*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } zero = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !zero ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ __mpu_sbrk( -(int)(3*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } one = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !one ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ __mpu_sbrk( -(int)(4*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } negone = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !negone ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ __mpu_sbrk( -(int)(5*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } two = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !two ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ __mpu_sbrk( -(int)(6*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ /*** Allocate memory for tx, ty . ***************************/ tx = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !tx ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ __mpu_sbrk( -(int)(7*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } ty = (EMUSHORT *)__mpu_sbrk( (int)(np*SIZE_OF_EMUSHORT) ); if( !ty ) { /* fatal error */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ __mpu_sbrk( -(int)(8*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /************************************************************/ _gen_zero( zero, nb ); _gen_one ( one, nb ); _gen_two ( two, nb ); _gen_one ( negone, nb ); ei_neg ( negone, nb ); if( ei_cmp( y, zero, nb ) == 0 ) { /* #( 1) */ ei_copy( eic, one, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( (ei_cmp( y, one, nb ) == 0) || ei_isnans( x, nb ) ) { /* #( 2) || #( 4) */ ei_copy( eic, x, nb ); /* Not set the FLAGS */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_isnans( y, nb ) ) { /* #( 3) */ ei_copy( eic, y, nb ); /* Not set the FLAGS */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_isinfin( eiy, nb ) ) { ei_copy( t, x, nb ); if( ei_cmp( t, one, nb ) == 0 ) { /* #(11) */ /* +-1^+-Inf = InD; */ /* "argument domain error" */ /* return: InD */ ei_ind( eic, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __DOMAIN__, eic, x, y, nb ); __STDOM; /* InD - produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_cmp( t, one, nb ) > 0 ) { if( ei_cmp( y, zero, nb ) < 0 ) { /* #( 7) */ ei_copy( eic, zero, nb ); /* Not set the FLAGS */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_cmp( x, zero, nb ) < 0 ) { /* #( 6) */ /* raise inexact flag */ ei_nanmax( eic, (unsigned)0, nb ); /* Not set the FLAGS */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else /* ( x >= zero ) */ { /* #( 5) */ ei_copy( eic, y, nb ); /* Not set the FLAGS */ /* NOTE: y == +Infinity */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( |x| > one ) */ else /* ( |x| < one ) */ { if( ei_cmp( y, zero, nb ) > 0 ) { /* #( 8) */ ei_copy( eic, zero, nb ); /* Not set the FLAGS */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_cmp( x, zero, nb ) < 0 ) { /* #(10) */ /* raise inexact flag */ ei_nanmax( eic, (unsigned)0, nb ); /* Not set the FLAGS */ /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, /* for No change set (EMUSHORT *)0 */ (__mpu_char8_t *)"pow", __INEXACT__, eic, x, y, nb ); __STINX; /* produsing Inexact Flag */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else /* ( x >= zero ) */ { /* #( 9) */ ei_copy( eic, y, nb ); /* Not set the FLAGS */ /* NOTE: y == +Infinity */ ei_neg( eic, nb ); /* return( -y ); */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End ( |x| < one ) */ } /* End if( eiy == +-Infinity ) */ else /* ( eiy != +-Infinity ) */ { if( ei_cmp( y, two, nb ) == 0 ) { ei_mul( eic, x, x, nb ); /* return( x*x ); */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else if( ei_cmp( y, negone, nb ) == 0 ) { ei_div( eic, one, x, nb ); /* return( 1.0/x ); */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( y == -1.0 ) */ else /* ( y != -1.0 ) */ { ei_copysign( tx, one, x, nb ); if( ei_cmp( tx, one, nb ) == 0 ) { /* x > 0, my be x == +0 */ __ei_pow__P( eic, x, y, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( x >= +0 ) */ else /* ( x < +0 ) */ { ei_drem( t, y, two, nb ); if( ei_cmp( t, zero, nb ) == 0 ) { /* y - четное целое (y is an even integer) */ ei_neg( x, nb ); __ei_pow__P( eic, x, y, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( y == even integer ) */ else /* ( y != even integer ) */ { ei_copysign( ty, t, one, nb ); if( ei_cmp( ty, one, nb ) == 0 ) { /* y - нечетное целое (y is an odd integer) */ ei_neg( x, nb ); __ei_pow__P( eic, x, y, nb ); ei_neg( eic, nb ); /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End if( y == odd integer ) */ else /* ( y != odd integer ) */ { /* y is not an integer */ if( ei_cmp( x, zero, nb ) == 0 ) { /* x is -0. NOTE: for ei_cmp() -0.0 equals +0.0. ***************************************/ if( ei_cmp( y, zero, nb ) > 0 ) { /* #(13) */ ei_neg( x, nb ); /* return +0 */ ei_copy( eic, x, nb ); /* return( -x ); */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } else { /* #(15) */ /////////////////////////////////////// //|unix|//ei_neg( x, nb ); // //|unix|///* return( 1/(-x) ); */ // //|unix|//ei_div( eic, one, x, nb ); // /////////////////////////////////////// /* "function singularity" */ /* -0^(-(y!=0,NaN,odd_int)) */ /* return: +Infinity */ ei_infin( eic, (unsigned)0, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __SING__, eic, x, y, nb ); __STSNG; /* produsing Sing Flag */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } } /* End if( x == -0 ) */ else /* ( x != -0 ) */ { /* #(23) */ /* "argument domain error" */ /* return: InD */ ei_ind( eic, nb ); /*************************************************** NOTE: Функция _mtherr() является переходником между внутренним форматом и внешней, переопределяемой пользователем, функцией _math_error(). Кроме основных действий она выставляет системную переменную errno следующим образом. errno = __mpu_math_errnotab[type]; ***************************************************/ _mtherr( eic, (__mpu_char8_t *)"pow", __DOMAIN__, eic, x, y, nb ); __STDOM; /* produsing Domain Flag */ /* FREE x *****************/ /* FREE y *****************/ /* FREE t *****************/ /* FREE zero **************/ /* FREE one ***************/ /* FREE negone ************/ /* FREE two ***************/ /* FREE tx ****************/ /* FREE ty ****************/ __mpu_sbrk( -(int)(9*np*SIZE_OF_EMUSHORT) ); /**************************/ return; } /* End ( x != -0 ) */ } /* End ( y != odd integer ) */ } /* End ( y != even integer ) */ } /* End ( x < +0 ) */ } /* End ( y != -1.0 ) */ } /* End ( eiy != +-Infinity ) */ } /* End of ei_pow() */ /*************************************************************** Hide internal symbols: ***************************************************************/ __mpu_hidden_decl(ei_trunc); __mpu_hidden_decl(ei_sin); __mpu_hidden_decl(ei_cos); __mpu_hidden_decl(ei_tan); __mpu_hidden_decl(ei_log1p); __mpu_hidden_decl(ei_log); __mpu_hidden_decl(ei_log10); __mpu_hidden_decl(ei_log2); __mpu_hidden_decl(ei_expm1); __mpu_hidden_decl(ei_exp); __mpu_hidden_decl(ei_atan2); __mpu_hidden_decl(ei_atan); __mpu_hidden_decl(ei_sinh); __mpu_hidden_decl(ei_cosh); __mpu_hidden_decl(ei_tanh); __mpu_hidden_decl(ei_asinh); __mpu_hidden_decl(ei_acosh); __mpu_hidden_decl(ei_atanh); __mpu_hidden_decl(ei_asin); __mpu_hidden_decl(ei_acos); __mpu_hidden_decl(ei_pow); /* End of hide internal symbols. ***************************************************************/