/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>


#include <stdio.h>
#include <math.h>
#include <string.h>
#include <dx/dx.h>
#include "_compute.h"
#include "_compputils.h"
#include "_compoper.h"

complexFloat
_dxfComputeAddComplexFloat(complexFloat x, complexFloat y) 
{
    complexFloat result;
    REAL(result) = REAL(x) + REAL(y);
    IMAG(result) = IMAG(x) + IMAG(y);

    return result;
}
complexFloat
_dxfComputeSubComplexFloat(complexFloat x, complexFloat y) 
{
    complexFloat result;
    REAL(result) = REAL(x) - REAL(y);
    IMAG(result) = IMAG(x) - IMAG(y);

    return result;
}
complexFloat
_dxfComputeNegComplexFloat(complexFloat x)
{
    complexFloat result;
    REAL(result) = -REAL(x);
    IMAG(result) = -IMAG(x);

    return result;
}
complexFloat
_dxfComputeMulComplexFloat(complexFloat x, complexFloat y) 
{
    complexFloat result;
    REAL(result) = REAL(x) * REAL(y) - IMAG(x) * IMAG(y);
    IMAG(result) = REAL(x) * IMAG(y) + IMAG(x) * REAL(y);

    return result;
}
complexFloat
_dxfComputeDivComplexFloat(complexFloat x, complexFloat y) 
{
    complexFloat result;
    float        denom;

    denom = REAL(y) * REAL(y) + IMAG(y) * IMAG(y);
    REAL(result) = REAL(x) * REAL(y) + IMAG(x) * IMAG(y);
    IMAG(result) = IMAG(x) * REAL(y) - REAL(x) * IMAG(y);
    REAL(result) /= denom;
    IMAG(result) /= denom;

    return result;
}
float
_dxfComputeAbsComplexFloat(complexFloat x)
{
    float result;

    result = REAL(x) * REAL(x) + IMAG(x) * IMAG(x);
    result = sqrt(result);

    return result;
}
float
_dxfComputeArgComplexFloat(complexFloat x)
{
    float result;

    result = atan2(IMAG(x), REAL(x));

    return result;
}
complexFloat
_dxfComputeSqrtComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = REAL(x) + _dxfComputeAbsComplexFloat(x);
    REAL(result) = sqrt(REAL(result) / 2.0);
    IMAG(result) = -REAL(x) + _dxfComputeAbsComplexFloat(x);
    IMAG(result) = sqrt(IMAG(result) / 2.0);
    if (IMAG(x) < 0)
	IMAG(result) = -IMAG(result);

    return result;
}
complexFloat
_dxfComputeLnComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = log(_dxfComputeAbsComplexFloat(x));
    IMAG(result) = _dxfComputeArgComplexFloat(x);

    return result;
}
complexFloat
_dxfComputePowComplexFloatFloat(complexFloat x, float y)
{
    complexFloat result;
    float        k;

    k = pow(_dxfComputeAbsComplexFloat(x), y);

    REAL(result) = k * cos(y * _dxfComputeArgComplexFloat(x));
    IMAG(result) = k * sin(y * _dxfComputeArgComplexFloat(x));

    return result;
}
complexFloat
_dxfComputeExpComplexFloat(complexFloat x)
{
    complexFloat result;
    float        k;

    k = exp(REAL(x));

    REAL(result) = k * cos(IMAG(x));
    IMAG(result) = k * sin(IMAG(x));

    return result;
}
complexFloat
_dxfComputePowComplexFloat(complexFloat x, complexFloat y)
{
    complexFloat result;

    result = _dxfComputeMulComplexFloat(y, _dxfComputeLnComplexFloat(x));
    result = _dxfComputeExpComplexFloat(result);

    return result;
}
complexFloat
_dxfComputeSinComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = sin(REAL(x)) * cosh(IMAG(x));
    IMAG(result) = cos(REAL(x)) * sinh(IMAG(x));

    return result;
}
complexFloat
_dxfComputeCosComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = cos(REAL(x)) * cosh(IMAG(x));
    IMAG(result) = -sin(REAL(x)) * sinh(IMAG(x));

    return result;
}
complexFloat
_dxfComputeSinhComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = sinh(REAL(x)) * cos(IMAG(x));
    IMAG(result) = cosh(REAL(x)) * sin(IMAG(x));

    return result;
}
complexFloat
_dxfComputeCoshComplexFloat(complexFloat x)
{
    complexFloat result;

    REAL(result) = cosh(REAL(x)) * cos(IMAG(x));
    IMAG(result) = sinh(REAL(x)) * sin(IMAG(x));

    return result;
}
complexFloat
_dxfComputeTanComplexFloat(complexFloat x)
{
    complexFloat result;
    float        denom;

    denom = cos(2.0*REAL(x)) + cosh(2.0*IMAG(x));
    REAL(result) =  sin(2.0*REAL(x)) / denom;
    IMAG(result) = sinh(2.0*IMAG(x)) / denom;

    return result;
}
complexFloat
_dxfComputeTanhComplexFloat(complexFloat x)
{
    complexFloat result;

    result = _dxfComputeCoshComplexFloat(x);
    if (REAL(result) == 0 && IMAG(result) == 0)
	REAL(result) = IMAG(result) = 1e20;
    else
	result = _dxfComputeDivComplexFloat(_dxfComputeSinhComplexFloat(x), result);

    return result;
}
complexFloat
_dxfComputeAsinComplexFloat(complexFloat x)
{
    complexFloat result;
    complexFloat temp;

    REAL(result) = 1;
    IMAG(result) = 0;
    temp = _dxfComputeMulComplexFloat(x, x);
    result = _dxfComputeSubComplexFloat(result, temp);
    result = _dxfComputeSqrtComplexFloat(result);
    REAL(temp) = 0;
    IMAG(temp) = 1;
    temp = _dxfComputeMulComplexFloat(temp, x);
    result = _dxfComputeAddComplexFloat(temp, result);
    result = _dxfComputeLnComplexFloat(result);
    REAL(temp) = 0;
    IMAG(temp) = -1;
    result = _dxfComputeMulComplexFloat(temp, result);

    return result;
}
complexFloat
_dxfComputeAcosComplexFloat(complexFloat x)
{
    complexFloat result;
    complexFloat temp;

    REAL(result) = 1;
    IMAG(result) = 0;
    temp = _dxfComputeMulComplexFloat(x, x);
    result = _dxfComputeSubComplexFloat(temp, result);
    result = _dxfComputeSqrtComplexFloat(result);
    result = _dxfComputeAddComplexFloat(x, result);
    result = _dxfComputeLnComplexFloat(result);
    REAL(temp) = 0;
    IMAG(temp) = -1;
    result = _dxfComputeMulComplexFloat(temp, result);

    return result;
}
/* 1/2 i ln[(1+ix) / (1-ix)] */
complexFloat
_dxfComputeAtanComplexFloat(complexFloat x)
{
    complexFloat result;
    complexFloat ix;
    complexFloat temp;

    REAL(temp) = 0;
    IMAG(temp) = 1;
    ix = _dxfComputeMulComplexFloat(temp, x);
    REAL(temp) = 1;
    IMAG(temp) = 0;
    result = _dxfComputeAddComplexFloat(temp, ix);
    temp = _dxfComputeSubComplexFloat(temp, ix);
    if (REAL(temp) != 0 || IMAG(temp) != 0)
	result = _dxfComputeDivComplexFloat(result, temp);
    result = _dxfComputeLnComplexFloat(result);
    REAL(temp) = 0;
    IMAG(temp) = -0.5;
    result = _dxfComputeMulComplexFloat(result, temp);

    return result;
}