nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal pow(final BigDecimal x, final BigDecimal y) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            throw new ArithmeticException("Cannot power negative " + x.toString());
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else {
            /* return x^y = exp(y*log(x)) ;
             */
            BigDecimal logx = log(x);
            BigDecimal ylogx = y.multiply(logx);
            BigDecimal resul = exp(ylogx);
            /* The estimation of the relative error in the result is |log(x)*err(y)|+|y*err(x)/x|
             */
            double errR = Math.abs(logx.doubleValue() * y.ulp().doubleValue() / 2.)
                            + Math.abs(y.doubleValue() * x.ulp().doubleValue() / 2. / x.doubleValue());
            MathContext mcR = new MathContext(err2prec(1.0, errR));
            return resul.round(mcR);
        }
    }
static public BigDecimal pow(final BigDecimal x, final BigDecimal y) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            throw new ArithmeticException("Cannot power negative " + x.toString());
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else {
            /* return x^y = exp(y*log(x)) ;
             */
            BigDecimal logx = log(x);
            BigDecimal ylogx = y.multiply(logx);
            BigDecimal resul = exp(ylogx);
            /* The estimation of the relative error in the result is |log(x)*err(y)|+|y*err(x)/x|
             */
            double errR = Math.abs(logx.doubleValue() * y.ulp().doubleValue() / 2.)
                            + Math.abs(y.doubleValue() * x.ulp().doubleValue() / 2. / x.doubleValue());
            MathContext mcR = new MathContext(err2prec(1.0, errR));
            return resul.round(mcR);
        }
    }
Power function.

@param x Base of the power.
@param y Exponent of the power.
@return x^y.
The estimation of the relative error in the result is |log(x)*err(y)|+|y*err(x)/x|