Arbitrary precision decimal numbers

[Dom Disc]

On Friday, 5 August 2022 at 14:25:39 UTC, Ruby The Roobster wrote:
> I'm currently working on an arbitrarily precise division 
> algortihm based off of the done-by-hand standard algorithm, but 
> I need to get it to work for complex numbers, [which it's just 
> not](https://forum.dlang.org/post/czidpbdywsohstyviitk@forum.dlang.org).

I once did a completely inline implementation of xlcmplx based on 
an arbitrary precision float (in this context called xlfloat) in 
C++.
You only need to exchange xlfloat with your floatingpoint type:

```C
class xlcmplx
{
	xlfloat Re, Im;
public:
	friend inline void swap(xlcmplx& a, xlcmplx& b) { swap(a.Re, 
b.Re); swap(a.Im, b.Im); }

	xlcmplx(const xlfloat& re =0, const xlfloat& im =0) : Re(re), 
Im(im) { }
	xlcmplx(const xlcmplx& c) : Re(c.Re), Im(c.Im) { }
	// two special "constructors":
	friend inline xlcmplx i(const xlfloat& x =1) { return 
xlcmplx(0,x); } // make pure imaginary number
	friend inline xlcmplx polar(const xlfloat& arg, const xlfloat& 
norm =1) { return (!norm || !arg) ? norm : 
xlcmplx(cos(arg),sin(arg)) *= norm; }

	uint32 Decimal(char* s, uint32 max) const; // string 
representation

	inline xlcmplx& operator=(const xlcmplx& c) { if(this != &c) Re 
= c.Re, Im = c.Im; return *this; }
	inline xlcmplx& operator=(const xlfloat& f) { Re = f, Im = 0; 
return *this; }

	inline xlcmplx& operator+() { return *this; } // dummy operator
	inline xlcmplx operator-() const { return xlcmplx(-Re, -Im); }
	inline xlcmplx& operator+=(const xlcmplx& c) { Re += c.Re; Im += 
c.Im; return *this; }
	inline xlcmplx& operator+=(const xlfloat& f) { Re += f; return 
*this; }
	inline xlcmplx& operator-=(const xlcmplx& c) { Re -= c.Re; Im -= 
c.Im; return *this; }
	inline xlcmplx& operator-=(const xlfloat& f) { Re -= f; return 
*this; }
	inline xlcmplx& operator*=(const xlcmplx& c) { xlfloat t = Re; t 
*= c.Re; t -= Im*c.Im; Im *= c.Re; Im += Re*c.Im; Re = t; return 
*this; }
	inline xlcmplx& operator*=(const xlfloat& f) { Re *= f; Im *= f; 
return *this; }
	inline xlcmplx& operator/=(const xlcmplx& c) { return *this *= 
inv(c); }
	inline xlcmplx& operator/=(const xlfloat& f) { Re /= f; Im /= f; 
return *this; }

	inline bool operator==(const xlcmplx& c) const { return Im == 
c.Im && Re == c.Re; }
	inline bool operator==(const xlfloat& f) const { return !Im && 
Re == f; }
	inline bool operator!=(const xlcmplx& c) const { return Im != 
c.Im || Re != c.Re; }
	inline bool operator!=(const xlfloat& f) const { return !!Im || 
Re != f; }

	friend inline const xlfloat& real(const xlcmplx& c) { return 
c.Re; }
	friend inline const xlfloat& imag(const xlcmplx& c) { return 
c.Im; }
	friend inline xlfloat norm(const xlcmplx& c) { return c.Im ? 
c.Re ? sqr(quad(c.Re)+quad(c.Im)) : abs(c.Im) : abs(c.Re); }
	friend inline xlfloat arg(const xlcmplx& c) { return atan2(c.Re, 
c.Im); } // an angle in [0..2*pi[

	friend inline xlcmplx conj(xlcmplx c) { c.Im.neg(); return c; }
	friend inline xlcmplx inv(const xlcmplx& c) { xlfloat t = 
quad(c.Re); t += quad(c.Im); return conj(c) /= t; }
	friend inline xlcmplx quad(const xlcmplx& c) { return 
xlcmplx(quad(c.Re)-quad(c.Im),(c.Re*c.Im)<<1); }
	friend inline xlcmplx exp(const xlcmplx& c) { return c.Im ? 
xlcmplx(cos(c.Im), sin(c.Im)) *= exp(c.Re) : exp(c.Re); }
	friend inline xlcmplx sin(const xlcmplx& c) { return c.Im ? c.Re 
? xlcmplx(sin(c.Re)*cosh(c.Im), cos(c.Re)*sinh(c.Im)) : 
i(sinh(c.Im)) : sin(c.Re); }
	friend inline xlcmplx cos(const xlcmplx& c) { return c.Im ? c.Re 
? xlcmplx(cos(c.Re)*cosh(c.Im), sin(c.Re)*sinh(c.Im)) : 
cosh(c.Im) : cos(c.Re); }
	friend inline xlcmplx ln(const xlcmplx& c) { return 
xlcmplx(ln(norm(c)), arg(c)); } // main value - multiples of 
2*pi*i added are also valid results
};

inline bool operator==(const xlfloat& f, const xlcmplx& c) { 
return c == f; }
inline bool operator!=(const xlfloat& f, const xlcmplx& c) { 
return c != f; }

inline xlcmplx operator+(xlcmplx a, const xlcmplx& b) { return a 
+= b; }
inline xlcmplx operator+(xlcmplx a, const xlfloat& b) { return a 
+= b; }
inline xlcmplx operator+(const xlfloat& a, xlcmplx b) { return b 
+= a; }
inline xlcmplx operator-(xlcmplx a, const xlcmplx& b) { return a 
-= b; }
inline xlcmplx operator-(xlcmplx a, const xlfloat& b) { return a 
-= b; }
inline xlcmplx operator-(const xlfloat& a, const xlcmplx& b) { 
return -b += a; }
inline xlcmplx operator*(xlcmplx a, const xlfloat& b) { return a 
*= b; }
inline xlcmplx operator*(xlcmplx a, const xlcmplx& b) { return a 
*= b; }
inline xlcmplx operator*(const xlfloat& a, xlcmplx b) { return b 
*= a; }
inline xlcmplx operator/(xlcmplx a, const xlcmplx& b) { return a 
/= b; }
inline xlcmplx operator/(xlcmplx a, const xlfloat& b) { return a 
/= b; }
inline xlcmplx operator/(const xlfloat& a, const xlcmplx& b) { 
return xlcmplx(a) /= b; }

inline xlcmplx pow(const xlcmplx& c, const int exp) { return 
polar(arg(c)*exp, pow(norm(c), exp)); }
inline xlcmplx pow(const xlcmplx& c, const xlfloat& exp) { return 
polar(arg(c)*exp, pow(norm(c), exp)); }
inline xlcmplx log(const xlcmplx& c, const uint32 base = 10) { 
return ln(c) /= ln(number(base)); }
```

You can also get the implementation of xlfloat, if you like - but 
that one is NOT completely inline and much longer :-D