# Approximating a special case of the Riemann Theta function

This challenge is to write fast code that can perform a computationally difficult infinite sum.

Input

An n by n matrix P with integer entries that are smaller than 100 in absolute value. When testing I am happy to provide input to your code in any sensible format your code wants. The default will be one line per row of the matrix, space separated and provided on standard input.

P will be positive definite which implies it will always be symmetric. Other than that you don't really need to know what positive definite means to answer the challenge. It does however mean that there will actually be an answer to the sum defined below.

You do however need to know what a matrix-vector product is.

Output

Your code should compute the infinite sum:

to within plus or minus 0.0001 of the correct answer. Here Z is the set of integers and so Z^n is all possible vectors with n integer elements and e is the famous mathematical constant that approximately equals 2.71828. Note that the value in the exponent is simply a number. See below for an explicit example.

How does this relate to the Riemann Theta function?

In the notation of this paper on approximating the Riemann Theta function we are trying to compute . Our problem is a special case for at least two reasons.

• We set the initial parameter called z in the linked paper to 0.
• We create the matrix P in such a way that the minimum size of an eigenvalue is 1. (See below for how the matrix is created.)

Examples

P = [[ 5.,  2.,  0.,  0.],
[ 2.,  5.,  2., -2.],
[ 0.,  2.,  5.,  0.],
[ 0., -2.,  0.,  5.]]

Output: 1.07551411208


In more detail, let us see just one term in the sum for this P. Take for example just one term in the sum:

and x^T P x = 30. Notice thate^(-30) is about 10^(-14) and so is unlikely to be important for getting the correct answer up to the given tolerance. Recall that the infinite sum will actually use every possible vector of length 4 where the elements are integers. I just picked one to give an explicit example.

P = [[ 5.,  2.,  2.,  2.],
[ 2.,  5.,  4.,  4.],
[ 2.,  4.,  5.,  4.],
[ 2.,  4.,  4.,  5.]]

Output = 1.91841190706

P = [[ 6., -3.,  3., -3.,  3.],
[-3.,  6., -5.,  5., -5.],
[ 3., -5.,  6., -5.,  5.],
[-3.,  5., -5.,  6., -5.],
[ 3., -5.,  5., -5.,  6.]]

Output = 2.87091065342

P = [[6., -1., -3., 1., 3., -1., -3., 1., 3.],
[-1., 6., -1., -5., 1., 5., -1., -5., 1.],
[-3., -1., 6., 1., -5., -1., 5., 1., -5.],
[1., -5., 1., 6., -1., -5., 1., 5., -1.],
[3., 1., -5., -1., 6., 1., -5., -1., 5.],
[-1., 5., -1., -5., 1., 6., -1., -5., 1.],
[-3., -1., 5., 1., -5., -1., 6., 1., -5.],
[1., -5., 1., 5., -1., -5., 1., 6., -1.],
[3., 1., -5., -1., 5., 1., -5., -1., 6.]]

Output: 8.1443647932

P = [[ 7.,  2.,  0.,  0.,  6.,  2.,  0.,  0.,  6.],
[ 2.,  7.,  0.,  0.,  2.,  6.,  0.,  0.,  2.],
[ 0.,  0.,  7., -2.,  0.,  0.,  6., -2.,  0.],
[ 0.,  0., -2.,  7.,  0.,  0., -2.,  6.,  0.],
[ 6.,  2.,  0.,  0.,  7.,  2.,  0.,  0.,  6.],
[ 2.,  6.,  0.,  0.,  2.,  7.,  0.,  0.,  2.],
[ 0.,  0.,  6., -2.,  0.,  0.,  7., -2.,  0.],
[ 0.,  0., -2.,  6.,  0.,  0., -2.,  7.,  0.],
[ 6.,  2.,  0.,  0.,  6.,  2.,  0.,  0.,  7.]]

Output = 3.80639191181


Score

I will test your code on randomly chosen matrices P of increasing size.

Your score is simply the largest n for which I get a correct answer in less than 30 seconds when averaged over 5 runs with randomly chosen matrices P of that size.

What about a tie?

If there is a tie, the winner will be the one whose code runs fastest averaged over 5 runs. In the event that those times are also equal, the winner is the first answer.

How will the random input be created?

