Jacobi theta functions

[stla]

Hello,

I've implemented the Jacobi theta functions in Asymptote. This seems to work fine. Are you interested in including this implementation in Asymptote, or in a package? (I don't know how to do an Asymptote package). The Jacobi theta functions are the basis of numerous other special functions, e.g. the Dedekind eta function and the Neville theta functions.

real modulo(real a, real p) {
    real i = a > 0 ? floor(a / p) : ceil(a / p);
    return a - i * p;
}

pair _i_ = (0.0, 1.0);

pair calctheta3(pair z, pair tau) {
    pair out = (1.0, 0.0);
    int n = 0;
    while(true) {
        n = n + 1;
        real nn = n;
        pair qweight = exp(nn * _i_ * pi * (nn * tau + 2.0 * z)) +
            exp(nn * _i_ * pi * (nn * tau - 2.0 * z));
        out += qweight;
//      if(length(out) == 0) {
//          error("log(0)");
//      } else 
        if(n >= 3 && (out + qweight) == out) {
            break;
        }
    }
    return log(out);
}

pair argtheta3(pair z, pair tau, int pass_in) {
    int passes = pass_in + 1;
//  if(passes > 800) {
//      error("passes > 800 (argtheta3)");
//  }
    real z_img = z.y;
    real h = tau.y / 2.0;
    pair zuse = (modulo(z.x, 1.0), z_img);
    pair out;
    if(z_img < -h) {
        out = argtheta3(-zuse, tau, passes);
    } else if(z_img >= h) {
        pair zmin = zuse - tau;
        out = -2.0 * pi * _i_ * zmin + argtheta3(zmin, tau, passes) -
            _i_ * pi * tau;
    } else {
        out = calctheta3(zuse, tau);
    }
    return out;
}

pair dologtheta3(pair z, pair tau, int pass_in) {
    int passes = pass_in + 1;
//  if(passes > 800) {
//      error("passes > 800 (dologtheta3)");
//  }
    pair tau2;
    real rl = tau.x;
    if(rl >= 0) {
        tau2 = modulo(rl + 1.0, 2) - 1.0 + _i_ * tau.y;
    } else {
        tau2 = modulo(rl - 1.0, 2) + 1.0 + _i_ * tau.y;
    }
    rl = tau2.x;
    pair out;
    if(rl > 0.6) {
        out = dologtheta3(z + 0.5, tau2 - 1.0, passes);
    } else if(rl < -0.6) {
        out = dologtheta3(z + 0.5, tau2 + 1.0, passes);
    } else if(length(tau2) < 0.98 && tau2.y < 0.98) {
        pair tauprime = (-1.0, 0.0) / tau2;
        out = _i_ * pi * tauprime * z * z +
            dologtheta3(-z * tauprime, tauprime, passes) -
            log(sqrt(-_i_ * tau2));
    } else {
        out = argtheta3(z, tau2, passes);
    }
    return out;
}

pair M(pair z, pair tau) {
    return _i_ * pi * (z + tau / 4.0);
}

pair ljtheta1(pair z, pair tau) {
    return M(z/pi - 0.5, tau) + dologtheta3(z/pi -0.5 + 0.5 * tau, tau, 0);
}

pair ljtheta2(pair z, pair tau) {
    return M(z/pi, tau) + dologtheta3(z/pi + 0.5 * tau, tau, 0);
}

pair ljtheta3(pair z, pair tau) {
    return dologtheta3(z/pi, tau, 0);
}

pair ljtheta4(pair z, pair tau) {
    return dologtheta3(z/pi + 0.5, tau, 0);
}

pair jtheta1(pair z, pair tau) {
    return exp(ljtheta1(z, tau));
}

pair jtheta2(pair z, pair tau) {
    return exp(ljtheta2(z, tau));
}

pair jtheta3(pair z, pair tau) {
    return exp(ljtheta3(z, tau));
}

pair jtheta4(pair z, pair tau) {
    return exp(ljtheta4(z, tau));
}

// checks:
write(ljtheta1((2, 2), (0.5, 0.5))); // 2.347537-4.509974i
write(ljtheta2((2, 2), (0.5, 0.5))); // 2.707104-2.15378i
write(ljtheta3((2, 2), (0.5, 0.5))); // 2.732256-2.15378i
write(ljtheta4((2, 2), (0.5, 0.5))); // 2.199908-6.080771i

[stla]

Finally I took an example to do a package. Is it possible to submit it somewhere?

[johncbowman]

It would probably be best to port your code to C++ and build it into asy, for greater efficiency. An even better option would be to get these functions included into the GNU Scientific Library, which we can then access from Asymptote.

You should make use of the builtin function expi(theta)=(cos(theta),sin(theta)).

We will soon set up an Asymptote-aware wiki where third-party packages like this can be submitted and maintained.

[stla]

I don't know how to build C++ code in Asymptote. Is there a tuto somewhere? I didn't find by googling.