Problem mit Fixed-Point-Arithmetik

*Q* 1

Hallo zusammen,

ich arbeite gerade an Funktionen für Fixed-Point-Arithmetik, um damit später das Mandelbrotfraktal genauer berechnen zu können als mit Floating-Point-Arithmetik.
Bisher habe ich schon addition und subtraktion zwischen 2 fp zahlen und umwandlung von und zu double.
Jetzt wage ich mich an die Multiplikation ran, aber irgendwo ist noch ein Fehler, sodass die ergebnisse falsch sind.
Ich habe mich ein wenig an diesem sourcecode hier orientiert: http://www.bealto.com/mp-mandelbrot-src.html
(Wichtige dateien MPBase.h und FPReal.h).
Allerdings will ich das so modifizieren, dass in den rechnungen NUR uint32 und int32 vorkommen. Der Hintergrund dafür ist, dass das ganze im fragment shader laufen soll. Ich schreibe das dann später mit glsl und soweit ich weiß gibt es da weder uint64 noch uint16, sondern nur int32 und uint32.
Wundert euch nicht, dass das so nach C aussieht und ohne operatorüberladung usw. ist, das ist auch weil ichs später nach glsl portieren will.

Also hier mal ein wenig code:

typedef unsigned int uint;
const int numMantissaUints = 7;
struct FixedPoint
{
	int sign; //0, -1 or 1
	uint mantissa[numMantissaUints];
};

// result = a * b
void mul(FixedPoint& result, const FixedPoint& a, const FixedPoint& b)
{
	result.sign = a.sign * b.sign;

	uint aux[2* numMantissaUints];
	memset(aux, 0, sizeof(aux));
	for(int x = 0; x < numMantissaUints; x++)
	{
		for(int y = 0; (x + y) < numMantissaUints; y++)
		{
			uint a1 = a.mantissa[x] & 0x0000FFFFu;
			uint a2 = (a.mantissa[x] & 0xFFFF0000u) >> 16;
			uint b1 = b.mantissa[y] & 0x0000FFFFu;
			uint b2 = (b.mantissa[y] & 0xFFFF0000u) >> 16;
			uint prodHinten = a1 * b1;
			uint prodVorne = a2 * b2;
			int idx = 2 * (x + y);
			aux[idx] += prodVorne;
			aux[idx + 1] += prodHinten;
		}
	}
	//write aux to result
	uint carry = 0;
	for(int i = numMantissaUints - 1; i >= 0; --i)
	{
		int idx = 2 * i;
		carry += aux[idx + 1];
		result.mantissa[i] = (carry & 0x0000FFFFu);
		carry >>= 16;
		carry += aux[idx];
		result.mantissa[i] |= ((carry & 0x0000FFFFu) << 16);
		carry >>= 16;
	}
}

Ich hoffe, dass das ganze halbwegs verständlich ist.
Weil ich nur uint32 hab, zerlege ich die uint32'er in 2 teile von je 16 bit und rechne damit.
Die Rechnung erfolgt quasi wie schriftliche multiplikationm.
Die Variable aux speichert die zwischenergebnisse und überträge und am ende wird alles zusammengefügt.
Allerdings ist wie gesagt da noch ein Fehler drin, sodass die Ergebnis etwas daneben liegen...
Also falls jemand vorschläge hat oder denm Fehler sieht, wäre das sehr gut

SeppJ

Darf ich dir dazu mal 2 Fragen stellen? (Oder 3 mit dieser hier)

1. Warum wird es genauer, wenn du fixed-point benutzt?
2. Wieso implementierst du das selber?

*Q* 1

1. Weil ich mehr bytes benutzen kann, ansonsten bin ich auf 4 oder mit double 8 bytes beschränkt. Hier kann ichs wählen wie ich will.

2. Weil das auf glsl laufen soll und ich deswegen nur wenige Datentypen benutzen kann. Ich habe ja eine Vorlage (auch verlinkt), aber diese benutzt uint64 und deswegen schreibe ich das ein wenig um.
Wenn du eine gute implementierung hast die mit int32 und uint32 auskommt, dann her damit bitte

Edit:
Hier mal zum Vergleich die Vorlage die mit 64 bit arbeitet.

// Z = X * Y (truncated to N words)
  friend void mul(FPReal & z,const FPReal & x,const FPReal & y)
  {
    if (x.sign == 0 || y.sign == 0) { z.zero(); return; }
    z.sign = x.sign * y.sign;

    // Multiply (the trivial way) and accumulate in AUX
    uint64 aux[N];
    memset(aux,0,N*sizeof(aux[0]));
    for (int i=0;i<N;i++) for (int j=0;j<N;j++)
    {
      int k = i+j;
      if (k > N) continue; // ignored
      uint64 u1 = (uint64)(x.m[i]) * (uint64)(y.m[j]);
      uint64 u0 = u1 & 0xFFFFFFFFULL; // lower 32 bits, index K
      u1 >>= (uint64)32; // higher 32 bits, index K-1
      if (k < N) aux[k] += u0;
      if (k > 0) aux[k-1] += u1;
    }
    // Propagate carry into result
    uint64 c = 0;
    for (int i=N-1;i>=0;i--)
    {
      c += aux[i];
      z.m[i] = (uint32)(c & 0xFFFFFFFFULL);
      c >>= (uint64)32;
    }
  }