1. Let M be a random m by n matrix with m<=n and entries which are -1 or 1. In Python/numpy M = np.random.choice([0,1], size = (m,n))*2-1 . In practice I will set m to be about n/2.
2. Let P be the identity matrix + M^T M. In Python/numpy P =np.identity(n)+np.dot(M.T,M). We are now guaranteed that P is positive definite and the entries are in a suitable range.

Note that this means that all eigenvalues of P are at least 1, making the problem potentially easier than the general problem of approximating the Riemann Theta function.

Languages and libraries

You can use any language or library you like. However, for the purposes of scoring I will run your code on my machine so please provide clear instructions for how to run it on Ubuntu.

My Machine The timings will be run on my machine. This is a standard Ubuntu install on an 8GB AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

• n = 47 in C++ by Ton Hospel
• n = 8 in Python by Maltysen
• It may be worth mentioning that a positive definite matrix is symmetric by definition. – 2012rcampion Apr 4 '16 at 13:03
• @2012rcampion Thanks. Added. – user9206 Apr 4 '16 at 13:06
• Ok, maybe this is a dumb question, but I've stared at this for ages and I cannot figure out how you got an x of [-1,0,2,1]. Can you elaborate on this? (Hint: I am not a math guru) – wnnmaw Apr 4 '16 at 19:07
• @wnnmaw Sorry for being confusing. The sum has one term for every possible vector x of length 4 in this case. [-1,0,2,1] is just one I picked at random to show explicitly what the term would be in that case. – user9206 Apr 4 '16 at 19:09
• @Lembik The way you generate the SPD matrices implies that no singular value has ever an absolute value below 1. Can we use that knowledge? – flawr Apr 5 '16 at 18:43

# C++

No more naive approach. Only evaluate inside the ellipsoid.

Uses the armadillo, ntl, gsl and pthread libraries. Install using

apt-get install libarmadillo-dev libntl-dev libgsl-dev


Compile the program using something like:

g++ -Wall -std=c++11 -O3 -fno-math-errno -funsafe-math-optimizations -ffast-math -fno-signed-zeros -fno-trapping-math -fomit-frame-pointer -march=native -s infinity.cpp -larmadillo -lntl -lgsl -lpthread -o infinity


On some systems you may need to add -lgslcblas after -lgsl.

Run with the size of the matrix followed by the elements on STDIN:

./infinity < matrix.txt


matrix.txt:

4
5  2  0  0
2  5  2 -2
0  2  5  0
0 -2  0  5


Or to try a precision of 1e-5:

./infinity -p 1e-5 < matrix.txt


infinity.cpp:

// Based on http://arxiv.org/abs/nlin/0206009

#include <iostream>
#include <vector>
#include <stdexcept>
#include <cstdlib>
#include <cmath>
#include <string>
#include <future>
#include <chrono>

using namespace std;

#include <getopt.h>

using namespace arma;

#include <NTL/mat_ZZ.h>
#include <NTL/LLL.h>

using namespace NTL;

#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

double const EPSILON = 1e-4;       // default precision
double const GROW    = 2;          // By how much we grow the ellipsoid volume
double const UPSCALE = 1e9;        // lattice reduction, upscale real to integer
double const THREAD_SEC = 0.1;     // Use threads if need more time than this
double const RADIUS_MAX = 1e6;     // Maximum radius used in root finding
double const RADIUS_INTERVAL = 1e-6; // precision of target radius
int const ITER_MAX = 1000;         // Maximum iterations in root finding
unsigned long POINTS_MIN = 1000;   // Minimum points before getting fancy

struct Result {
Result& operator+=(Result const& add) {
return *this;
}

friend Result operator-(Result const& left, Result const& right) {
return Result{left.sum - right.sum,
left.elapsed - right.elapsed,
left.points - right.points};
}

double sum, elapsed;
unsigned long points;
};

struct Params {
double half_rho, half_N, epsilon;
};

double fill_factor_error(double r, void *void_params) {
auto params = static_cast<Params*>(void_params);
r -= params->half_rho;
return gsl_sf_gamma_inc(params->half_N, r*r) - params->epsilon;
}

