[sc-users] Sources of NaN in a UGen?

["Dan Stowell" <danstowell@xxxxxxxxx>]

Hi -

Writing a UGen to represent the "Chua's Circuit" chaotic oscillator, I
find that the UGen tends very quickly to outputting NaN, although I
can't tell why. The meat of the UGen only does multiplication and
addition (there are a couple of trivial divisions). Can anyone suggest
why this might be happening?

The UGen code is:

void ChuaL_next(ChuaL *unit, int inNumSamples)
{
       float *out = ZOUT(0);
       float freq = ZIN0(0);
       double a = ZIN0(1);
       double b = ZIN0(2);
       double c = ZIN0(3);
       double d = ZIN0(4);
       double rr = ZIN0(5); // Reciprocal of r - to save on a
division operation
       double h = ZIN0(6);
       double x0 = ZIN0(7);
       double y0 = ZIN0(8);
       double z0 = ZIN0(9);

       double xn = unit->xn;
       double yn = unit->yn;
       double zn = unit->zn;
       float counter = unit->counter;
       double xnm1 = unit->xnm1;
       double ynm1 = unit->ynm1;
       double znm1 = unit->znm1;
       double frac = unit->frac;

       float samplesPerCycle;
       double slope;
       if(freq < unit->mRate->mSampleRate){
               samplesPerCycle = unit->mRate->mSampleRate /
sc_max(freq, 0.001f);
               slope = 1.f / samplesPerCycle;
       }
       else {
               samplesPerCycle = 1.f;
               slope = 1.f;
       }

       if((unit->x0 != x0) || (unit->y0 != y0) || (unit->z0 != z0)){
               xnm1 = xn;
               ynm1 = yn;
               znm1 = zn;
               unit->x0 = xn = x0;
               unit->y0 = yn = y0;
               unit->z0 = zn = z0;
       }

       double dx = xn - xnm1;

       for (int i=0; i<inNumSamples; ++i) {
               if(counter >= samplesPerCycle){
                       counter -= samplesPerCycle;
                       frac = 0.f;

                       xnm1 = xn;
                       ynm1 = yn;
                       znm1 = zn;

                       double k1x, k2x, k3x, k4x,
                               k1y, k2y, k3y, k4y,
                               k1z, k2z, k3z, k4z,
                               kxHalf, kyHalf, kzHalf,
                               ymx, ymx3, ymx3b, hOverR, xnm1plus;

                       hOverR = h * rr;

                       // 4th order Runge-Kutta approximate solution
for differential equations

                       ymx = ynm1 - xnm1;
                       ymx3 = ymx * ymx * ymx;
                       ymx3b = ymx3 * b;
                       k1x = hOverR * (d * ymx - c * xnm1 - (a * xnm1
* xnm1 * xnm1) - ymx3b);
                       k1y = - (h * (d * ymx - znm1 - ymx3b));
                       k1z = h * ynm1;
                       kxHalf = k1x * 0.5;
                       kyHalf = k1y * 0.5;
                       kzHalf = k1z * 0.5;

                       ymx = ynm1 + kyHalf - xnm1 - kxHalf;
                       ymx3 = ymx * ymx * ymx;
                       ymx3b = ymx3 * b;
                       xnm1plus = xnm1 + kxHalf;
                       k2x = hOverR * (d * ymx - c * (xnm1 + kxHalf)
- (a * xnm1plus *
xnm1plus * xnm1plus) - ymx3b);
                       k2y = - (h * (d * ymx - znm1 - kzHalf - ymx3b));
                       k2z = h * (ynm1 + kyHalf);
                       kxHalf = k2x * 0.5;
                       kyHalf = k2y * 0.5;
                       kzHalf = k2z * 0.5;

                       ymx = ynm1 + kyHalf - xnm1 - kxHalf;
                       ymx3 = ymx * ymx * ymx;
                       ymx3b = ymx3 * b;
                       xnm1plus = xnm1 + kxHalf;
                       k3x = hOverR * (d * ymx - c * (xnm1 + kxHalf)
- (a * xnm1plus *
xnm1plus * xnm1plus) - ymx3b);
                       k3y = - (h * (d * ymx - znm1 - kzHalf - ymx3b));
                       k3z = h * (ynm1 + kyHalf);

                       xnm1plus = xnm1 + k3x;
                       k4x = hOverR * (d * ymx - c * (xnm1 + k3x) -
(a * xnm1plus *
xnm1plus * xnm1plus) - ymx3b);
                       k4y = - (h * (d * ymx - znm1 - k3z - ymx3b));
                       k4z = h * (ynm1 + k3y);

                       xn = xn + (k1x + 2.0*(k2x + k3x) + k4x) * ONESIXTH;
                       yn = yn + (k1y + 2.0*(k2y + k3y) + k4y) * ONESIXTH;
                       zn = zn + (k1z + 2.0*(k2z + k3z) + k4z) * ONESIXTH;

                       dx = xn - xnm1;
               }
               counter++;
               ZXP(out) = (xnm1 + dx * frac) * 0.5f;
               frac += slope;
       }

       unit->xn = xn;
       unit->yn = yn;
       unit->zn = zn;
       unit->counter = counter;
       unit->xnm1 = xnm1;
       unit->ynm1 = ynm1;
       unit->znm1 = znm1;
       unit->frac = frac;
}

void ChuaL_Ctor(ChuaL* unit){
       SETCALC(ChuaL_next);

       unit->x0 = unit->xn = unit->xnm1 = ZIN0(7);
       unit->y0 = unit->yn = unit->ynm1 = ZIN0(8);
       unit->z0 = unit->zn = unit->znm1 = ZIN0(9);
       unit->counter = 0.f;
       unit->frac = 0.f;

       ChuaL_next(unit, 1);
}


and the SC class:

ChuaL : ChaosGen {
       *ar { arg freq=22050, a=0.3286, b=0.9336, c= -0.8126, d=0.399,
rr=0.6.reciprocal, h=0.05, xi=0.1, yi=0, zi=0, mul=1.0, add=0.0;
               ^this.multiNew('audio', freq, a, b, c, d, rr, h, xi,
yi, zi).madd(mul, add)
       }
}

Thanks in advance for any advice -
Dan

P.S. Not sure if this is an "sc-users" or "sc-dev" issue so sorry for
sending to both lists, for those who read both.

--
http://www.mcld.co.uk