atanFast function

double atanFast(double xa, [ double xb = 0.0, bool leftPlane = false ])

Internal helper function to compute arctangent.

• xa number from which arctangent is requested.
• xb extra bits for x (may be 0.0).
• leftPlane if true, result angle must be put in the left half plane. Returns atan(xa + xb) (or angle shifted by PI if leftPlane is true)

Implementation

double atanFast(double xa, [double xb = 0.0, bool leftPlane = false]) {
  if (xa == 0.0) {
    // Matches +/- 0.0; return correct sign
    return leftPlane ? copySign(math.pi, xa) : xa;
  }

  bool negate;
  if (xa < 0) {
    // negative
    xa = -xa;
    xb = -xb;
    negate = true;
  } else {
    negate = false;
  }

  if (xa > 1.633123935319537E16) {
    // Very large input
    return (negate != leftPlane) ? (-math.pi * F_1_2) : (math.pi * F_1_2);
  }

  /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
  int idx;
  if (xa < 1) {
    idx = (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5).toInt();
  } else {
    final oneOverXa = 1 / xa;
    idx = (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) +
            13.07)
        .toInt();
  }

  final ttA = TANGENT_TABLE_A[idx];
  final ttB = TANGENT_TABLE_B[idx];

  var epsA = xa - ttA;
  var epsB = -(epsA - xa + ttA);
  epsB += xb - ttB;

  var temp = epsA + epsB;
  epsB = -(temp - epsA - epsB);
  epsA = temp;

  /* Compute eps = eps / (1.0 + xa*tangent) */
  temp = xa * HEX_40000000;
  var ya = xa + temp - temp;
  var yb = xb + xa - ya;
  xa = ya;
  xb += yb;

  //if (idx > 8 || idx == 0)
  if (idx == 0) {
    /// If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice
    //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
    final denom = 1.0 / (1.0 + (xa + xb) * (ttA + ttB));
    //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
    ya = epsA * denom;
    yb = epsB * denom;
  } else {
    var temp2 = xa * ttA;
    var za = 1.0 + temp2;
    var zb = -(za - 1.0 - temp2);
    temp2 = xb * ttA + xa * ttB;
    temp = za + temp2;
    zb += -(temp - za - temp2);
    za = temp;

    zb += xb * ttB;
    ya = epsA / za;

    temp = ya * HEX_40000000;
    final yaa = (ya + temp) - temp;
    final yab = ya - yaa;

    temp = za * HEX_40000000;
    final zaa = (za + temp) - temp;
    final zab = za - zaa;

    /* Correct for rounding in division */
    yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;

    yb += -epsA * zb / za / za;
    yb += epsB / za;
  }

  epsA = ya;
  epsB = yb;

  // Evaluate polynomial
  final epsA2 = epsA * epsA;

  /*
  yb = -0.09001346640161823;
  yb = yb * epsA2 + 0.11110718400605211;
  yb = yb * epsA2 + -0.1428571349122913;
  yb = yb * epsA2 + 0.19999999999273194;
  yb = yb * epsA2 + -0.33333333333333093;
  yb = yb * epsA2 * epsA;
  */

  yb = 0.07490822288864472;
  yb = yb * epsA2 - 0.09088450866185192;
  yb = yb * epsA2 + 0.11111095942313305;
  yb = yb * epsA2 - 0.1428571423679182;
  yb = yb * epsA2 + 0.19999999999923582;
  yb = yb * epsA2 - 0.33333333333333287;
  yb = yb * epsA2 * epsA;

  ya = epsA;

  temp = ya + yb;
  yb = -(temp - ya - yb);
  ya = temp;

  /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
  yb += epsB / (1.0 + epsA * epsA);

  final eighths = EIGHTHS[idx];

  //result = yb + eighths[idx] + ya;
  var za = eighths + ya;
  var zb = -(za - eighths - ya);
  temp = za + yb;
  zb += -(temp - za - yb);
  za = temp;

  var result = za + zb;

  if (leftPlane) {
    // Result is in the left plane
    final resultb = -(result - za - zb);
    final pia = 1.5707963267948966 * 2;
    final pib = 6.123233995736766E-17 * 2;

    za = pia - result;
    zb = -(za - pia + result);
    zb += pib - resultb;

    result = za + zb;
  }

  if (negate != leftPlane) {
    result = -result;
  }

  return result;
}