// Calculate radius needed for target precision
double radius(int N, double rho, double lat_det, double epsilon) {
Params params;

params.half_rho = rho / 2.;
params.half_N   = N   / 2.;
params.epsilon = epsilon*lat_det*gsl_sf_gamma(params.half_N)/pow(M_PI, params.half_N);

// Calculate minimum allowed radius
auto r = sqrt(params.half_N)+params.half_rho;
auto val = fill_factor_error(r, &params);
cout << "Minimum R=" << r << " -> " << val << endl;

if (val > 0) {
// The minimum radius is not good enough. Work out a better one by
// finding the root of a tricky function
auto low  = r;
auto high = RADIUS_MAX * 2 * params.half_rho;
auto val = fill_factor_error(high, &params);
if (val >= 0)
throw(logic_error("huge RADIUS_MAX is still not big enough"));

gsl_function F;
F.function = fill_factor_error;
F.params   = &params;

auto T = gsl_root_fsolver_brent;
auto s = gsl_root_fsolver_alloc (T);
gsl_root_fsolver_set (s, &F, low, high);

int status = GSL_CONTINUE;
for (auto iter=1; status == GSL_CONTINUE && iter <= ITER_MAX; ++iter) {
gsl_root_fsolver_iterate (s);
low  = gsl_root_fsolver_x_lower (s);
high = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval(low, high, 0, RADIUS_INTERVAL  * 2 * params.half_rho);
}
r = gsl_root_fsolver_root(s);
gsl_root_fsolver_free(s);
if (status == GSL_CONTINUE)
throw(logic_error("Search for R did not converge"));
}
return r;
}

// Recursively walk down the ellipsoids in each dimension
void ellipsoid(int d, mat const& A, double const* InvD, mat& Accu,
Result& result, double r2) {
auto r = sqrt(r2);
auto offset = Accu(d, d);
// InvD[d] = 1/ A(d, d)
auto from = ceil((-r-offset) * InvD[d]);
auto to   = floor((r-offset) * InvD[d]);
for (auto v = from; v <= to; ++v) {
auto value  = v * A(d, d)+offset;
auto residu = r2 - value*value;
if (d == 0) {
result.sum += exp(residu);
++result.points;
} else {
for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + v * A(d, i);
ellipsoid(d-1, A, InvD, Accu, result, residu);
}
}
}

// Specialised version of ellipsoid() that will only process points an octant
void ellipsoid(int d, mat const& A, double const* InvD, mat& Accu,
Result& result, double r2, unsigned int octant) {
auto r = sqrt(r2);
auto offset = Accu(d, d);
// InvD[d] = 1/ A(d, d)
long from = ceil((-r-offset) * InvD[d]);
long to   = floor((r-offset) * InvD[d]);
auto points = to-from+1;
auto base = from + points/2;
if (points & 1) {
auto value = base * A(d, d) + offset;
auto residu = r2 - value * value;
if (d == 0) {
if ((octant & (octant - 1)) == 0) {
result.sum += exp(residu);
++result.points;
}
} else {
for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + base * A(d, i);
ellipsoid(d-1, A, InvD, Accu, result, residu, octant);
}
++base;
}
if ((octant & 1) == 0) {
to = from + points / 2 - 1;
base = from;
}
octant /= 2;
for (auto v = base; v <= to; ++v) {
auto value = v * A(d,d)+offset;
auto residu = r2 - value*value;
if (d == 0) {
if ((octant & (octant - 1)) == 0) {
result.sum += exp(residu);
++result.points;
}
} else {
for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + v * A(d, i);
if (octant == 1)
ellipsoid(d-1, A, InvD, Accu, result, residu);
else
ellipsoid(d-1, A, InvD, Accu, result, residu, octant);
}
}
}

// Prepare call to ellipsoid()
Result sym_ellipsoid(int N, mat const& A, const vector<double>& InvD, double r,
unsigned int octant = 1) {
auto start = chrono::steady_clock::now();
auto r2 = r*r;

mat Accu(N, N);
Accu.row(N-1).zeros();

Result result{0, 0, 0};
// 2*octant+1 forces the points into the upper half plane, skipping 0
// This way we use the lattice symmetry and calculate only half the points
ellipsoid(N-1, A, &InvD[0], Accu, result, r2, 2*octant+1);
// Compensate for the extra factor exp(r*r) we always add in ellipsoid()
result.sum /= exp(r2);
auto end = chrono::steady_clock::now();
result.elapsed = chrono::duration<double>{end-start}.count();

return result;
}

// Prepare multithreaded use of sym_ellipsoid(). Each thread gets 1 octant
Result sym_ellipsoid_t(int N, mat const& A, const vector<double>& InvD, double r, unsigned int nr_threads) {

vector<future<Result>> results;
results.emplace_back(async(launch::async, sym_ellipsoid, N, ref(A), ref(InvD), r, i));
auto result = sym_ellipsoid(N, A, InvD, r, nr_threads);
for (auto i=0U; i<nr_threads-1; ++i) result += results[i].get();
return result;
}

