// Written in the D programming language

/* This module is obsolete and will be removed eventually */

/*
 * Copyright (c) 2002
 * Pavel "EvilOne" Minayev
 *
 * Permission to use, copy, modify, distribute and sell this software
 * and its documentation for any purpose is hereby granted without fee,
 * provided that the above copyright notice appear in all copies and
 * that both that copyright notice and this permission notice appear
 * in supporting documentation.  Author makes no representations about
 * the suitability of this software for any purpose. It is provided
 * "as is" without express or implied warranty.
 */
 
module std.math2;
private import std.math, std.string, std.c.stdlib, std.c.stdio;

//debug=math2;

/****************************************
 * compare floats with given precision
 */

bool feq(real a, real b)
{
	return feq(a, b, 0.000001);
} 
 
bool feq(real a, real b, real eps)
{
	return abs(a - b) <= eps;
}

/*********************************
 * Modulus
 */
 
int abs(int n)
{
	return n > 0 ? n : -n;
}
 
long abs(long n)
{
	return n > 0 ? n : -n;
}

real abs(real n)
{
	// just let the compiler handle it
	return std.math.fabs(n);
}

/*********************************
 * Square
 */

int sqr(int n)
{
	return n * n;
}
 
long sqr(long n)
{
	return n * n;
} 
 
real sqr(real n)
{
	return n * n;
}

unittest
{
	assert(sqr(sqr(3)) == 81);
}

private ushort fp_cw_chop = 7999;

/*********************************
 * Integer part
 */
 
real trunc(real n)
{
	ushort cw;
	asm
	{
		fstcw cw;
		fldcw fp_cw_chop;
		fld n;
		frndint;
		fldcw cw;
	}	
}

unittest
{
	assert(feq(trunc(+123.456), +123.0L));
	assert(feq(trunc(-123.456), -123.0L));
}

/*********************************
 * Fractional part
 */
 
real frac(real n)
{
	return n - trunc(n);
}

unittest
{
	assert(feq(frac(+123.456), +0.456L));
	assert(feq(frac(-123.456), -0.456L));
}

/*********************************
 * Sign
 */

int sign(int n)
{
	return (n > 0 ? +1 : (n < 0 ? -1 : 0));
}

unittest
{
	assert(sign(0) == 0);
	assert(sign(+666) == +1);
	assert(sign(-666) == -1);
}
 
int sign(long n)
{
	return (n > 0 ? +1 : (n < 0 ? -1 : 0));
}

unittest
{
	assert(sign(0) == 0);
	assert(sign(+666L) == +1);
	assert(sign(-666L) == -1);
}
 
int sign(real n)
{
	return (n > 0 ? +1 : (n < 0 ? -1 : 0));
}

unittest
{
	assert(sign(0.0L) == 0);
	assert(sign(+123.456L) == +1);
	assert(sign(-123.456L) == -1);
}

/**********************************************************
 * Cycles <-> radians <-> grads <-> degrees conversions
 */
 
real cycle2deg(real c)
{
	return c * 360;
}

real cycle2rad(real c)
{
	return c * PI * 2;
}

real cycle2grad(real c)
{
	return c * 400;
}

real deg2cycle(real d)
{
	return d / 360;
}

real deg2rad(real d)
{
	return d / 180 * PI;
}

real deg2grad(real d)
{
	return d / 90 * 100;
}

real rad2deg(real r)
{
	return r / PI * 180;
}

real rad2cycle(real r)
{
	return r / (PI * 2);
}

real rad2grad(real r)
{
	return r / PI * 200;
}

real grad2deg(real g)
{
	return g / 100 * 90;
}

real grad2cycle(real g)
{
	return g / 400;
}

real grad2rad(real g)
{
	return g / 200 * PI;
}

unittest
{
	assert(feq(cycle2deg(0.5), 180));
	assert(feq(cycle2rad(0.5), PI));
	assert(feq(cycle2grad(0.5), 200));
	assert(feq(deg2cycle(180), 0.5));
	assert(feq(deg2rad(180), PI));
	assert(feq(deg2grad(180), 200));
	assert(feq(rad2deg(PI), 180));
	assert(feq(rad2cycle(PI), 0.5));
	assert(feq(rad2grad(PI), 200));
	assert(feq(grad2deg(200), 180));
	assert(feq(grad2cycle(200), 0.5));
	assert(feq(grad2rad(200), PI));
}

