pow(a, b) precision loss

[depler]

This code has precision loss:

define pow(a, b) {
    if (scale(b) == 0) {
        return a ^ b;
    }
    return e(b*l(a));
}

For example pow(32, 0.2) should be equal to 2, but calculated as 1.99999999...

[idealvin]

It is design by bc, we can do nothing for that, maybe..

[depler]

There is similar problem with cbrt(27) which shows 2.999999... instead of 3. Probably we can use the fact, that divider is integer (5 for 1/5, 3 for 1/3), somehow...

[depler]

Here is what I've made for roots:

define root(a, n)
{
    auto x, xk, diff, limit, nn;

    if (n == 0)
        return 1;

    if (n < 0)
        return 1 / root(a, -n);

    x = 1;
    diff = 1;
    limit = 10^(-scale);
    nn = n - 1;
    
    while (diff > limit)
    {
        xk = (nn * x + a / x ^ nn) / n;
        diff = abs(xk - x);
        x = xk;
    }
    
    return x;
}

This is based on https://www.geeksforgeeks.org/n-th-root-number/. Usage:

echo "root(32,5)" | bc -ql bc.txt
2.00000000000000000000

echo "root(7^19,2)" | bc -ql bc.txt
106765608.62643524030500013498

[depler]

And here is code for rational power:

define int(x) 
{
    auto tmp;
    
    tmp = scale;
    scale = 0;
    
    x /= 1;
    
    scale = tmp;
    return x;
}

define mod(x, y)
{
    return x - (y * int(x / y))
}

define rpow(a, x, y)
{
    auto ax, ay;

    if (scale(x) != 0 || scale(y) != 0)
    {
        print "x or y is not integer";
        halt;
    }
    
    if (y == 0)
    {
        print "y is zero";
        halt;
    }

    if (x == 0)
        return 1;

    if (x == y)
        return a;
    
    ax = abs(x);
    ay = abs(y);

    if (ax > ay && mod(ax, ay) == 0)
    {
        x = int(x / y);
        return a ^ x;
    }

    if (ay > ax && mod(ay, ax) == 0)
    {
        y = int(y / x);
        return root(a, y);
    }
    
    return root(a, y) ^ x;
}

Usage:

echo "rpow(64, 1, 6)" | bc -ql bc.txt
2.00000000000000000000

echo "rpow(64, 2, 12)" | bc -ql bc.txt
2.00000000000000000000

echo "rpow(2, -1, 2)" | bc -ql bc.txt
.70710678118654752440

echo "rpow(2, 1, -2)" | bc -ql bc.txt
.70710678118654752440

[idealvin]

@depler Thanks. Could you send a PR?

[depler]

I could, but I have found this: https://github.com/gavinhoward/bc/blob/master/gen/lib2.bc. It is already contains nearly everything, including fractional power. And that all integrated into binary. Feel free if you still want to add my code by yourself.