#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define _x(i)       (p[i])
#define _y(i)       (q[i])
#define _x2(i)      (p[i]*p[i])
#define _y2(i)      (q[i]*q[i])
#define _xy(i)      (p[i]*q[i])

#define EXPR36(x)   x( 0)+x( 1)+x( 2)+x( 3)+x( 4)+x( 5)+x( 6)+x( 7)+x( 8)+x( 9)+x(10)+x(11)+\
                    x(12)+x(13)+x(14)+x(15)+x(16)+x(17)+x(18)+x(19)+x(20)+x(21)+x(22)+x(23)+\
                    x(24)+x(25)+x(26)+x(27)+x(28)+x(29)+x(30)+x(31)+x(32)+x(33)+x(34)+x(35)

#define n           36


static inline float
atan2i ( float x, float y )       // maximum 0.915° error between -0.82765 ... +2.39845 radian
{
    if      ( x < y ) {
        x  = x / y;
        y  = x * x;
        x *= 1. - (1 - M_PI/4) * y;
    }
    else if ( x > y ) {
        x  = y / x;
        y  = x * x;
        x *= 1. - (1 - M_PI/4) * y;
        x  = M_PI/2 - x;
    }
    else {
        x  = M_PI/4;
    }

    return x;
}


void
Regression ( float* const  _r,
             float* const  _b,
             const float*  p ,
             const float*  q )
{
    float  x  = EXPR36 (_x );
    float  y  = EXPR36 (_y );
    float  x2 = EXPR36 (_x2);
    float  y2 = EXPR36 (_y2);
    float  xy = EXPR36 (_xy);
    float  r;
    float  sx;
    float  sy;
    float  b;

    r  = (x2*n - x*x) * (y2*n - y*y);
    r  = r > 0.  ?  (xy*n - x*y) / sqrt (r)  :  1.;
    sx = sqrt ( (x2 - x*x/n) / (n - 1) );
    sy = sqrt ( (y2 - y*y/n) / (n - 1) );
    x  = x/n;
    y  = y/n;

    b  = atan2 ( sy, sx );
    if ( r < 0 )
        b = M_PI - b;

    *_r = r;
    *_b = b;

    printf ( "r=%6.3f  sx=%6.3f  sy=%6.3f  b=%6.3f\n", r, sx, sy, b );
}




/*
int
main ( void )
{
    float  p[36];
    float  q[36];
    float  r;
    float  b;
    int    i;
    int    j;

    float  w;
    float  c;
    float  s;
    float  rnd;
    float  rnd1;
    float  rnd2;

    for ( i = 0; i <= 200; i++ ) {
        w = 2 * M_PI / 200. * i;
        c = cos (w);
        s = sin (w);
        for ( j = 0; j < 36; j++ ) {
            rnd  = 2. * rand() / RAND_MAX - 1;
            p[j] = c * rnd;
            q[j] = s * rnd;
        }
        printf ( "i=%3u: ", i );
        Regression ( &r, &b, p, q );
    }
    printf ( "\n" );

    for ( i = 0; i <= 200; i++ ) {
        w = M_PI / 200. * i;
        c = cos(w);
        s = sin(w);
        for ( j = 0; j < 36; j++ ) {
            rnd1 = 2. * rand() / RAND_MAX - 1;
            rnd2 = 2. * rand() / RAND_MAX - 1;
            p[j] = rnd1;
            q[j] = rnd1 * c + rnd2 * s;
        }
        printf ( "i=%3u: (%6.4f) ", i, s/c );
        Regression ( &r, &b, p, q );
    }
    printf ( "\n" );

    return 0;
}
*/