/************************************
 * Arithmetic average of values
 */
 
real avg(real[] n)
{
	real result = 0;
	for (uint i = 0; i < n.length; i++)
		result += n[i];
	return result / n.length;
}

unittest
{
	static real[4] n = [ 1, 2, 4, 5 ];
	assert(feq(avg(n), 3));
}

/*************************************
 * Sum of values
 */

int sum(int[] n)
{
	long result = 0;
	for (size_t i = 0; i < n.length; i++)
	    result += n[i];
	return cast(int)result;
}

unittest
{
	static int[3] n = [ 1, 2, 3 ];
	assert(sum(n) == 6);
}
 
long sum(long[] n)
{
	long result = 0;
	for (uint i = 0; i < n.length; i++)
		result += n[i];
	return result;
}

unittest
{
	static long[3] n = [ 1, 2, 3 ];
	assert(sum(n) == 6);
}
 
real sum(real[] n)
{
	real result = 0;
	for (uint i = 0; i < n.length; i++)
		result += n[i];
	return result;
}

unittest
{
	static real[3] n = [ 1, 2, 3 ];
	assert(feq(sum(n), 6));
}

/*************************************
 * The smallest value
 */

int min(int[] n)
{
	int result = int.max;
	for (uint i = 0; i < n.length; i++)
		if (n[i] < result)
			result = n[i];
	return result;
}

unittest
{
	static int[3] n = [ 2, -1, 0 ];
	assert(min(n) == -1);
}
 
long min(long[] n)
{
	long result = long.max;
	for (uint i = 0; i < n.length; i++)
		if (n[i] < result)
			result = n[i];
	return result;
}

unittest
{
	static long[3] n = [ 2, -1, 0 ];
	assert(min(n) == -1);
}

real min(real[] n)
{
	real result = real.max;
	for (uint i = 0; i < n.length; i++)
	{
		if (n[i] < result)
			result = n[i];
	}
	return result;
}

unittest
{
	static real[3] n = [ 2.0, -1.0, 0.0 ];
	assert(feq(min(n), -1));
}

int min(int a, int b)
{
	return a < b ? a : b;
}

unittest
{
	assert(min(1, 2) == 1);
}

long min(long a, long b)
{
	return a < b ? a : b;
}

unittest
{
	assert(min(1L, 2L) == 1);
}

real min(real a, real b)
{
	return a < b ? a : b;
}

unittest
{
	assert(feq(min(1.0L, 2.0L), 1.0L));
}

/*************************************
 * The largest value
 */

int max(int[] n)
{
	int result = int.min;
	for (uint i = 0; i < n.length; i++)
		if (n[i] > result)
			result = n[i];
	return result;
}

unittest
{
	static int[3] n = [ 0, 2, -1 ];
	assert(max(n) == 2);
}
 
long max(long[] n)
{
	long result = long.min;
	for (uint i = 0; i < n.length; i++)
		if (n[i] > result)
			result = n[i];
	return result;
}

unittest
{
	static long[3] n = [ 0, 2, -1 ];
	assert(max(n) == 2);
}

real max(real[] n)
{
	real result = real.min;
	for (uint i = 0; i < n.length; i++)
		if (n[i] > result)
			result = n[i];
	return result;
}

unittest
{
	static real[3] n = [ 0.0, 2.0, -1.0 ];
	assert(feq(max(n), 2));
}

int max(int a, int b)
{
	return a > b ? a : b;
}

unittest
{
	assert(max(1, 2) == 2);
}

long max(long a, long b)
{
	return a > b ? a : b;
}

unittest
{
	assert(max(1L, 2L) == 2);
}

real max(real a, real b)
{
	return a > b ? a : b;
}

unittest
{
	assert(feq(max(1.0L, 2.0L), 2.0L));
}

/*************************************
 * Arccotangent
 */

real acot(real x)
{
	return tan(1.0 / x);
}

unittest
{
	assert(feq(acot(cot(0.000001)), 0.000001));
}

/*************************************
 * Arcsecant
 */

real asec(real x)
{
	return std.math.cos(1.0 / x);
}


/*************************************
 * Arccosecant
 */

