...für einen Speziellen Typ (int, unsigned int, short, usw) und einer gewissan Anzahl
an Komma-Bits, möge man die Multiplikation und division selbständig für die Fixkomma-
Zahlen adaptieren. 
Ich habe dafür eine komplette Template-Klasse in C++ geschrieben, die mir alles übernimmt,
auch das statuische multiplizieren mit floats  usw. Alles fixkomma.

Das ist hochoptimeirt das Heron_Verfahren 2x ineinander eingesetzt, sodass man nur 1x dividieren muss,
da dividieren i.d.R. langsam ist.
Prinzipiell ist die Wurzel genau genug, um damit einen Software-3D-Renderer zu implementieren.
Dzurch die Optimierung kommen gewisse Magic-Numbers zu stande, die nunmal so sind, wie sie sind.
Getestet und Funktionsfähig. 
Ansonsten kann mans auch gerne selbst besser machen...

	FixedPoint<Typ, FracBits> fast_rsqrt() const
	{
		if (Data <= 0) return 0;
		int ld = (31-__builtin_clz(Data)) - FracBits;
		if (ld >= 0)
		{
			if (ld >= ((31-FracBits)/2)-1)
			{
				unsigned long long sq = (((unsigned long long)Data)*Data)>>FracBits;
				unsigned int a = 6U << ld ;
				unsigned long long b = 1ULL << ((ld<<1)+FracBits);
				unsigned int c = 1U << (ld>>1);
				unsigned long long d = 1ULL << (((ld*6)>>2)+FracBits);
				if (ld&1)
				{
					c *= 23; c >>= 4;
					d *= 22; d >>= 4;
				}
				FixedPoint<Typ, FracBits> x;
				x.Data =  div_64_64u(mac_u(c, Data, d)<<(FracBits+2), mac_u(a, Data, sq+b));
				return x;
			}
			else
			{
				FixedPoint<Typ, FracBits> sq = (*this)*(*this);
				unsigned int a = 6U << ld;
				unsigned int b = 1U << (ld<<1);
				FixedPoint<Typ, FracBits> c = (1U << (ld>>1));
				FixedPoint<Typ, FracBits> d = (1U << ((ld*6)>>2));
				if (ld&1)
				{
					c *= 1.4142136;
					d *= 1.4142136;
				}
				FixedPoint<Typ, FracBits> x = (FixedPoint<Typ, FracBits>(c * (*this)) + d) / (( sq + a*(*this) + b)>>2);
				return x;
			}
		}
		else
		{		
			ld = Abs(ld);
			FixedPoint<Typ, FracBits> sq = (*this)*(*this);
			FixedPoint<Typ, FracBits> a = 6; a >>= ld;
			FixedPoint<Typ, FracBits> b = 1;  b >>= ld<<1;
			FixedPoint<Typ, FracBits> c = 1;  c >>= ld>>1;
			FixedPoint<Typ, FracBits> d = 1;  d >>= (ld*6)>>2;
			if (ld&1)
			{
				c *= 0.70710578;
				d *= 0.70710578;
			}
			FixedPoint<Typ, FracBits> x =  (FixedPoint<Typ, FracBits>(c * (*this)) + d) / (( sq + a*(*this) + b)>>2);
			return x;
		}
	}