pow.c

来自「OMAP1030 处理器的ARM 侧硬件测试代码 OMAP1030 是TI」· C语言 代码 · 共 66 行

/****************************************************************************/
/*  pow    v2.24                                                           */
/*  Copyright (c) 1995-2002 Texas Instruments Incorporated                  */
/****************************************************************************/
#include <math.h>
#include <values.h>
#include <errno.h>

double _log(double x);                 /* Like log() with no error checking */

/****************************************************************************/
/*  POW() - Power							    */
/*									    */
/*  z = mantissa x							    */
/*  log2(z) = c1 * z ^ 9 + c2 * z ^ 7 + c3 * z ^ 5 + c4 * z ^ 3 + c5 * z    */
/*  log2(x) = exponent x + log2(z)					    */
/*  a = y * log2(x)							    */
/*  b = mantissa a							    */
/*  2 ^ b =  (((((d1 * b + d2) * b + d3) * b + d4) * b + d5) * b + d6)	    */
/*									    */
/*  result = 2 ^ exponent a * 2 ^ b					    */
/*	   = 2 ^ (y * log2(x))						    */
/****************************************************************************/
double pow(double x, double y)
{
    int z;

    if (x && y && y != 1.0)
    {
        int sign = 0;
	if (x < 0.0)
	{
	    z = (int) y;

	    /*****************************************************************/
	    /* if y is not an integer, a domain error occurs                 */
	    /*****************************************************************/
	    if (y - z) { errno = EDOM; return 0.0; }

	    /*****************************************************************/
	    /* if x < 0, compute the power of |x|                            */
	    /*****************************************************************/
	    x = -x;

	    /*****************************************************************/
	    /* for odd exponents, negate the answer                          */
	    /*****************************************************************/
	    sign = z % 2; 
	}

	x = exp(_log(x) * y);

	return (sign) ? -x : x;
    }

    if (y == 1.0) return (x);

    /**************************************************************************/
    /* if x = 0 and y < 0, return a domain error                              */
    /**************************************************************************/
    if (x == 0.0 && y <= 0.0) { errno = EDOM; return (0.0); }

    return (y == 0.0) ? 1.0 : 0.0;
}