Search code examples
cfunctiontestingpow

Checking limits of own pow function in C?


So I have am trying to implement my own pow function. To make things somewhat simpler I am allowed to use type int for the exponent. So the code is working...mostly. I have a file called xmath.c in which I implement my pow function and another test.c to test the limits. The only case in which it's not working is for pow(-111, INT_MAX). How could I improve my code ? Also tell me if I should post the test.c code as well. xfabs is another function that I had to implement on my own instead of using fabs.

Forgot to include xfabs in my post:

#define xfabs(x) (((x) > 0) ? (x) : (-(x)))

#include <stdio.h>
#include <errno.h>
#include <math.h>
#include <float.h>
#include "xmath.h"

double xpow(double x, int y)
{
    errno = 0;
    double product = 1;
    double prod = 1;
    int i;

    if (y == 0) {
        return 1.0;
    } else if (x == 0 && y < 0) {
        errno = EDOM;
        return 0.0;
    }

    if (xfabs(x) >= 1) {
        if (y > 0) {
            for (i = 0; i < y; i++) {
                if (xfabs(product) > DBL_MAX / xfabs(x)) {
                        errno = ERANGE;
                        return HUGE_VAL;
                    }
                    product *= x;
                }
            return product;
        } else if (y < 0) {
            for (i = 0; i > y; i--) {
                if (xfabs(product) < DBL_MIN / xfabs(x)) {
                    errno = ERANGE;
                    return -0.0;
                }
                product *= 1 / x;
            }
            return product;
        }

    } else if (xfabs(x) < 1) {
        if (y > 0) {
            for (i = 0; i < y; i++) {
                if (xfabs(prod) < DBL_MIN / xfabs(x)) {
                    errno = ERANGE;
                    return +0.0;
                }
                prod *= x;
            }
            return prod;
        } else if (y < 0) {
            for (i = 0; i > y; i--) {
                if (xfabs(prod) < DBL_MAX / xfabs(x)) {
                    errno = ERANGE;
                    return -HUGE_VAL;
                }
                prod *= 1 / x;
            }
            return prod;
        }
    }
}

Here is the test function. Also I am creating an object file for xmath.c which I use in test.c:

#include <stdio.h>
#include <float.h>
#include <math.h>
#include <errno.h>
#include <limits.h>

#include "xmath.h"

void xpow_tests(void)
{
    printf("XPOW TESTS (SUCCESS == 1 / FAILURE == 0)\n");
    printf("=================================================\n");
    printf("Test (xpow(2, 5)): \t\t%i\t", xpow(2, 5) == pow(2, 5));
    xpow(2, 5);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-2, -4)): \t\t%i\t", xpow(-2, -4) == pow(-2, -4));
    xpow(-2, -4);
    printf("errno: %i\n", errno);
    printf("Test (xpow(1, 15)): \t\t%i\t", xpow(1, 15) == pow(1, 15));
    xpow(1, 15);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-4, 414)): \t\t%i\t", xpow(-4, 414) == pow(-4, 414));
    xpow(-4, 414);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-5, 303)): \t\t%i\t", xpow(-5, 303) == pow(-5, 303));
    xpow(-5, 303);
    printf("errno: %i\n", errno);
    printf("Test (xpow(0, -3)): \t\t%i\t", xpow(0, -3) == 0);
    xpow(0, -3);
    printf("errno: %i\n", errno);
    printf("Test (xpow(1.0e-10, 100)): \t%i\t", xpow(1.0e-10, 100) == pow(1.0e-10, 100));
    xpow(1.0e-10, 100);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-1.0e-10, 101)): \t%i\t", xpow(-1.0e-10, 101) == pow(-1.0e-10, 101));
    xpow(-1.0e-10, 101);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-111, INT_MAX)): \t%i\t", xpow(-111, INT_MAX) == pow(-111, INT_MAX));
    xpow(-111, INT_MAX);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-111, INT_MAX-1)): \t%i\t", xpow(-111, INT_MAX-1) == pow(-111, INT_MAX-1));
    xpow(-111, INT_MAX-1);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-111, -INT_MAX)): \t%i\t", xpow(-111, -INT_MAX) == pow(-111, -INT_MAX));
    xpow(-111, -INT_MAX);
    printf("errno: %i\n", errno);
    printf("Test (xpow(-111, -INT_MAX+1)): \t%i\t", xpow(-111, -INT_MAX+1) == pow(-111, -INT_MAX+1));
    xpow(-111, -INT_MAX+1);
    printf("errno: %i\n", errno);
    printf("=================================================\n");
    printf("errno (EDOM: %i / ERANGE: %i)\n\n", EDOM, ERANGE);
}

int main(void)
{
    xpow_tests();
    return 0;
}