int main(int argc, char* const* argv) {
cout.precision(12);

double epsilon    = EPSILON; // Target absolute error
bool inv_modular  = true;    // Use modular transform to get the best matrix
bool lat_reduce   = true;    // Use lattice reduction to align the ellipsoid
bool conservative = false;   // Use provable error bound instead of a guess
bool eigen_values = false;   // Show eigenvalues

int option_char;
while ((option_char = getopt(argc, argv, "p:n:MRce")) != EOF)
switch (option_char) {
case 'p': epsilon      = atof(optarg); break;
case 'n': threads_max  = atoi(optarg); break;
case 'M': inv_modular  = false;        break;
case 'R': lat_reduce   = false;        break;
case 'c': conservative = true;         break;
case 'e': eigen_values = true;         break;
default:
cerr << "usage: " << argv[0] << " [-p epsilon] [-n threads] [-M] [-R] [-e] [-c]" << endl;
exit(EXIT_FAILURE);
}
if (optind < argc) {
cerr << "Unexpected argument" << endl;
exit(EXIT_FAILURE);
}
if (threads_max < 1) threads_max = 1;
cout << "Using up to " << threads_max << " threads" << endl;

int N;
cin >> N;

mat P(N, N);
for (auto& v: P) cin >> v;

if (eigen_values) {
vec eigval = eig_sym(P);
cout << "Eigenvalues:\n" << eigval << endl;
}

// Decompose P = A * A.t()
mat A = chol(P, "lower");

// Calculate lattice determinant
double lat_det = 1;
for (auto i=0; i<N; ++i) {
if (A(i,i) <= 0) throw(logic_error("Diagonal not Positive"));
lat_det *= A(i,i);
}
cout << "Lattice determinant=" << lat_det << endl;

auto factor = lat_det / pow(M_PI, N/2.0);
if (inv_modular && factor < 1) {
epsilon *= factor;
cout << "Lattice determinant is small. Using inverse instead. Factor=" << factor << endl;
P = M_PI * M_PI * inv(P);
A = chol(P, "lower");
// We could simple calculate the new lat_det as pow(M_PI,N)/lat_det
lat_det = 1;
for (auto i=0; i<N; ++i) {
if (A(i,i) <= 0) throw(logic_error("Diagonal not Positive"));
lat_det *= A(i,i);
}
cout << "New lattice determinant=" << lat_det << endl;
} else
factor = 1;

// Prepare for lattice reduction.
// Since the library works on integer lattices we will scale up our matrix
double min = INFINITY;
for (auto i=0; i<N; ++i) {
for (auto j=0; j<N;++j)
if (A(i,j) != 0 && abs(A(i,j) < min)) min = abs(A(i,j));
}

auto upscale = UPSCALE/min;
mat_ZZ a;
a.SetDims(N,N);
for (auto i=0; i<N; ++i)
for (auto j=0; j<N;++j) a[i][j] = to_ZZ(A(i,j)*upscale);

// Finally do the actual lattice reduction
mat_ZZ u;
auto rank = G_BKZ_FP(a, u);
if (rank != N) throw(logic_error("Matrix is singular"));
mat U(N,N);
for (auto i=0; i<N;++i)
for (auto j=0; j<N;++j) U(i,j) = to_double(u[i][j]);

// There should now be a short lattice vector at row 0
ZZ sum = to_ZZ(0);
for (auto j=0; j<N;++j) sum += a[0][j]*a[0][j];
auto rho = sqrt(to_double(sum))/upscale;
cout << "Rho=" << rho << " (integer square " <<
rho*rho << " ~ " <<
static_cast<int>(rho*rho+0.5) << ")" << endl;

// Lattice reduction doesn't gain us anything conceptually.
// The same number of points is evaluated for the same exponential values
// However working through the ellipsoid dimensions from large lattice
// base vectors to small makes ellipsoid() a *lot* faster
if (lat_reduce) {
mat B = U * A;
P = B * B.t();
A = chol(P, "lower");
if (eigen_values) {
vec eigval = eig_sym(P);
cout << "New eigenvalues:\n" << eigval << endl;
}
}

vector<double> InvD(N);;
for (auto i=0; i<N; ++i) InvD[i] = 1 / A(i, i);

// Calculate radius needed for target precision
auto r = radius(N, rho, lat_det, epsilon);
cout << "Safe R=" << r << endl;

Result result;
if (conservative) {
// Walk all points inside the ellipsoid with transformed radius r
result = sym_ellipsoid_t(N, A, InvD, r, nr_threads);
} else {
// First grow the radius until we saw POINTS_MIN points or reach the
double i = floor(N * log2(r/rho) / log2(GROW));
if (i < 0) i = 0;
auto R = r * pow(GROW, -i/N);
cout << "Initial R=" << R << endl;
result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
auto max_new_points = result.points;
while (--i >= 0 && result.points < POINTS_MIN) {
R = r * pow(GROW, -i/N);
auto change = result;
result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
change = result - change;

if (change.points > max_new_points) max_new_points = change.points;
}

// Now we have enough points that it's worth bothering to use threads
while (--i >= 0) {
R = r * pow(GROW, -i/N);
auto change = result;
result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
change = result - change;
// This is probably too crude and might misestimate the error
// I've never seen it fail though
if (change.points > max_new_points) {
max_new_points = change.points;
if (change.sum < epsilon/2) break;
}
}
cout << "Final R=" << R << endl;
}

// We calculated half the points and skipped 0.
result.sum = 2*result.sum+1;

// Modular transform factor
result.sum /= factor;

// Report result
cout <<
"Evaluated " << result.points << " points\n" <<
"Sum = " << result.sum << endl;
}

