Stimmt, da bin ich durcheinandergekommen.
Hab den Algorithmus jetzt mal neu aufgebaut und dabei noch einige "Fehler" im alten gefunden. U.a. hab ich zu wenig gecastet.
(ui16_t * ui16_t ist wieder ui16_t, obwohl z.B. eine ui32_t erwartet wird)

Jetzt funktionierts schon fast perfekt. Hier ist ein Video, in dem das Programm einfach gesagt nur folgendes durchführt:
- Hol die Rohwerte der Position.
- Berechne mit dem "mare_cirisum-Algorithmus" die echte Position
- Pinsel das Pixel an der entsprechenden Stelle weiss.

http://www.youtube.com/watch?v=fdpV7baRqgs

Wie man auf dem Video sieht, geht das in der linken Displayhälfte schon perfekt. Nur in der rechten ist noch ein kleiner Offset nach unten drin.
Ich vermute, dass ich da nochmal irgendwo nen Cast vergessen hab.
Lässt sich aber auch noch rausfinden.

Hier der neue Algorithmus:

Code:
// Those values are determined by a calibration procedure (touch all 4 corners)
// Here they are loaded with default values
ui16_t r0[2] = { 101, 145 };
ui16_t r1[2] = { 867, 150 };
ui16_t r2[2] = { 848, 802 };
ui16_t r3[2] = { 102, 759 };

// variables for precise algorithm; declared here to save stack memory
  ui16_t r_raw[2];    // raw value vector
  double W = 0;       // Iteration parameter; Start value = 0
  i16_t r1_r0[2];    // solution vector for r1 - r0
  i16_t r3_r0[2];    // solution vector for r3 - r0
  i16_t rraw_r0[2];  // solution vector for rraw - r0
  i32_t D;           // determinant
  i32_t D_part[2];   // partial solution for determinant
  double factor[2];  // solution vector factors
  i16_t r01_r32[2];  // solution vector for (r0 - r1) - (r3 - r2)
  double solVec[2];  // soluton vector
  double D_XY[2];    // determinant D_X, D_Y


#define LCD_Width   128
#define LCD_Height   64
#define Iterations   4
#define X 0
#define Y 1


void TPAD_GetPosPrecise(ui8_t *xp, ui8_t *yp)
{
  ui8_t step;

  // Check if constant 'Iteration' is defined already. If not, use default value of 5.
  #ifndef Iterations
    #define Iterations 5
  #endif

  TPAD_GetPos(&r_raw[X], &r_raw[Y]); // get raw values via 'old' method (ADC-Values)
//  r_raw[X] = 560; // test value
//  r_raw[Y] = 200; // test value
  // Constant values, needed for calculation; constant for a single iteration:
  
  // r1-r0
  r1_r0[X] = r1[X] - r0[X];
  r1_r0[Y] = r1[Y] - r0[Y];

  // r3-r0
  r3_r0[X] = r3[X] - r0[X];
  r3_r0[Y] = r3[Y] - r0[Y];

  // Determinant
  D_part[0] = (i32_t) r1_r0[X] * (i32_t) r3_r0[Y];
  D_part[1] = (i32_t) r1_r0[Y] * (i32_t) r3_r0[X];
  D = D_part[0] - D_part[1];

  // r_raw - r0
  rraw_r0[X] = r_raw[X] - r0[X];
  rraw_r0[Y] = r_raw[Y] - r0[Y];
  
  // (r0 - r1) - (r3 - r2)
  r01_r32[X] = (r0[X] - r1[X]) - (r3[X] - r2[X]);
  r01_r32[Y] = (r0[Y] - r1[Y]) - (r3[Y] - r2[Y]);

  // Iteration Step 0
  for (step = 0 ; step < Iterations ; step++)
  {  
  
    // Solution vector
    solVec[X] = (double) rraw_r0[X] - (W * (double) r01_r32[X]);
    solVec[Y] = (double) rraw_r0[Y] - (W * (double) r01_r32[Y]);

    // Determinant
    D_XY[X] = -(((double) r3_r0[X] * solVec[Y]) - ((double) r3_r0[Y] * solVec[X]));
    D_XY[Y] =   ((double) r1_r0[X] * solVec[Y]) - ((double) r1_r0[Y] * solVec[X]);

    // factor
    factor[X] = (double) D_XY[X] / (double) D;
    factor[Y] = (double) D_XY[Y] / (double) D;

    // New iteration parameter:
    W = factor[X] * factor[Y];
  }
  *xp = (ui16_t) (factor[X] * ((double) (LCD_Width - 1)));
  *yp = (ui16_t) (factor[Y] * ((double) (LCD_Height - 1)));
}