[comp.os.minix] A atan function

u31b3hs@cip-s01.informatik.rwth-aachen.de (Michael Haardt) (01/28/91)

From: u31b3hs%cip-s01.informatik.rwth-aachen.de@unido.de

Hello world!

A few weeks ago, someone asked for a arctan function.  Here it is.  I
use this function since the good old CP/M days, therefore I don't know
its origin.  It works with PC MINIX 1.5.10 floating point library, for
higher precision you should change some definitions.  I am still waiting
for a complete libm.a, perhaps other people find other functions in
their sources.  Post them, if they are missing in MINIX! Note to AST: I
heard that there will be official floating point support for 1.6.  What
about libm.a, is it complete or should we work to finish it?

-------- snip, snip, snip --------- not the listing!!!! ... oh boy ----------

#include <lib.h>
#include <math.h>

#define eps 1.49011611938477e-8
#define Smallest 1.0e-7
#define Biggest 1.0e15
#define MaxArctan 1.5707963267948966192
#define SqrtThree 1.7320508075688772935
#define SqrtThreeLessOne 0.7320508075688772935
#define TwoLessSqrtThree 0.26794919243112270647
#define ArcTanP0 -13.6887688941919
#define ArcTanP1 -20.5058551958617
#define ArcTanP2 -8.49462403513207
#define ArcTanP3 -0.837582993681501
#define ArcTanQ0 41.0663066825758
#define ArcTanQ1 86.1573495971302
#define ArcTanQ2 59.5784361425973
#define ArcTanQ3 15.0240011600286

static double ArctanA[4] =
{
  0.0,
  0.523598775598299,
  1.5707963267948966192,
  1.04719755119660
};

double atan(x) double x;
{
  double f,g,R,P,Q,Result;
  unsigned int n;
  
  f=fabs(x);
  if (f<Smallest) return x; /* Underflow! */
  else if (x>=Biggest) return MaxArctan; /* Overflow! */
  else if (x<=-Biggest) return -MaxArctan; /* Overflow! */
  if (f>1.0) { f=1.0/f; n=2; } else n=0;
  if (f>TwoLessSqrtThree)
  {
    f=(((SqrtThreeLessOne*f-0.5)-0.5)+f)/(SqrtThree+f);
    n++;
    /* ignore scale f step for floating-point implementation */
  }
  if (fabs(f)<eps) Result=f;
  else
  {
    g=f*f;
    /* evaluate R(g) */
    P=((ArcTanP3*g+ArcTanP2)*g+ArcTanP1)*g+ArcTanP0;
    Q=(((g+ArcTanQ3)*g+ArcTanQ2)*g+ArcTanQ1)*g+ArcTanQ0;
    R=g*P/Q;
    Result=f+f*R;
  }
  if (n>1) Result = -Result;
  Result=ArctanA[n]+Result;
  if (x<0.0) Result = -Result;
  return Result;
}

----------------------------- snip, snip again  ----------------------------

Michael Haardt

u31b3hs%cip-s01.informatik.rwth-aachen.de@unido.bitnet
--------------------------------- Namaskaar --------------------------------