Riemann Theta函数特例的逼近

Question

这个挑战是编写能够执行计算困难的无限和的快速代码。

输入

n by n矩阵P，其整数项在绝对值上小于100。在测试时，我很高兴以您的代码想要的任何合理格式向您的代码提供输入。默认情况是矩阵的每一行一行，在标准输入中分隔和提供空格。

P将是正定，这意味着它总是对称的。除此之外，你不需要真正知道什么是肯定的含义，以应对挑战。然而，这确实意味着，实际上将有一个答案的总和定义如下。

但是，您确实需要知道矩阵向量积是什么。

输出

您的代码应该计算无限和：

在正确答案的正负0.0001以内。这里，Z是整数集合，所以Z^n是所有可能的向量，n整数元素，e是大约等于2.71828的著名数学常数。注意，指数中的值只是一个数字。有关显式示例，请参阅下面。

这与Riemann Theta函数有什么关系？

在本文关于Riemann Theta函数的逼近的表示法中，我们试图计算

。我们的问题是一个特例，至少有两个原因。

• 我们将链接文件中名为z的初始参数设置为0。
• 我们以这样的方式创建矩阵P，使得特征值的最小大小是1。(有关如何创建矩阵，请参见下面。)

示例

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

Output: 1.07551411208

更详细地说，让我们只看到这个P的和中的一个项。例如，在和中的一个项：

和x^T P x = 30。请注意，e^(-30)是关于10^(-14)的，因此不太可能对获得正确的答案达到给定的容忍度很重要。回想一下，无限和实际上将使用长度为4的每个可能向量，其中元素是整数。我只是选了一个例子来明确说明。

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

评分

我会在随机选择的矩阵P上测试你的代码。

你的分数是最大的n，我在不到30秒的时间内得到了一个正确的答案，平均超过5次，随机选择的矩阵为P。

领带怎么样？

如果有一个平局，胜利者将是谁的代码运行最快的平均超过5次。如果这些时间也是平等的，胜利者就是第一个答案。

如何创建随机输入？

1. 设M是带有m<=n的n个矩阵的随机m，在Python/numpy M = np.random.choice([0,1], size = (m,n))*2-1中，条目是-1或1。在实践中，我将m设置为关于n/2。
2. 设P是Python/numpy P =np.identity(n)+np.dot(M.T,M)中的恒等矩阵+ M^T。现在我们保证P是正定的，并且条目在适当的范围内。

请注意，这意味着P的所有特征值至少为1，这使得问题可能比一般的Riemann Theta函数近似问题更容易解决。

语言和库

您可以使用任何您喜欢的语言或库。但是，为了打分，我将在我的机器上运行您的代码，所以请提供关于如何在Ubuntu上运行它的明确说明。

我的机器计时将在我的机器上运行。这是一个标准的Ubuntu安装在一个8GB的AMD FX-8350八核处理器上.这也意味着我需要能够运行您的代码。

领先答案

• n = 47 in C++
• 马尔蒂森在Python中的n = 8

Answer 1

C++

别再天真了。只在椭球内求值。

使用armadillo、ntl、gsl和p线程库。安装使用

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

使用以下内容编译程序：

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

在某些系统上，您可能需要在-lgslcblas之后添加-lgsl。

以矩阵的大小运行，后面跟着STDIN上的元素：

./infinity < matrix.txt

matrix.txt：

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

或者尝试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 <thread>
#include <future>
#include <chrono>

using namespace std;

#include <getopt.h>

#include <armadillo>

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) {
        sum     += add.sum;
        elapsed += add.elapsed;
        points  += add.points;
        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, ¶ms);
    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, ¶ms);
        if (val >= 0)
            throw(logic_error("huge RADIUS_MAX is still not big enough"));

        gsl_function F;
        F.function = fill_factor_error;
        F.params   = ¶ms;

        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) {
    nr_threads = pow(2, ceil(log2(nr_threads)));

    vector<future<Result>> results;
    for (auto i=nr_threads+1; i<2*nr_threads; ++i)
        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) {
    cin.exceptions(ios::failbit | ios::badbit);
    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  threads_max  = thread::hardware_concurrency();

    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;
    threads_max = pow(2, ceil(log2(threads_max)));
    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;

    auto nr_threads = threads_max;
    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
        // target radius
        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);
        nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
        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);
            nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
            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);
            nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
            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;
}

Answer 2

Python 3

12秒n=8在我的电脑上，ubuntu4core.

太天真了，不知道我在做什么。

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))

这将继续增加它使用的Z的范围，直到得到一个足够好的答案为止。我写了我自己的矩阵乘法，很明显应该用numpy。

Answer 3

Python，numpy

使用numpy的示例程序如下所示。

太天真了，不知道我在做什么。

import numpy as np
from itertools import product
from math import e
from time import time

n=9
is_randomized=False

if is_randomized:
    m = n // 2
    M = np.random.choice([0, 1], size=(m, n)) * 2 - 1
    P = np.identity(n) + np.dot(M.T, M)
    print(P)
else:
    if n==4:
        P = [[5., 2., 2., 2.],
             [2., 5., 4., 4.],
             [2., 4., 5., 4.],
             [2., 4., 4., 5.]]
        P = np.array(P)
        print(P)
    elif n==5:
        P = np.array([[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.]])
        P = np.array(P)
        print(P)
    elif n==9:
        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.]]
        P = np.array(P)
        print(P)


# estimate how much N should be chosen?  based on matrix P or size n ?
if n==9:
    N = 2
elif n==5:
    N = 4
elif n==4:
    N = 5
else:
    N = 3

n_list = []
start_time=time()
for x in product(range(-N, N+1), repeat=len(P)):
    x = np.array(x)
    summand = np.dot(np.dot(x, P), x.T)
    n_list.append(summand)
result = sum(e**-i for i in n_list)
print(result)
print(f"duration: {time()-start_time}")

8秒，n=9在我的电脑上。