float t = phase / period; //normalized position in period
float dt = deltaPhase / period; // normalized phasestep

float polyBLAMP = 0.0f;
if (t < dt) {
    float d = (t)/dt; // fractional delay
    polyBLAMP = (3.0f*powf(d,5.0f)-10.0f*powf(d,4.0f)+40.0f*powf(d,2.0f)-60.0f*d/2.0f+28.0f)/120.0f;
} else if (t < dt*2.0f) {
    float d = (t-dt)/dt;
    polyBLAMP = (-powf(d,5.0f)+5.0f*powf(d,4.0f)-10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)-5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt) {
    float d = (dt-1.0f+t)/dt;
    polyBLAMP = (-3.0f*powf(d,5.0f)+5.0f*powf(d,4.0f)+10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)+5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt*2.0f) {
    float d = (2.0f*dt-1.0f+t)/dt;
    polyBLAMP = powf(d,5.0f)/120.0f;
}
float output = osc(phase)+polyBLAMP;

public static double blep(double p, double dp) {
    double y = 0;
    if (0 <= p && p < 2 * dp) {
        final double x = p / dp;
        double u = 2 - x;
        u *= u;
        u *= u;
        y += u;
        if (p < dp) {
            double v = 1 - x;
            v *= v;
            v *= v;
            y -= 4 * v;
        }
    }
    p = 1 - p;
    if (0 <= p && p < 2 * dp) {
        final double x = p / dp;
        double u = 2 - x;
        u *= u;
        u *= u;
        y -= u;
        if (p < dp) {
            double v = 1 - x;
            v *= v;
            v *= v;
            y += 4 * v;
        }
    }
    return y / 12;
}

public static double saw(double p, double dp) {
    double v = 2 * p - 1;
    v += blep(p, dp);
    return v;
}

double y = saw(phi, dphi);
phi += d;
phi -= Math.floor(phi);
y -= saw(phi, dphi);