real acosec(real x)
{
	return std.math.sin(1.0 / x);
}

/*************************************
 * Tangent
 */

/+
real tan(real x)
{
	asm
	{
		fld x;
		fptan;
		fstp ST(0);
		fwait;
	}
}

unittest
{
	assert(feq(tan(PI / 3), std.math.sqrt(3)));
}
+/

/*************************************
 * Cotangent
 */

real cot(real x)
{
	asm
	{
		fld x;
		fptan;
		fdivrp;
		fwait;
	}
}

unittest
{
	assert(feq(cot(PI / 6), std.math.sqrt(3.0L)));
}

/*************************************
 * Secant
 */

real sec(real x)
{
	asm
	{
		fld x;
		fcos;
		fld1;
		fdivrp;
		fwait;
	}
}


/*************************************
 * Cosecant
 */

real cosec(real x)
{
	asm
	{
		fld x;
		fsin;
		fld1;
		fdivrp;
		fwait;
	}
}

/*********************************************
 * Extract mantissa and exponent from float
 */

/+
real frexp(real x, out int exponent)
{
	asm
	{
		fld x;
		mov EDX, exponent;
		mov dword ptr [EDX], 0;
		ftst;
		fstsw AX;
		fwait;
		sahf;
		jz done;
		fxtract;
		fxch;	
		fistp dword ptr [EDX];
		fld1;
		fchs;
		fxch;
		fscale;
		inc dword ptr [EDX];
		fstp ST(1);
done:
		fwait;
	}
}

unittest
{
	int exponent;
	real mantissa = frexp(123.456, exponent);
	assert(feq(mantissa * pow(2.0L, cast(real)exponent), 123.456));
}
+/

/*************************************
 * Hyperbolic cotangent
 */

real coth(real x)
{
	return 1 / tanh(x);
}

unittest
{
    // @@@BUG1637@@@
    //assert(feq(coth(1), cosh(1) / sinh(1)));
}

/*************************************
 * Hyperbolic secant
 */

real sech(real x)
{
	return 1 / cosh(x);
}

/*************************************
 * Hyperbolic cosecant
 */

real cosech(real x)
{
	return 1 / sinh(x);
}

/*************************************
 * Hyperbolic arccosine
 */

/+
real acosh(real x)
{
	if (x <= 1)
		return 0;
	else if (x > 1.0e10)
		return log(2) + log(x);
	else
		return log(x + std.math.sqrt((x - 1) * (x + 1)));
}

unittest
{
	assert(acosh(0.5) == 0);
	assert(feq(acosh(cosh(3)), 3));
}
+/

/*************************************
 * Hyperbolic arcsine
 */

/+
real asinh(real x)
{
	if (!x)
		return 0;
	else if (x > 1.0e10)
		return log(2) + log(1.0e10);
	else if (x < -1.0e10)
		return -log(2) - log(1.0e10);
	else
	{
		real z = x * x;
		return x > 0 ? 
			log1p(x + z / (1.0 + std.math.sqrt(1.0 + z))) :
			-log1p(-x + z / (1.0 + std.math.sqrt(1.0 + z)));
	}
}

unittest
{
	assert(asinh(0) == 0);
	assert(feq(asinh(sinh(3)), 3));
}
+/

/*************************************
 * Hyperbolic arctangent
 */
/+
real atanh(real x)
{
	if (!x)
		return 0;
	else
	{
		if (x >= 1)
			return real.max;
		else if (x <= -1)
			return -real.max;
		else
			return x > 0 ?
				0.5 * log1p((2.0 * x) / (1.0 - x)) :
				-0.5 * log1p((-2.0 * x) / (1.0 + x));
	}
}

unittest
{
	assert(atanh(0) == 0);
	assert(feq(atanh(tanh(0.5)), 0.5));
}
+/

/*************************************
 * Hyperbolic arccotangent
 */

real acoth(real x)
{
	return 1 / acot(x);
}

unittest
{
	assert(feq(acoth(coth(0.01)), 100));
}

/*************************************
 * Hyperbolic arcsecant
 */

real asech(real x)
{
	return 1 / asec(x);
}

/*************************************
 * Hyperbolic arccosecant
 */

real acosech(real x)
{
	return 1 / acosec(x);
}

