Accurate arctan()

[Duthomhas]

I've been playing around with bigfloats and Carlson's method of accurately calculating arctangent values.

So here are my questions:

- Does anyone have any reason to use arctan to more than 12 places (decimal) accuracy? (Because the algorithm gets exponential pretty quickly -- more on this in a moment.)

- Does anyone know of an online resource that gives N digits of arctan(x) (for whatever x, I'll work with it) up in at least the thousands range? (I could calculate them myself, but I don't know where the accuracy lies past what Calc can offer, and I don't trust Calc after about 15-20 digits.)


If you are interested, you can find Carlson's algorithm here:
http://dx.doi.org/10.1090/S0025-5718-1972-0307438-2

As listed, it is exponentially expensive in two ways:
(1) stack space (for the binary recursions)
(2) heap space (for the bigfloat allocations)
(Think computer trees blossoming. Or mushroom clouds rising over a city. Either thought will do.)


After memoizing, it remains exponentially expensive in only one way: think bubble-sort. The precision of the algorithm, n, directly affects the number of iterations to compute, which are n!. Yay, not.

(And additional heap space is n+c. Stack space is O(1).)

But at least it is relatively fast now. For a 100 bit big number it computes in about .385 seconds, using n=4 (24 iterations).


Thanks for reading.

[Duthomhas]

A reference implementation for y'all to play with. :O)

#include <ciso646>
#include <cmath>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <vector>
using namespace std;

//----------------------------------------------------------------------------
// The recursive way

long double an( long double x, unsigned n )
  {
  long double a0 = 1;
  long double g0 = sqrt( 1 + x*x );
  while (n--)
    {
    a0 = 0.5 * (a0 + g0);
    g0 = sqrt( a0*g0 );
    }
  return a0;
  }

long double d( long double x, unsigned k, unsigned n )
  {
  if (k == 0) return an( x, n );
  long double d0 = d( x, k-1, n   );
  long double d1 = d( x, k-1, n-1 );
  long double p = 1/pow( 2.0, 2 * k );
  return (d0 - (p * d1)) / (1 - p);
  }

long double arctan_r( long double x, unsigned n = 10 )
  {
  return x / d( x, n, n );
  }

//----------------------------------------------------------------------------
// Now to tweak it up fine

long double arctan( long double x, unsigned N = 10 )
  {
  if (N < 2) N = 1;

  long double a0, g0;            // ongoing calculation of a(n)
  vector <long double> d( N+1 ); // results of d(k,n), where d[k] == d(k,n) for largest n

  g0 = sqrt( 1 + (x * x) );
  a0 = 0.5 * (1 + g0);
  g0 = sqrt( a0 * g0 );
  d[ 0 ] = a0;                  // d[0] == d(0,1) == a(1)
  d[ 1 ] = (a0 - 0.25) / 0.75;  // d[1] == d(1,1)

  // for n = 2 to N
  for (unsigned n = 1; n++ < N; )
    {
    // calc a(n)
    a0 = 0.5 * (a0 + g0);
    g0 = sqrt( a0 * g0 );

    long double d0 = a0;

    // for k = 1 to n
    for (unsigned k = 0; k++ < n; )
      {
      // d0     = d(k-1,n)
      // d[k-1] = d(k-1,n-1)
      // p      = 1/(2**(2*k))
      // d(k,n) = (d0 - (p * d[k-1])) / (1 - p)
      long double p = 1 / pow( 2.0, 2*k );
      long double d1 = (d0 - (p * d[ k-1 ])) / (1 - p);

      d[ k-1 ] = d0;
      d0       = d1;
      }
    d[ n ] = d0;
    }

  return x / d[ N ];
  }

//----------------------------------------------------------------------------
int main( int argc, char** argv )
  {
  double x;
  int n = 8;
  if (argc == 3)
    {
    string s( argv[ 2 ] );
    istringstream ss( s );
    ss >> n;
    argc = 2;
    }
  if (argc == 2)
    {
    string s( argv[ 1 ] );
    istringstream ss( s );
    ss >> x;
    if (!ss) argc = 1;
    }
  if (argc != 2)
    {
    cout << "calculate atan(x); input x> ";
    cin >> x;
    if (!cin) { cout << "fooey!\n"; return 1; }
    cin >> n;
    }

  cout << setprecision( numeric_limits <long double> ::digits10 + 1 );
  cout << "recursive = " << arctan_r( x, n ) << "\n";
  cout << "iterative = " << arctan( x, n ) << "\n";
  cout << "standard  = " << atan( x ) << "\n";

  cout << "test iterations? ";
  unsigned long iterations;
    {
    string s;
    getline( cin, s );
    istringstream ss( s );
    ss >> iterations;
    if (ss and iterations)
      {
      clock_t time;

      time = clock();
      for (unsigned long n = 0; n < iterations; n++)
        arctan( x, n );
      time = clock() - time;
      cout << "iterative * " << iterations << " = "
           << ((long double)time/CLOCKS_PER_SEC) << " seconds\n";

      time = clock();
      for (unsigned long n = 0; n < iterations; n++)
        arctan_r( x, n );
      time = clock() - time;
      cout << "recursive * " << iterations << " = "
           << ((long double)time/CLOCKS_PER_SEC) << " seconds\n";
      }
    }

  return 0;
  }

Oh yeah, g++ -Wall -ansi -pedantic arctan.cpp -lm

[chrisname]

WolframAlpha will be able to help you, although I don't know if you can get thousands of digits without having a paid subscription.

[ResidentBiscuit]

WolframAlpha seems to be able to do anything math super quick. Kind of blows my mind. Seems to have no issues getting hundreds of arctan(x) digits.

[chrisname]

I love WolframAlpha. My favourite searches are for nutrition information of foods in ridiculous quantities, for example "one lunar mass of cheese nutrition" tells you how many calories you would get if the moon was made of cheese and you ate it (2.147x10²⁶ kcal).

[N36fSL3A]

chrisname wrote:

I love WolframAlpha. My favourite searches are for nutrition information of foods in ridiculous quantities, for example "one lunar mass of cheese nutrition" tells you how many calories you would get if the moon was made of cheese and you ate it (2.147x1026 kcal).

http://prntscr.com/1gv7pm

Seems cool. What cheese did you specify?

[ResidentBiscuit]

chrisname wrote:

I love WolframAlpha. My favourite searches are for nutrition information of foods in ridiculous quantities, for example "one lunar mass of cheese nutrition" tells you how many calories you would get if the moon was made of cheese and you ate it (2.147x1026 kcal).

This is just hilarious. Might have to show this around...

[Duthomhas]

Dang... I need to pay money to get something other than a picture of the first 100000 digits of arctan(1).

[helios]

Does anyone have any reason to use arctan to more than 12 places (decimal) accuracy? (Because the algorithm gets exponential pretty quickly -- more on this in a moment.)

Maybe for measuring distances to stars through parallax?

[chrisname]

@Duoas
I guess you could always transcribe them by hand? Or hire some third world children to do it?

[Duthomhas]

Maybe for measuring distances to stars through parallax?

Crap. Well, I don't suppose anyone will be using my code for that.

I guess you could always transcribe them by hand? Or hire some third world children to do it?

Alas...