[PD-dev] question about q8_sqrt

[Christof Ressi]

Hi, I have a question about q8_sqrt (out of curiousity):

in d_math.c, the exported functions q8_sqrt and q8_rsqrt are defined like that:

t_float q8_rsqrt(t_float f0)
{
union {
float f;
long l;
} u;
u.f=f0;
if (u.f < 0) return (0);
else return (rsqrt_exptab[(u.l >> 23) & 0xff] *
rsqrt_mantissatab[(u.l >> 13) & 0x3ff]);
}

t_float q8_sqrt(t_float f0)
{
union {
float f;
long l;
} u;
u.f=f0;
if (u.f < 0) return (0);
else return (u.f * rsqrt_exptab[(u.l >> 23) & 0xff] *
rsqrt_mantissatab[(u.l >> 13) & 0x3ff]);
}

but the signal objects dont' simply use these functions but do this instead:

// sigrsqrt_perform
while (n--)
{
t_sample f = *in++;
union {
float f;
long l;
} u;
u.f = f;
if (f < 0) *out++ = 0;
else
{
t_sample g = rsqrt_exptab[(u.l >> 23) & 0xff] *
rsqrt_mantissatab[(u.l >> 13) & 0x3ff];
*out++ = 1.5 * g - 0.5 * g * g * g * f;
}
}

// sigsqrt_perform
while (n--)
    {
        t_sample f = *in++;
        union {
          float f;
          long l;
        } u;
        u.f = f;
        if (f < 0) *out++ = 0;
        else
        {
            t_sample g = rsqrt_exptab[(u.l >> 23) & 0xff] *
                rsqrt_mantissatab[(u.l >> 13) & 0x3ff];
            *out++ = f * (1.5 * g - 0.5 * g * g * g * f);
        }
    }

the last two versions strongly resemble the fast inverse sqrt from Quake: https://en.wikipedia.org/wiki/Fast_inverse_square_root.
my understanding is that the function version only give the first approximation whereas the perform routines apply Newton's method (once). please correct me if I'm wrong.

I guess the key difference is the first approximation: in Quake it's a magic number but in Pd there are two lookup tables. how do the lookup tables work and what's the advantage?

Christof