Solution

  • Caveat: This isn't a solution to your failing test case per se.

    But, I simplified your function a bit--this may help.

    And I greatly simplified/enhanced your test case function [using some preprocessor trickery].

    However, in addition to your failed test -111,INT_MAX, when I ran it here, the 0,-3 test reports failure.


    Here's the modified source:

    #include <stdio.h>
    #include <string.h>
    #include <errno.h>
    #include <math.h>
    #include <float.h>
    #include <limits.h>
    //#include "xmath.h"
    
    double
    xfabs(double x)
    {
    
        if (x < 0)
            x = -x;
    
        return x;
    }
    
    double
    xpow(double x, int y)
    {
        errno = 0;
        int i;
        double lim;
        double prod = 1;
        double abs_x;
    
        if (y == 0)
            return 1.0;
        if ((x == 0) && (y < 0)) {
            errno = EDOM;
            return 0.0;
        }
    
        abs_x = xfabs(x);
    
        if (abs_x >= 1) {
            if (y > 0) {
                lim = DBL_MAX / abs_x;
                for (i = 0; i < y; i++) {
                    if (xfabs(prod) > lim) {
                        errno = ERANGE;
                        return HUGE_VAL;
                    }
                    prod *= x;
                }
            }
    
            if (y < 0) {
                lim = DBL_MIN / abs_x;
                for (i = 0; i > y; i--) {
                    if (xfabs(prod) < lim) {
                        errno = ERANGE;
                        return -0.0;
                    }
                    prod *= 1 / x;
                }
            }
    
            return prod;
        }
    
        if (y > 0) {
            lim = DBL_MIN / abs_x;
            for (i = 0; i < y; i++) {
                if (xfabs(prod) < lim) {
                    errno = ERANGE;
                    return +0.0;
                }
                prod *= x;
            }
        }
    
        if (y < 0) {
            lim = DBL_MAX / abs_x;
            for (i = 0; i > y; i--) {
                if (xfabs(prod) < lim) {
                    errno = ERANGE;
                    return -HUGE_VAL;
                }
                prod *= 1 / x;
            }
        }
    
        return prod;
    }
    
    int passcnt = 0;
    int failcnt = 0;
    
    void
    dotest(const char *name,double x,int y)
    {
        int sverr_pow;
        int sverr_xpow;
        double expected;
        double actual;
        int passflg;
    
        printf("\n");
        printf("%s\n",name);
    
        errno = 0;
        expected = pow(x,y);
        sverr_pow = errno;
    
        printf("%.16g (pow)\n",expected);
        if (sverr_pow)
            printf("errno: %d (%s)\n",sverr_pow,strerror(sverr_pow));
    
        errno = 0;
        actual = xpow(x,y);
        sverr_xpow = errno;
    
        printf("%.16g (xpow)\n",actual);
        if (sverr_xpow)
            printf("errno: %d (%s)\n",sverr_xpow,strerror(sverr_xpow));
    
        if ((actual == expected) && (sverr_xpow == sverr_pow))
            passcnt += 1;
        else
            failcnt += 1;
    
        printf("Result: Value:%s errno:%s\n",
            (actual == expected) ? "PASS" : "FAIL",
            (sverr_xpow == sverr_pow) ? "PASS" : "FAIL");
    }
    
    #define TEST(x,y) \
        dotest(#x "," #y,x,y)
    
    void xpow_tests(void)
    {
        printf("XPOW TESTS (SUCCESS == 1 / FAILURE == 0)\n");
        printf("=================================================\n");
    
        TEST(2, 5);
        TEST(-2, -4);
        TEST(1, 15);
        TEST(-4, 414);
        TEST(-5, 303);
        TEST(0, -3);
        TEST(1.0e-10, 100);
        TEST(-1.0e-10, 101);
        TEST(-111, INT_MAX);
        TEST(-111, INT_MAX-1);
        TEST(-111, -INT_MAX);
        TEST(-111, -INT_MAX+1);
    
        printf("=================================================\n");
        printf("errno (EDOM: %i / ERANGE: %i)\n\n", EDOM, ERANGE);
    
        printf("PASSED: %d, FAILED: %d\n",passcnt,failcnt);
    }
    
    int
    main(void)
    {
        xpow_tests();
        return 0;
    }
    

    Here's the output of the tests:

    XPOW TESTS (SUCCESS == 1 / FAILURE == 0)
    =================================================
    
    2,5
    32 (pow)
    32 (xpow)
    Result: Value:PASS errno:PASS
    
    -2,-4
    0.0625 (pow)
    0.0625 (xpow)
    Result: Value:PASS errno:PASS
    
    1,15
    1 (pow)
    1 (xpow)
    Result: Value:PASS errno:PASS
    
    -4,414
    1.789931494904685e+249 (pow)
    1.789931494904685e+249 (xpow)
    Result: Value:PASS errno:PASS
    
    -5,303
    -6.136366831622158e+211 (pow)
    -6.136366831622158e+211 (xpow)
    Result: Value:PASS errno:PASS
    
    0,-3
    inf (pow)
    errno: 34 (Numerical result out of range)
    0 (xpow)
    errno: 33 (Numerical argument out of domain)
    Result: Value:FAIL errno:FAIL
    
    1.0e-10,100
    0 (pow)
    errno: 34 (Numerical result out of range)
    0 (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:PASS errno:PASS
    
    -1.0e-10,101
    -0 (pow)
    errno: 34 (Numerical result out of range)
    0 (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:PASS errno:PASS
    
    -111,INT_MAX
    -inf (pow)
    errno: 34 (Numerical result out of range)
    inf (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:FAIL errno:PASS
    
    -111,INT_MAX-1
    inf (pow)
    errno: 34 (Numerical result out of range)
    inf (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:PASS errno:PASS
    
    -111,-INT_MAX
    -0 (pow)
    errno: 34 (Numerical result out of range)
    -0 (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:PASS errno:PASS
    
    -111,-INT_MAX+1
    0 (pow)
    errno: 34 (Numerical result out of range)
    -0 (xpow)
    errno: 34 (Numerical result out of range)
    Result: Value:PASS errno:PASS
    =================================================
    errno (EDOM: 33 / ERANGE: 34)
    
    PASSED: 10, FAILED: 2