// Interesting...
// In VC++ 5.0, if you statically link with the multi-threaded C run
// time, your executable (well, this one) gets 50K bigger if you
// include <iostream> instead of <iostream.h>. Stupid!
// I'll be stdio.h would be even smaller, but I'm not going to check.
#include <iostream.h>
#include <cassert>
#include <cmath>

namespace std {}
using namespace std;

/*
	Written by:

Bruce Dawson, Cygnus Software      Author of Fractal eXtreme for Win32
http://www.cygnus-software.com/           comments@cygnus-software.com

Last modified, March 2001


	This source code is supplied on an as-is basis.  You are free
to use it and modify it, as long as this notice remains intact, and
as long as Cygnus Software receives credit for any use where it is
a substantial portion of the information involved.

	This source code has little in common with the high precision
math code used for Fractal eXtreme, which has been optimized much
more.


	This source file declares and defines a partially functional
high precision math class, then uses it to calculate PI to high
precision.


For more information on the algorithms see:


http://www.cygnus-software.com/misc/pidigits.htm


	In a real application, the class declaration, definition,
and usage would be in three separate files.

	This code was compiled with MSVC5.0.  It should compile on
MSVC6.0, and many other compilers.

	The platform specific code has been isolated at the top of
the file, to allow adjusting to different sizes of ints, and
little endian versus big endian byte ordering.

	This is not a complete math class.  In particular, it does
not have full support for negative numbers.

Exercises:

	1) Fully implement this class, to properly handle negative numbers,
multiplication and division of InfPrec numbers, cos, sin, etc. of
InfPrec numbers.

	2) Optimize this class.

	3) Template this class on its precision.
*/

#if	(_MSC_VER < 1100)
	#pragma warning(disable : 4237)	// Shutup warnings about bool.
	typedef	int bool;

	#define	true	(bool)1
	#define	false	(bool)0
#endif

// Platform specific!
// Are we little endian?  Comment this out if we are not.
#define	LITTLE_ENDIAN

// Platform specific!
#define	BITSPERLONG	32

// uint32 should be twice as big as uint16.
// Platform specific!
typedef	unsigned int uint32;
// Platform specific!
typedef	unsigned short uint16;

#define	LOWORD(l)		((uint16)(l))
#define HIWORD(l)		((uint16)(((uint32)(l) >> 16) & 0xFFFF))
#define	MAKELONG(h,l)	((uint32)((h << 16) | l))



// For simplicity, statically determine the accuracy.
// Do a small enough number of digits that it fits on
// a 25x80 command prompt.
#define	NUMLONGS	190
#define	NUMWORDS	(NUMLONGS*2)

class InfPrec
{
public:
	InfPrec(int Value = 0) {Init(Value);}
	InfPrec(const InfPrec& rhs) {Init(rhs);}
	InfPrec& operator=(const InfPrec& rhs) {Init(rhs); return *this;}

	bool operator==(int rhs) const;
	bool operator<(int rhs) const;
	bool operator!=(int rhs) const {return !(*this == rhs);}
	bool operator<=(int rhs) const {return *this < rhs || *this == rhs;}
	bool operator>=(int rhs) const {return !(*this < rhs);}
	bool operator>(int rhs) const {return !(*this <= rhs);}

	InfPrec& operator*=(int rhs);
	InfPrec operator*(int rhs) const {return InfPrec(*this) *= rhs;}

	InfPrec& operator/=(int rhs);
	InfPrec operator/(int rhs) const {return InfPrec(*this) /= rhs;}

	InfPrec& operator+=(const InfPrec& rhs);
	InfPrec operator+(const InfPrec& rhs) const {return InfPrec(*this) += rhs;}

	InfPrec operator-() const;
	InfPrec& operator-=(const InfPrec& rhs);
	InfPrec operator-(const InfPrec& rhs) const {return InfPrec(*this) -= rhs;}

	friend ostream& operator<<(ostream& output, const InfPrec& rhs);

private:
	void Init(int Value);
	void Init(const InfPrec& rhs);