*Q* 1

Ich glaub ich weiß jetzt wo das Problem liegt:

uint a1 = a.mantissa[x] & 0x0000FFFFu; 
uint a2 = a.mantissa[x] & 0xFFFF0000u; 
uint b1 = b.mantissa[y] & 0x0000FFFFu; 
uint b2 = b.mantissa[y] & 0xFFFF0000u;
uint prodHinten = a1 * b1;
uint prodVorne = a2 * b2;

Diese komponentenweise Multiplikation ist einfach grundlegend falsch glaube ich...
Das müsste noch irgendwie über kreuz gehen.

SeppJ

Q schrieb:

1. Weil ich mehr bytes benutzen kann, ansonsten bin ich auf 4 oder mit double 8 bytes beschränkt. Hier kann ichs wählen wie ich will.

Nein, das ist ein Argument, warum du nicht die eingebauten Datentypen benutzt, aber kein Argument für fixed-point. Es gibt haufenweise arbitrary precision floating point libraries.

2. Weil das auf glsl laufen soll und ich deswegen nur wenige Datentypen benutzen kann. Ich habe ja eine Vorlage (auch verlinkt), aber diese benutzt uint64 und deswegen schreibe ich das ein wenig um.
Wenn du eine gute implementierung hast die mit int32 und uint32 auskommt, dann her damit bitte

Die floating point libraries sollten größtenteils unabhängig von der unterliegenden Architektur sein.

*Q* 1

Ich wusste nicht das es noch alternativen zu fixed point gibt. Was würdest du denn empfehlen (ne bestimmte bibliothek oder so?)? Arbitrary precision hört sich für mich aber nach overkill an, weils ja unendlich genau ist...
Wichtig ist das es genauer als double ist (am besten einstellbar wie genau), trotzudem so schnell wie möglich und auf glsl benutzbar ist (heißt konkret nur uint32/int32).

Und ich hab nen video von jemandem gesehen der ins apfelmännchen reingezoomt ist und er hat das mit 832 (oder sowas) bit fixed point gemacht, deswegen dachte ich das wär ne gute wahl dafür.

Ich kann mir eig nicht vorstellen, dass arbitrary precision irgendwie besser sein kann. Das muss dann ja eigentlich mit dynamischen arrays gemacht werden und in glsl kann man soweit ich weiß nichtmal dynamischen speicher anfordern.
Und wenns nicht dynamisch ist dann ists ja eigentlich wieder fixed point außer wenn man noch einen exponenten speichert.

SeppJ

Q schrieb:

Ich wusste nicht das es noch alternativen zu fixed point gibt. Was würdest du denn empfehlen (ne bestimmte bibliothek oder so?)? Arbitrary precision hört sich für mich aber nach overkill an, weils ja unendlich genau ist...
Wichtig ist das es genauer als double ist (am besten einstellbar wie genau), trotzudem so schnell wie möglich und auf glsl benutzbar ist (heißt konkret nur uint32/int32).

Die basieren intern in der Regel auf Integertypen, da kannst du bestimmt auch wählen, was du möchtest. Und arbitrary heißt nicht unendlich. Es heißt, du kannst aussuchen, wie genau es sein soll.

Das muss dann ja eigentlich mit dynamischen arrays gemacht werden

Das könnte wiederum durchaus der Fall sein. Guck dir doch einfach mal eine gängige Implementierung an.

*Q* 1

Hab gerade mal bei GMP reingeguckt und ich muss sagen... ich versteh nur bahnhof.

Ich bin mit meiner implementierung ja auch fast fertig, mul ist das einzige, was noch nicht klappt. Ich hab jetzt mal versucht das mit dem über kreuz reinzumachen, aber auch das klappt nicht richtig.

Sieht jetzt so auss:

uint a1 = a.mantissa[x] & 0x0000FFFFu; 
uint a2 = (a.mantissa[x] & 0xFFFF0000u) >> 16;
uint b1 = b.mantissa[y] & 0x0000FFFFu;
uint b2 = (b.mantissa[y] & 0xFFFF0000u) >> 16;
uint cross = a1 * b2 + a2 * b1;
uint prodHinten = a1 * b1 + ((cross & 0x0000FFFF) << 16);
uint prodVorne = a2 * b2 + ((cross & 0xFFFF0000));
int idx = 2 * (x + y);
aux[idx] += prodVorne;
aux[idx + 1] += prodHinten;

Bei 36.84 * 0.15 kommt 5.40000759843457520000 raus statt 5.526 (nur als Beispiel).