/* 
 * Math2.c - misc. mathmatical functions
 * 
 * Author:  L. Van Warren
 * Date:    April 1983
 * Copyright (c) 1984 L. Van Warren
 * 
 * Last Modified: April 1989
 * Copyright (c) 1989 The California Institute of Technology (CIT).
 *           All Rights Reserved. U.S. Government Sponsorship under
 *           NASA Contract NAS7-918 is acknowledged.
 */

#include <std.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "VMisc.h"
#include "ComplexMath.h"
#include "Math2.h"

/*****************************************************************
 * TAG( xToN )
 * 
 * Raise a REALing poInt number to an Integer power.
 * Inputs:
 *  A REALing poInt number and an Integer power.
 * Outputs:
 *  A REALing poInt number result.
 * Assumptions:
 *  [None]
 * Algorithm:
 *  [None]
 */

Real
xToN (Real x, Int n)
{
    Real  answer;
    Real  xp;

    if(n == 0) return(1.0);
    for (answer = 1.0, xp = x; n > 1; n >>= 1)
    {
    if (n & 1)
        answer *= xp;
    xp *= xp;
    }

    answer *= xp;

    return (answer);
}

/*****************************************************************
 * TAG( iToN )
 * 
 * Raise an Integer number to an Integer power.
 * Inputs:
 *  A Integer number and an Integer power.
 * Outputs:
 *  A Integer number result.
 */

Int
iToN (Int i, Int n)
{
    Int  answer;
    Int  xp;

    if(n == 0) return(1);
    for (answer = 1, xp = i; n > 1; n >>= 1)
    {
    if (n & 1)
        answer *= xp;
    xp *= xp;
    }

    answer *= xp;

    return (answer);
}

Int
signOf(Real a)
{
    if (a >= 0.0)
    return (1);
    else
    return (0);
}

/*****************************************************************
 * TAG( findParam )
 * 
 * Find the value of the parameter t given three scalars a, b, and c.
 * If c = b then t = 1.0, else if c = a then t = 0.0.
 * 
 * Inputs:
 *  Two bounding scalars "a" and "b" and an Interior scalar "c".
 * Outputs:
 *  A value of t between 0.0 and 1.0.
 * Assumptions:
 *  [None]
 * Algorithm:
 *  [None]
 */

Real
findParam(Real a, Real b, Real c)
{
    return((c - a)/(b - a));
}

/* ***************************** */
/* Random numbers.               */
/* ***************************** */

/*****************************************************************
 * TAG( BooleanRand )
 * TAG( uniformRandom )
 * TAG( GaussianRand )
 * 
 * Random number generators.
 * Inputs:
 *  [None]
 * Outputs:
 *  BooleanRand: a uniformly distributed boolean random number
 *      x == 0 || x == 1
 *  uniformRandom: a uniformly distributed random number with 0 mean between
 *      -0.5 <= x <= 0.5.
 *  Gaussian: a Gaussian random number with 0 mean between
 *      -0.5 <= x <= 0.5.
 * Assumptions:
 *  Depends on random().
 *      Gaussian uses the central limit thereom with n = 6.
 * Algorithm:
 *  [None]
 */

Int
BooleanRand(void)
{
    return (IRND(uniformRandom() + 0.5));
}

Real
uniformRandom(void)
{
    Real result, numer, denom;

    numer = REAL(rand());
    denom = 32767.0;

    result = numer/denom - 0.5;

    return (result);
}

Real
GaussianRand(void)
{
    Int i;
    Real result;

    for (result = 0.0, i = 0; i < 6; i++)
    {
        result += uniformRandom();
    }

    return(result/6.0);
}

/* ************************************************************* */
/* A couple of useful functions  for reallocation.               */
/* ************************************************************* */

/*****************************************************************
 * TAG( isPwrOfTwo )
 * 
 * Is a number a power of two predicate.
 * Inputs:
 *  A number.
 * Outputs:
 *  TRUE if number is a power of two.
 *  FALSE otherwise.
 * Assumptions:
 *  [None]
 * Algorithm:
 *  Easier done than said.
 */

Int
isPwrOfTwo (Int number)
{
    Int     testval;

    testval = 1;
    while (testval < number)
    testval <<= 1;

    return (testval == number ? 1 : 0);
}

/*****************************************************************
 * TAG( nextPwrOfTwo )
 * 
 * Compute the next greater power of two.
 * Inputs:
 *  A number.
 * Outputs:
 *  The next greater power of two.
 * Assumptions:
 *  [None]
 * Algorithm:
 *  Easier done than said.
 */
 
Int
nextPwrOfTwo (Int number)
{
    Int     testval;

    testval = 1;
    while (testval <= number)
    testval <<= 1;

    return (testval);
}


/* ****************** */
/* Quadratic Equation */
/* ****************** */

QuadRoots
solveQuadratic(Real A, Real B, Real C)
{
    Real discriminant, discriminantRoot, twoA, doubleRoot;
    QuadRoots roots;
    
    if(APX_EQ(A, 0.0))
    {
        fprintf(stderr, "solveQuadratic: A is 0.0, not a quadratic form\n");
        roots.type = NOT_QUADRATIC_FORM;
    }
    else
    {
        discriminant = B*B - 4.0*A*C;
        twoA = (2.0*A);

        if(discriminant < 0.0)
        {
            roots.type = TWO_COMPLEX_ROOTS;
            roots.x[0].Re = -B/(twoA);
            roots.x[0].Im = sqrt(-discriminant);

            roots.x[1].Re = -B/(twoA);
            roots.x[1].Im = -sqrt(-discriminant);
        }
        else if(APX_EQ(discriminant, 0.0))
        {
            roots.type = DOUBLE_ROOT;
            doubleRoot = -B/(twoA);
            roots.x[0].Re = doubleRoot;
            roots.x[0].Im = 0.0;
            roots.x[1].Re = doubleRoot;
            roots.x[1].Im = 0.0;
        }
        else
        {
            roots.type = TWO_REAL_ROOTS;
            discriminantRoot = sqrt(discriminant);
            roots.x[0].Re = (-B + discriminantRoot)/twoA;
            roots.x[0].Im = 0.0;
            roots.x[1].Re = (-B - discriminantRoot)/twoA;
            roots.x[1].Im = 0.0;
        }
    }

    return(roots);
}

QuadRoots fscanQuadRoots(FILE * f)
{
    QuadRoots ret ;

    fscanf (f, "%d", &ret.type) ;
    {
    Int i ;

    for (i=0; i<2; i++)
    {
        fscanf (f, "%lf", &ret.x[i]) ;
    }
    }

    return ret ;
}

void fprIntQuadRoots (FILE* f, QuadRoots inst)
{
    fprintf (f, "%d", inst.type) ;
    fprintf (f, "\t") ;
    {
    Int j ;

    for (j=0; j<2; j++)
    {
        fprintf (f, "%lf", inst.x[j]) ;
        if (j < (2 - 1)) fprintf (f, "\t") ;
    }
    }
    NEWLINE(f);
}

QuadRoots freadQuadRoots(FILE* f)
{
    QuadRoots ret ;

    fread (&ret, sizeof(ret), 1, f) ;
    return ret ;
}

void fwriteQuadRoots (FILE * f, QuadRoots inst)
{
    fwrite (&inst, sizeof(inst), 1, f) ;
}