/*************************************
 * Convert string to float
 */

real atof(string s)
{
	if (!s.length)
		return real.nan;
	real result = 0;
	uint i = 0;
	while (s[i] == '\t' || s[i] == ' ')
		if (++i >= s.length)
			return real.nan;
	bool neg = false;
	if (s[i] == '-')
	{
		neg = true;
		i++;
	}
	else if (s[i] == '+')
		i++;
	if (i >= s.length)
		return real.nan;
	bool hex;
	if (s[s.length - 1] == 'h')
	{
		hex = true;
		s.length = s.length - 1;
	}
	else if (i + 1 < s.length && s[i] == '0' && s[i+1] == 'x')
	{
		hex = true;
		i += 2;
		if (i >= s.length)
			return real.nan;
	}
	else
		hex = false;
	while (s[i] != '.')
	{
		if (hex)
		{
			if ((s[i] == 'p' || s[i] == 'P'))
				break;
			result *= 0x10;
		}
		else
		{
			if ((s[i] == 'e' || s[i] == 'E'))
				break;
			result *= 10;
		}
		if (s[i] >= '0' && s[i] <= '9')
			result += s[i] - '0';
		else if (hex)
		{
			if (s[i] >= 'a' && s[i] <= 'f')
				result += s[i] - 'a' + 10;
			else if (s[i] >= 'A' && s[i] <= 'F')
				result += s[i] - 'A' + 10;
			else
				return real.nan;
		}
		else
			return real.nan;
		if (++i >= s.length)
			goto done;
	}
	if (s[i] == '.')
	{
		if (++i >= s.length)
			goto done;
		ulong k = 1;
		while (true)
		{
			if (hex)
			{
				if ((s[i] == 'p' || s[i] == 'P'))
					break;
				result *= 0x10;
			}
			else
			{
				if ((s[i] == 'e' || s[i] == 'E'))
					break;
				result *= 10;
			}
			k *= (hex ? 0x10 : 10);
			if (s[i] >= '0' && s[i] <= '9')
				result += s[i] - '0';
			else if (hex)
			{
				if (s[i] >= 'a' && s[i] <= 'f')
					result += s[i] - 'a' + 10;
				else if (s[i] >= 'A' && s[i] <= 'F')
					result += s[i] - 'A' + 10;
				else
					return real.nan;
			}
			else
				return real.nan;
			if (++i >= s.length)
			{
				result /= k;
				goto done;
			}
		}
		result /= k;
	}
	if (++i >= s.length)
		return real.nan;
	bool eneg = false;
	if (s[i] == '-')
	{
		eneg = true;
		i++;
	}
	else if (s[i] == '+')
		i++;
	if (i >= s.length)
		return real.nan;
	int e = 0;
	while (i < s.length)
	{
		e *= 10;
		if (s[i] >= '0' && s[i] <= '9')
			e += s[i] - '0';
		else
			return real.nan;
		i++;
	}
	if (eneg)
		e = -e;
	result *= pow(hex ? 2.0L : 10.0L, cast(real)e);
done:	
	return neg ? -result : result;
}

unittest
{
	assert(feq(atof("123"), 123));
	assert(feq(atof("+123"), +123));
	assert(feq(atof("-123"), -123));
	assert(feq(atof("123e2"), 12300));
	assert(feq(atof("123e+2"), 12300));
	assert(feq(atof("123e-2"), 1.23));
	assert(feq(atof("123."), 123));
	assert(feq(atof("123.E-2"), 1.23));
	assert(feq(atof(".456"), .456));
	assert(feq(atof("123.456"), 123.456));
	assert(feq(atof("1.23456E+2"), 123.456));
	//assert(feq(atof("1A2h"), 1A2h));
	//assert(feq(atof("1a2h"), 1a2h));
	assert(feq(atof("0x1A2"), 0x1A2));
	assert(feq(atof("0x1a2p2"), 0x1a2p2));
	assert(feq(atof("0x1a2p+2"), 0x1a2p+2));
	assert(feq(atof("0x1a2p-2"), 0x1a2p-2));
	assert(feq(atof("0x1A2.3B4"), 0x1A2.3B4p0));
	assert(feq(atof("0x1a2.3b4P2"), 0x1a2.3b4P2));
}