Since the rational polynome is a branch cut function, the result has no interest in Mb3d. I report the code and formula here for the interested people.
#include "utility_MATH.cpp" // include this always 1st - see my testing thread
void __attribute__((fastcall)) Formula(
double* x, // [eax]
double* y, // [edx]
double* arg, // [ebp+8], points to TIteration3Dext.C1
void* dummy // so we end w/ ret 8 as delphi expects
);
__attribute__((packed)) struct TIteration3Dext {
double Cw,Rold,RStopD,x,y,z,w; // use with neg indizes before C1. w is also used for 3d ifs and abox analytic DE
double Px, Py, Pz; // actual position, never change these! Can be used as input.
double Cx, Cy, Cz; //+24 these are the constants for adding. Pxyz or the julia seed. Cw is @-56 (begin of struct)
void* PVar; //+48 the actual formulas adress of constants and vars, constants at offset 0 increasing, user vars below -8
float SmoothItD; //+52
double Rout; //+56 the square of the vector length, can be used as input
int ItResultI; //+64
int maxIt; //+68
float RStop; //+72
int nHybrid[6]; //+76 all formulas iteration counts or single weights in interpol-hybrid
void* fHPVar[6]; //+100 all formulas pointer to constants+vars, PVars-8=0.5, use PVar for the actual formula
void* fHybrid[6]; //+124 the formulas adresses
int CalcSIT; //+148
int DoJulia; //+152
float LNRStop; //+156
int DEoption; //+160
float fHln[6]; //+164 for SmoothIts
int iRepeatFrom; //+188
double OTrap; //+192
double VaryScale; //+200 to use in vary by its
int bFirstIt; //+208 used also as iteration count, is set to 0 on it-start
int bTmp; //+212 tmpBuf, free of use.
double Dfree1; //+216
double Dfree2; //+224
double Deriv1; //+232 for 4D first deriv or as full derivs
double Deriv2; //+240
double Deriv3; //+248
float SMatrix4[4][4]; //+256 for 4d rotation, used like most other values only by the programs iteration loop procedure
};
// fastcall is not quite delphi fastcall.
// first two args are ok, third is in ecx in delphi, on stack here.
void __attribute__((fastcall)) Formula(
double* x, // [eax]
double* y, // [edx]
double* arg, // [ebp+8], points to TIteration3Dext.C1
void* dummy // so we end w/ ret 8 as delphi expects
) {
// Compute ptr to proper start of TIteration3Dext struct.
struct TIteration3Dext* cfg = (struct TIteration3Dext*)(arg-7);
double* cnst = (double*)(cfg->PVar);
int i;
double vx[8], vy[8], vz[8];
double OM = length(cfg->x,cfg->y,cfg->z);
double MO = recip(OM);
OM = intpow(OM,8);
vx[0] = MO*cfg->x; vy[0] = MO*cfg->y; vz[0] = MO*cfg->z;
for (i = 1; i < 8 ; i++)
{
vx[i] = vx[i-1] * vx[0];
vy[i] = vy[i-1] * vy[0];
vz[i] = vz[i-1] * vz[0];
}
cfg->x = OM*cnst[0]*(vx[0]+cnst[1]*vx[2]+cnst[2]*vx[4]-vx[6]) / (vx[7]+cnst[3]*vx[5]+cnst[4]*vx[3]+cnst[3]*vx[1]+load1()) + cfg->Cx;
cfg->y = OM*cnst[0]*(vy[0]+cnst[1]*vy[2]+cnst[2]*vy[4]-vy[6]) / (vy[7]+cnst[3]*vy[5]+cnst[4]*vy[3]+cnst[3]*vy[1]+load1()) + cfg->Cy;
cfg->z = OM*cnst[0]*(vz[0]+cnst[1]*vz[2]+cnst[2]*vz[4]-vz[6]) / (vz[7]+cnst[3]*vz[5]+cnst[4]*vz[3]+cnst[3]*vz[1]+load1()) + cfg->Cz;
}
edit;
an interesting variant, nothing exceptional though, can be got adding a fabs() to the denom so it never reach a negative value. I am adding this to the list...
example for x
cfg->x = OM*cnst[0]*(vx[0]+cnst[1]*vx[2]+cnst[2]*vx[4]-vx[6]) / (fabs(vx[7]+cnst[3]*vx[5]+cnst[4]*vx[3]+cnst[3]*vx[1])+load1()) + cfg->Cx;