/* Last edited on 2025-02-20 06:44:29 by stolfi */
    
    auto void fill_points(int32_t ip0, int32_t np, double A, double B, double C, double D);
      /* Fills the data points {vP[ip],wP[ip]}, for {ip = ip0+k} and {k}
        in {0..np-1}. In most cases, the weights {wP[ip]} will be
        randomly chosen in {[0_1]}, and the value {vP[ip]} will be taken
        from the polynomial {A*z^3 + B*z^2 + C*z + D} where {z} ranges
        from 0 to 1 as {k} ranges from 0 to {np-1}. Occasionally there
        will be gaps of one or more points. where {wP[ip]} is 0 and
        {vP[ip]} is {NAN}.
        
        Also fills {vH[ip]}, for {ip = ip0+k} and {k} in {0..np},
        with the same polynomial evaluated at half a step before {vP[ip]}. */
     
    int32_t NL = N/5;          /* Num of points in linear part. */
    int32_t NQ = N/4;          /* Num of points in quadratic part. */
    int32_t NK = N - NL - NQ;  /* Num of points in cubic part. */ 
    
    fill_points(0, NL,     00.0, 00.0, +1.0, +2.0);
    fill_points(NL, NQ,    00.0, -1.0, +3.0, +2.0);
    fill_points(NL+NQ, NK, +4.0, -6.0, +3.0, +2.0);

    
    void fill_points(int32_t ip0, int32_t np, double A, double B, double C, double D)
      { int32_t nskip = 0; /* Points to skip. */
        int32_t ngaps = 0; /* Gaps generated so far. */
        for (int32_t k = 0; k < np; k++)
          { double zPk = ((double)k)/((double)np -1);
            double vPk, wPk;
            if ((k <= 1) || (k >= np-2) || (nskip > 0))
              { vPk = NAN; wPk = 0;
                if (nskip > 0) { nskip--; }
              }
            else
              { vPk = ((A*zPk + B)*zPk + C)*zPk + D;
                wPk = dabrandom(0, 1);
              }
            /* Just in case: */
            if (wPk == 0) { vPk = NAN; }
            /* Save point: */
            int32_t ip = ip0 + k;
            vP[ip] = vPk; wP[ip] = wPk;
            /* Define random bursts of faults: */
            if ((nskip == 0) && (drandom() < 0.10))
              { nskip = ngaps+1;
                ngaps++;
              }
          }

        for (int32_t k = 0; k <= np; k++)
          { double zHk = ((double)k-0.5)/((double)np -1);
            double vHk = ((A*zHk + B)*zHk + C)*zHk + D;
            /* Save point: */
            int32_t ip = ip0 + k;
            vH[ip] = vHk;
          }
      }