	// These accessors make it easy to reverse the endianness
	// of the number representation.  The current implementation
	// assumes Intel style format, with the most significant
	// word coming at the end.
	// Return a particular 'digit' of the number.  GetLong(0)
	// returns the most significant long word.
	#ifdef	LITTLE_ENDIAN
	uint32& GetLong(int Index) {return mNumber[NUMLONGS - 1 - Index];}
	const uint32& GetLong(int Index) const {return mNumber[NUMLONGS - 1 - Index];}
	uint16& GetWord(int Index) {return reinterpret_cast<uint16 *>(mNumber)[NUMWORDS - 1 - Index];}
	const uint16& GetWord(int Index) const {return reinterpret_cast<const uint16 *>(mNumber)[NUMWORDS - 1 - Index];}
	#else
	uint32& GetLong(int Index) {return mNumber[Index];}
	const uint32& GetLong(int Index) const {return mNumber[Index];}
	uint16& GetWord(int Index) {return reinterpret_cast<uint16 *>(mNumber)[Index];}
	const uint16& GetWord(int Index) const {return reinterpret_cast<const uint16 *>(mNumber)[Index];}
	#endif

	uint32	mNumber[NUMLONGS];
};

bool InfPrec::operator==(int rhs) const
{
	if ((int)GetLong(0) != rhs)
		return false;
	for (int i = 1; i < NUMLONGS; i++)
		if (GetLong(i) != 0)
			return false;
	return true;
}

bool InfPrec::operator<(int rhs) const
{
	// Surprisingly, this all you need to check.
	// The least significant bits cannot affect
	// the result.
	return (int)GetLong(0) < rhs;
}

InfPrec& InfPrec::operator*=(int rhs)
{
	assert(rhs >= 0);
	assert(*this >= 0);
	uint32 accum = 0;
	for (int i = NUMWORDS - 1; i >= 0; i--)
	{
		accum += GetWord(i) * rhs;
		GetWord(i) = LOWORD(accum);
		accum = HIWORD(accum);
	}
	return *this;
}

InfPrec& InfPrec::operator/=(int rhs)
{
	assert(rhs > 0);
	assert(*this >= 0);
	uint16 remainder = 0;
	for (int i = 0; i < NUMWORDS; i++)
	{
		uint32 divisor = MAKELONG(remainder, GetWord(i));
		GetWord(i) = (uint16)(divisor / rhs);
		remainder = (uint16)(divisor % rhs);
	}
	return *this;
}

InfPrec& InfPrec::operator+=(const InfPrec& rhs)
{
	uint32 accum = 0;
	for (int i = NUMWORDS - 1; i >= 0; i--)
	{
		accum += GetWord(i) + rhs.GetWord(i);
		GetWord(i) = LOWORD(accum);
		accum = HIWORD(accum);
	}
	return *this;
}

InfPrec InfPrec::operator-() const
{
	InfPrec Result = *this;
	int i;
	// Negating a number is done by doing a bitwise NOT
	// then adding one at the least significant position.
	for (i = 0; i < NUMLONGS; i++)
		Result.GetLong(i) = ~Result.GetLong(i);
	bool carry = true;
	for (i = NUMLONGS - 1; i >= 0 && carry; i--)
		if (++Result.GetLong(i) != 0)
			carry = false;
	return Result;
}

InfPrec& InfPrec::operator-=(const InfPrec& rhs)
{
	*this += -rhs;
	return *this;
}

void InfPrec::Init(int Value)
{
	GetLong(0) = Value;
	for (int i = 1; i < NUMLONGS; i++)
		GetLong(i) = 0;
}

void InfPrec::Init(const InfPrec& rhs)
{
	for (int i = 0; i < NUMLONGS; i++)
		mNumber[i] = rhs.mNumber[i];
}

ostream& operator<<(ostream& output, const InfPrec& rhs)
{
	InfPrec Copy = rhs;
	output << Copy.GetLong(0) << '.';
	Copy.GetLong(0) = 0;
	int NumDigits = int((NUMLONGS - 1) * BITSPERLONG * (log(2.0) / log(10.0)));
	while (Copy != 0 && NumDigits-- > 0)
	{
		Copy *= 10;
		output << Copy.GetLong(0);
		Copy.GetLong(0) = 0;
	}
	return output;
}

// Calculate the arc tangent of 1/x

InfPrec ataninvint(int x)
{
    InfPrec Result = InfPrec(1) / x;
    int XSquared = x * x;

    int Divisor = 1;
    InfPrec Term = Result;
    InfPrec TempTerm;
    while (Term != 0)
    {
        Divisor += 2;
        Term /= XSquared;
        Result -= Term / Divisor;
    
        Divisor += 2;
        Term /= XSquared;
        Result += Term / Divisor;
    }
    return Result;
}

int main(int argc, char *argv[])
{
	InfPrec PI = (ataninvint(5) * 4 - ataninvint(239)) * 4;

	cout << "PI = " << PI << endl;

	cout << "Press enter to continue...";
	char temp;
	cin.read(&temp, 1);
	return 0;
}