nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal atan(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return atan(x.negate()).negate();


        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else if (x.doubleValue() > 0.7 && x.doubleValue() < 3.0) {
            /* Abramowitz-Stegun 4.4.34 convergence acceleration
             * 2*arctan(x) = arctan(2x/(1-x^2)) = arctan(y). x=(sqrt(1+y^2)-1)/y
             * This maps 0<=y<=3 to 0<=x<=0.73 roughly. Temporarily with 2 protectionist digits.
             */
            BigDecimal y = scalePrec(x, 2);
            BigDecimal newx = divideRound(hypot(1, y).subtract(BigDecimal.ONE), y);
            /* intermediate result with too optimistic error estimate*/
            BigDecimal resul = multiplyRound(atan(newx), 2);
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp. */


            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        } else if (x.doubleValue() < 0.71) {
            /* Taylor expansion around x=0; Abramowitz-Stegun 4.4.42 */
            final BigDecimal xhighpr = scalePrec(x, 2);
            final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
            BigDecimal resul = xhighpr.plus();
            /* signed x^(2i+1) */
            BigDecimal xpowi = xhighpr;
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
             */


            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));


            for (int i = 1;; i++) {
                xpowi = multiplyRound(xpowi, xhighprSq);
                BigDecimal c = divideRound(xpowi, 2 * i + 1);
                resul = resul.add(c);


                if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                    break;
                }


            }
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        } else {
            /* Taylor expansion around x=infinity; Abramowitz-Stegun 4.4.42 */
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
             */
            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
            /* start with the term pi/2; gather its precision relative to the expected result
             */
            MathContext mc = new MathContext(2 + err2prec(3.1416, eps));
            BigDecimal onepi = pi(mc);
            BigDecimal resul = onepi.divide(new BigDecimal(2));
            final BigDecimal xhighpr = divideRound(-1, scalePrec(x, 2));
            final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
            /* signed x^(2i+1) */
            BigDecimal xpowi = xhighpr;


            for (int i = 0;; i++) {
                BigDecimal c = divideRound(xpowi, 2 * i + 1);
                resul = resul.add(c);


                if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                    break;
                }
                xpowi = multiplyRound(xpowi, xhighprSq);


            }
            mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        }
    }
static public BigDecimal atan(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return atan(x.negate()).negate();


        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else if (x.doubleValue() > 0.7 && x.doubleValue() < 3.0) {
            /* Abramowitz-Stegun 4.4.34 convergence acceleration
             * 2*arctan(x) = arctan(2x/(1-x^2)) = arctan(y). x=(sqrt(1+y^2)-1)/y
             * This maps 0<=y<=3 to 0<=x<=0.73 roughly. Temporarily with 2 protectionist digits.
             */
            BigDecimal y = scalePrec(x, 2);
            BigDecimal newx = divideRound(hypot(1, y).subtract(BigDecimal.ONE), y);
            /* intermediate result with too optimistic error estimate*/
            BigDecimal resul = multiplyRound(atan(newx), 2);
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp. */


            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        } else if (x.doubleValue() < 0.71) {
            /* Taylor expansion around x=0; Abramowitz-Stegun 4.4.42 */
            final BigDecimal xhighpr = scalePrec(x, 2);
            final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
            BigDecimal resul = xhighpr.plus();
            /* signed x^(2i+1) */
            BigDecimal xpowi = xhighpr;
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
             */


            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));


            for (int i = 1;; i++) {
                xpowi = multiplyRound(xpowi, xhighprSq);
                BigDecimal c = divideRound(xpowi, 2 * i + 1);
                resul = resul.add(c);


                if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                    break;
                }


            }
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        } else {
            /* Taylor expansion around x=infinity; Abramowitz-Stegun 4.4.42 */
            /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
             */
            double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
            /* start with the term pi/2; gather its precision relative to the expected result
             */
            MathContext mc = new MathContext(2 + err2prec(3.1416, eps));
            BigDecimal onepi = pi(mc);
            BigDecimal resul = onepi.divide(new BigDecimal(2));
            final BigDecimal xhighpr = divideRound(-1, scalePrec(x, 2));
            final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
            /* signed x^(2i+1) */
            BigDecimal xpowi = xhighpr;


            for (int i = 0;; i++) {
                BigDecimal c = divideRound(xpowi, 2 * i + 1);
                resul = resul.add(c);


                if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                    break;
                }
                xpowi = multiplyRound(xpowi, xhighprSq);


            }
            mc = new MathContext(err2prec(resul.doubleValue(), eps));


            return resul.round(mc);


        }
    }
The inverse trigonometric tangent.

@param x the argument.
@return the principal value of arctan(x) in radians in the range -pi/2 to +pi/2.