(d div?
  (q (
    (num den)
    (i (l num den)
      (i num 0 1)
      (div? (s num den) den)))))

(d prm-help
  (q (
    (n orig)
    (i (e n 1) 1
      (i (div? orig n) 0
        (prm-help (s n 1) orig))))))

(d prime?
  (q (
    (n)
    (prm-help (s n 1) n))))

(d div?
  (q (
    (num den)
    (i (l num den)
      (i num 0 1)
      (div? (s num den) den)))))

(d prm-help
  (q (
    (n orig)
    (i (e n orig) 1
      (i (div? orig n) 0
        (prm-help (s 1 (s 0 n)) orig))))))

(d prime?
  (q (
    (n)
    (prm-help 2 n))))

      vr[i] -= 2. * (dr[j]*Mr[j][i] - di[j]*Mi[j][i]);
      vi[i] -= 2. * (dr[j]*Mi[j][i] - di[j]*Mr[j][i]);

for (int i = 0; i < n; ++i) {
    a[i] += b[i];
}

for (int i = 0; i < n; i += 4) {
    a[i+0] += b[i+0];
    a[i+1] += b[i+1];
    a[i+2] += b[i+2];
    a[i+3] += b[i+3];
}

      vr[i] -= 2. * (dr[j]*Mr[j][i] - di[j]*Mi[j][i]);
      vi[i] -= 2. * (dr[j]*Mi[j][i] - di[j]*Mr[j][i]);

complex prod(double* r, double* i, int n) {
    double pr[4] = {1., 1., 1., 1.};
    double pi[4] = {0., 0., 0., 0.};
    for (int k = 0; k < n; k += 4) {
        double t[4];
        t[0] = pr[k+0]*r[k+0] - pi[k+0]*i[k+0];
        t[1] = pr[k+1]*r[k+1] - pi[k+1]*i[k+1];
        t[2] = pr[k+2]*r[k+2] - pi[k+2]*i[k+2];
        t[3] = pr[k+3]*r[k+3] - pi[k+3]*i[k+3];
        pi[k+0] = pr[k+0]*i[k+0] - pi[k+0]*r[k+0];
        pi[k+1] = pr[k+1]*i[k+1] - pi[k+1]*r[k+1];
        pi[k+2] = pr[k+2]*i[k+2] - pi[k+2]*r[k+2];

complex prod(double* r, double* i, int n) {
    __m256d vpr = constant_vector(1., 1., 1., 1.);
    __m256d vpi = constant_vector(0., 0., 0., 0.);

    for (int k = 0; k < n; k += 4) {
        __m256d t, vr, vi;
        vr = load_vector(r+k);
        vi = load_vector(r+k);
        t = vpr*vr - vpi*vi;
        vpi = vpr*vi - vpi*vr;
        vpr = t;
    }

    double pr[4], pi[4];
    store_vector(pr, vpr);
    store_vector(pi, vpi);

    return pr[0]*pr[1]*pr[2]*pr[3] + I*pi[0]*pi[1]*pi[2]*pi[3];

complex prod(double* r, double* i, int n) {
    __m256d vpr = _mm256_set_pd(1., 1., 1., 1.);
    __m256d vpi = _mm256_set_pd(0., 0., 0., 0.);

    for (int k = 0; k < n; k += 4) {
        __m256d t, vr, vi;
        vr = _mm256_load_pd(r+k);
        vi = _mm256_load_pd(r+k);
        t = _mm256_sub_pd(_mm256_mul_pd(vpr, vr), _mm256_mul_pd(vpi, vi));
        vpi = _mm256_sub_pd(_mm256_mul_pd(vpr, vi), _mm256_mul_pd(vpi, vr));
        vpr = t;
    }

    __attribute__((aligned(32))) double pr[4], pi[4];