END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

kernel float4 quaternion_multiply(float4 p, float4 q)
{
    float w = p.w*q.w - dot(p.xyz,q.xyz);
    float3 v = float3(p.w,p.w,p.w)*q.xyz+
        float3(q.w,q.w,q.w)*p.xyz+
        cross(p.xyz,q.xyz);

    return float4(v.x,v.y,v.z,w);
}

kernel float2 bulb_iterate(float3 cval, int maxiter)
{
    float3 zval = cval;
    float4 qt2;
    float4 qt4;
    float4 qt8;

    float4 qp2;
    float4 qp4;
    float4 qp8;

    float4 quat_pre;
    float4 quat;
    float4 quat_post;
    float4 quat_half = float4(0.0f, -0.70710678f, 0.0f, 0.7071067811f);

    float2 res;
    float r2;
    float r;
    int i;

    r = sqrt(dot(zval,zval));

    for (i = 0; i < maxiter ; i++) {
        qt2.zw = normalize(zval.xy);
        qt2.x = 0.0f;
        qt2.y = 0.0f;
        qt4 = quaternion_multiply( qt2, qt2 );
        qt8 = quaternion_multiply( qt4, qt4 );

        qp2.x = 0.0f;
        qp2.z = 0.0f;
        qp2.w = zval.z / r;
        qp2.y = sqrt( 1.0f - qp2.w * qp2.w);
        qp4 = quaternion_multiply( qp2, qp2 );
        qp8 = quaternion_multiply( quat_half,
           quaternion_multiply( qp4, qp4 ));

        quat = float4(1.0f,0.0f,0.0f,0.0f);

        quat_pre = quaternion_multiply(qt8, qp8);
        quat_post.xyz = -quat_pre.xyz;
        quat_post.w = quat_pre.w;

        r = pow(r, 8.0f);
        quat = quaternion_multiply(
            quaternion_multiply( quat_pre, quat) , quat_post);

        zval = quat.xyz * float3(r,r,r) + cval;

        r2 = dot(zval,zval);
        r = sqrt(r2);

        if (r2>2.0f) break;

    }

    if (r<2.0f) {
        r = 0.0f;
    }

    if (i==0) {
        r = 0.0f;
    }

    res.x = (float)i;
    res.y = r;

    return res;
}