• This is very impressive and much better than the naive approach in my view. I look forward to the documentation :) – user9206 Apr 6 '16 at 13:57
• @TonHospel Can you tell us a bit more about how you come up with the bounds? – flawr Apr 7 '16 at 18:08
• I am using Arch Linux and needed the -lgslcblas flag to compile. Amazing answer by the way! – Rhyzomatic Apr 12 '16 at 5:26

# Python 3

12 secs n=8 on my computer, ubuntu 4 core.

Really naive, have no clue what I'm doing.

from itertools import product
from math import e

P = [[ 6., -3.,  3., -3.,  3.],
[-3.,  6., -5.,  5., -5.],
[ 3., -5.,  6., -5.,  5.],
[-3.,  5., -5.,  6., -5.],
[ 3., -5.,  5., -5.,  6.]]

N = 2

n = [1]

while e** -n[-1] > 0.0001:
n = []
for x in product(list(range(-N, N+1)), repeat = len(P)):
n.append(sum(k[0] * k[1] for k in zip([sum(j[0] * j[1] for j in zip(i, x)) for i in P], x)))
N += 1

print(sum(e** -i for i in n))


This will keep on increasing the range of Z that it uses until it gets a good enough answer. I wrote my own matrix multiplication, prolly should use numpy.

• Thanks ! Can you show some outputs and timings on your computer? – user9206 Apr 5 '16 at 6:41
• Your code runs in pypy which is great and is fast. Unfortunately, [[6.0, -1.0, -3.0, 1.0, 3.0, -1.0, -3.0, 1.0, 3.0], [-1.0, 6.0, -1.0, -5.0, 1.0, 5.0, -1.0, -5.0, 1.0], [-3.0, -1.0, 6.0, 1.0, -5.0, -1.0, 5.0, 1.0, -5.0], [1.0, -5.0, 1.0, 6.0, -1.0, -5.0, 1.0, 5.0, -1.0], [3.0, 1.0, -5.0, -1.0, 6.0, 1.0, -5.0, -1.0, 5.0], [-1.0, 5.0, -1.0, -5.0, 1.0, 6.0, -1.0, -5.0, 1.0], [-3.0, -1.0, 5.0, 1.0, -5.0, -1.0, 6.0, 1.0, -5.0], [1.0, -5.0, 1.0, 5.0, -1.0, -5.0, 1.0, 6.0, -1.0], [3.0, 1.0, -5.0, -1.0, 5.0, 1.0, -5.0, -1.0, 6.0]] gives just the wrong answer. – user9206 Apr 5 '16 at 7:22
• 8.1443647932-8.14381938863 = 0.00054540457 > 0.0001. – user9206 Apr 5 '16 at 7:22
• @Maltysen Your program does only check whether the last term is smaller than the given accuracy. But the error you make is by far bigger, as you'd have to consider the sum of all the other terms for the error too! – flawr Apr 5 '16 at 18:01