问:
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

```<code>P = [[ 5.,  2.,  0.,  0.],
     [ 2.,  5.,  2., -2.],
     [ 0.,  2.,  5.,  0.],
     [ 0., -2.,  0.,  5.]]

Output: 1.07551411208
</code>```

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.

```<code>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
</code>```

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?

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.

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.

Leading answers

n = 47 in C++ by Ton Hospel

n = 8 in Python by Maltysen

答:
C++

No more naive approach. Only evaluate inside the ellipsoid.

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

```<code>apt-get install libarmadillo-dev libntl-dev libgsl-dev
</code>```

Compile the program using something like:

```<code>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
</code>```

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

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

```<code>./infinity < matrix.txt
</code>```

matrix.txt:

```<code>4
5  2  0  0
2  5  2 -2
0  2  5  0
0 -2  0  5
</code>```

Or to try a precision of 1e-5:

```<code>./infinity -p 1e-5 < matrix.txt
</code>```

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, &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) {
    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;
}
```

Body:
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 <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, &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) {
    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;
}


问:
Are the brackets fully matched?
You must write a program or function that takes a string of brackets and outputs whether or not that string is fully matched. Your program should print a truthy or falsy value, and IO can be in any reasonable format.

Rules and definitions:

For the purpose of this challenge, a "bracket" is any of these characters: ()[]{}<>.

A pair of brackets is considered "matched" if the opening and closing brackets are in the right order and have no characters inside of them, such as

```<code>()
[]{}
</code>```

Or if every subelement inside of it is also matched.

```<code>[()()()()]
{<[]>}
(()())
</code>```

Subelements can also be nested several layers deep.

```<code>[(){<><>[()]}<>()]
<[{((()))}]>
</code>```

A string is considered "Fully matched" if and only if:

Every single character is a bracket,

Each pair of brackets has the correct opening and closing bracket and in the right order, and

Each bracket is matched.

You may assume the input will only contain printable ASCII.

Test IO

Here are some inputs that should return a truthy value:

```<code>()
[](){}<>
(((())))
({[<>]})
[{()<>()}[]]
[([]{})<{[()<()>]}()>{}]
</code>```

And here are some outputs that should return a falsy value:

```<code>(               Has no closing ')'
}{              Wrong order
(<)>            Each pair contains only half of a matched element
(()()foobar)    Contains invalid characters
[({}<>)>        The last bracket should be ']' instead of '>'
(((()))         Has 4 opening brackets, but only 3 closing brackets.
</code>```

As usual, this is code-golf, so standard loopholes apply, and shortest answer in bytes wins.

答:
Regex (PCRE flavor), <s>41</s> 37 bytes

```<code>^((<(?1)>|{(?1)}|\[(?1)]|\((?1)\))*)$
</code>```

Just a standard solution with recursive regex.

Thanks jimmy23013 for shaving off 4 bytes

ContentLicense: CC BY-SA 3.0

Title: Are the brackets fully matched?

Body:
Regex (PCRE flavor), <s>41</s> 37 bytes

^((<(?1)>|{(?1)}|\[(?1)]|\((?1)\))*)$

Just a standard solution with recursive regex.

Thanks jimmy23013 for shaving off 4 bytes


问:
Double-duty quine: quine!dlroW ,olleH
Write a program that will output its own source code when run, and nothing else. Seems easy, right?

The catch is that when the source code is reversed, your program must output "Hello, World!" exactly, without the quotes.

This is code-golf, so lowest byte count wins.

Edit: Your quine must be a proper quine.

答:
GolfScript, 29 28 bytes

```"`'.~]'#\"!dlroW ,olleH".~]

```

It has one trailing newline. Try it here.

Reversed:

```

]~."Hello, World!"\#']~.'`"
```

Try it here.

Title: Double-duty quine: quine!dlroW ,olleH

Body:
GolfScript, 29 28 bytes

"`'.~]'#\"!dlroW ,olleH".~]

It has one trailing newline. Try it here.

Reversed:

]~."Hello, World!"\#']~.'`"

Try it here.


问:
Magic number of a given length
Your program must take an input (n for the purpose of description) and output all permutations of a number that is n digits long with no repeating digits, where each of the digits preceding and including its index are divisible by the place in the number that it falls.

You can read about magic numbers here.

Rules:

1 <= n <= 10

No digits may be repeated

The leading 0 must be present (if applicable)

The 1st through xth digit of the number (starting with the first character as 1) must be divisible by x, i.e. in 30685, 3 is divisible by 1, 30 is divisible by 2, 306 is divisible by 3, 3068 is divisible by 4, and 30685 is divislbe by 5.

The program must take an integer as input (through the command line, as a function argument, etc.) and print all permutations that satisfy the rules.

Output must be separated by 1 or more white space character

Permutations may start and with zero (so they're not technically magic numbers).

The order of output does not matter

You do not need to handle unexpected input

Least characters in bytes wins

Examples

Given 1:

```<code>0
1
2
3
4
5
6
7
8
9
</code>```

Given 2:

```<code>02
04
06
08
10
12
14
16
18
20
24
26
28
30
32
34
36
38
40
42
46
48
50
52
54
56
58
60
62
64
68
70
72
74
76
78
80
82
84
86
90
92
94
96
98
</code>```

Given 10:

```<code>3816547290
</code>```

Credit to Pizza Hut & John H. Conway for the original puzzle (Option A).
Thanks to @Mego and @sp3000 for their links.

答:
Python, 132 bytes

lambda n:[x for x in map(("{:0%s}"%n).format,(range(10**n)))if all(int(x[:i])%i<1and len(set(x))==len(x)for i in range(1,len(x)+1))]
```

Dropped 26 bytes by getting rid of itertools, thanks to Sp3000 for making me realize I shouldn't be using it.

Dropped 2 bytes by using a list comprehension rather than filter, thanks again to Sp3000 for the tip.

Try it online: Python 2, Python 3

Title: Magic number of a given length

Body:
Python, 132 bytes

lambda n:[x for x in map(("{:0%s}"%n).format,(range(10**n)))if all(int(x[:i])%i<1and len(set(x))==len(x)for i in range(1,len(x)+1))]

Dropped 26 bytes by getting rid of itertools, thanks to Sp3000 for making me realize I shouldn't be using it.

Dropped 2 bytes by using a list comprehension rather than filter, thanks again to Sp3000 for the tip.

Try it online: Python 2, Python 3


问:
Rectangle Detection
Write a program or function that takes in a multiline string of 0's and 1's. No other characters will be in the string and the string will always be rectangular (all lines will have same number of characters), with dimensions as small as 1×1, but otherwise the 0's and 1's may be arranged arbitrarily.

You may assume the string has an optional trailing newline, and if desired you may use any two distinct printable ASCII characters in place of 0 and 1.

Print or return a truthy value if all of the path connected regions of both 0's and 1's in the string are solid rectangles, else output a falsy value.

A path connected region of 0's means that from any one 0 in the region, all the other 0's can be reached by only moving up, down, left and right to other 0's (and not moving diagonally, not moving to any 1, and not moving outside the string bounds). The same idea applies to 1 path connected regions.

A solid rectangle of 0's means the entire area of the rectangle is filled with 0's and no 1's. The same idea applies to 1 solid rectangles.

The shortest code in bytes wins. Tiebreaker is earlier answer.

(Note that the string does not wrap around with toroidal boundary conditions.)

Examples

1) This input string has 3 path connected regions (2 for 0 and 1 for 1). Only the bottom right 00 region is a solid rectangle though, so the output would be falsy.

```<code>0011
0111
0100
</code>```

2) This input string has 4 path connected regions (2 for both 0 and 1). All of them are solid rectangles so the output would be truthy.

```<code>0011
0011
1100
</code>```

3) This input has 2 path connected regions, but only one of them is a solid rectangle, so the output would be falsy.

```<code>00000000
01111110
00000000
</code>```

4) This input has only 1 path connected region and is trivially a solid rectangle, so the output is truthy.

```<code>11111111
11111111
11111111
</code>```

Test Cases

A T just below the input string means truthy, F means falsy.

```<code>0
T

1
T

00
T

01
T

10
T

11
T

0000000
T

1111111
T

011100100100101100110100100100101010100011100101
T

00
11
T

01
10
T

01
11
F

00
01
F

11
11
T

110
100
F

111
000
T

111
101
111
F

101
010
101
T

1101
0010
1101
0010
T

1101
0010
1111
0010
F

0011
0111
0100
F

0011
0011
1100
T

00000000
01111110
00000000
F

11111111
11111111
11111111
T

0000001111
0000001111
T

0000001111
0000011111
F

0000001111
1000001111
F

1000001111
1000001111
T

1110100110101010110100010111011101000101111
1010100100101010100100010101010101100101000
1110100110010010110101010111010101010101011
1010100100101010010101010110010101001101001
1010110110101010110111110101011101000101111
F</code>```

答:
JavaScript (ES6), 64 bytes

s=>(a=s.split`
`).every(l=>l==a[0]|l==a[0].replace(/./g,n=>n^1))
```

Based on @LevelRiverSt's observation that each line must either be the same or the opposite of the first.

Title: Rectangle Detection

Body:
JavaScript (ES6), 64 bytes

s=>(a=s.split`
`).every(l=>l==a[0]|l==a[0].replace(/./g,n=>n^1))

Based on @LevelRiverSt's observation that each line must either be the same or the opposite of the first.


问:
Palindromizing the strings

Introduction

For the ones who don't know, a palindrome is when a string is equal to the string backwards (with exception to interpunction, spaces, etc.). An example of a palindrome is:

```<code>abcdcba
</code>```

If you reverse this, you will end up with:

```<code>abcdcba
</code>```

Which is the same. Therefore, we call this a palindrome. To palindromize things, let's take a look at an example of a string:

```<code>adbcb
</code>```

This is not a palindrome. To palindromize this, we need to merge the reversed string into the initial string at the right of the initial string, leaving both versions intact. The shorter, the better.

The first thing we can try is the following:

```<code>adbcb
bcbda
^^ ^^
</code>```

Not all characters match, so this is not the right position for the reversed string. We go one step to the right:

```<code>adbcb
 bcbda
 ^^^^
</code>```

This also doesn't match all characters. We go another step to the right:

```<code>adbcb
  bcbda
</code>```

This time, all characters match. We can merge both string leaving the intact. The final result is:

```<code>adbcbda
</code>```

This is the palindromized string.

The Task

Given a string (with at least one character) containing only lowercase letters (or uppercase, if that fits better), output the palindromized string.

Test cases

```<code>Input     Output

abcb      abcba
hello     hellolleh
bonobo    bonobonob
radar     radar
hex       hexeh
</code>```

This is code-golf, so the submission with the least amount of bytes wins!

答:
Oracle SQL 11.2, 195 bytes

```<code>SELECT MIN(p)KEEP(DENSE_RANK FIRST ORDER BY LENGTH(p))FROM(SELECT:1||SUBSTR(REVERSE(:1),LEVEL+1)p FROM DUAL WHERE SUBSTR(:1,-LEVEL,LEVEL)=SUBSTR(REVERSE(:1),1,LEVEL)CONNECT BY LEVEL<=LENGTH(:1));
</code>```

Un-golfed

```<code>SELECT MIN(p)KEEP(DENSE_RANK FIRST ORDER BY LENGTH(p))
FROM (
       SELECT :1||SUBSTR(REVERSE(:1),LEVEL+1)p
       FROM   DUAL
       WHERE  SUBSTR(:1,-LEVEL,LEVEL)=SUBSTR(REVERSE(:1),1,LEVEL)
       CONNECT BY LEVEL<=LENGTH(:1)
     );
</code>```

Title: Palindromizing the strings

Body:
Oracle SQL 11.2, 195 bytes

SELECT MIN(p)KEEP(DENSE_RANK FIRST ORDER BY LENGTH(p))FROM(SELECT:1||SUBSTR(REVERSE(:1),LEVEL+1)p FROM DUAL WHERE SUBSTR(:1,-LEVEL,LEVEL)=SUBSTR(REVERSE(:1),1,LEVEL)CONNECT BY LEVEL<=LENGTH(:1));

Un-golfed

SELECT MIN(p)KEEP(DENSE_RANK FIRST ORDER BY LENGTH(p))
FROM (
       SELECT :1||SUBSTR(REVERSE(:1),LEVEL+1)p
       FROM   DUAL
       WHERE  SUBSTR(:1,-LEVEL,LEVEL)=SUBSTR(REVERSE(:1),1,LEVEL)
       CONNECT BY LEVEL<=LENGTH(:1)
     );


问:
Tips for golfing in Racket / Scheme
What general tips do you have for golfing in Racket / Scheme? I'm looking for ideas that can be applied to code golf problems in general that are at least somewhat specific to Racket / Scheme (e.g. "remove comments" is not an answer).

I'm aware Scheme and Racket (formerly PLT Scheme) are technically different languages but are quite similar in many ways and much code (I suspect) will run mostly as intended in either. If your tip only applies to one of the aforementioned languages, note as such.

答:
When using Racket, bind variables using λ to shave off a few bytes. In Scheme, lambda makes this trick not applicable, unless one is binding four or more variables.

Example: One variable saves 2 bytes over let/define

```<code>(define n 55)(* n n) ; 20 bytes

(let([n 55])(* n n)) ; 20 bytes

((λ(n)(* n n))55) ; 18 bytes
</code>```

Title: Tips for golfing in Racket / Scheme

Body:
When using Racket, bind variables using λ to shave off a few bytes. In Scheme, lambda makes this trick not applicable, unless one is binding four or more variables.

Example: One variable saves 2 bytes over let/define

(define n 55)(* n n) ; 20 bytes

(let([n 55])(* n n)) ; 20 bytes

((λ(n)(* n n))55) ; 18 bytes


问:
Diamondize a Matrix
Given a matrix, output a representation of the matrix where the top left element is on top, the anti-diagonal is the central row and the bottom right element is at the bottom.

For example, consider the following matrix:

```<code>1 2 3
4 5 6
7 8 9
</code>```

The diamond version of this matrix is:

```<code>  1
 4 2
7 5 3
 8 6
  9
</code>```

Inputs and outputs

An input matrix will be given as a list of lists (or anything similar in your language of choice). The output shall be a list of lists as well.

The matrices will only contain positive integers.

The input matrix will not necessarily be square.

The input matrix will be at least 1×1.

Test Cases

```<code>Input:  [[1]]
Output: [[1]]

Input:  [[1,2],[3,4]]
Output: [[1],[3,2],[4]]

Input:  [[1,2,3],[4,5,6]]
Output: [[1],[4,2],[5,3],[6]]

Input:  [[11,2,5],[3,99,3],[4,8,15],[16,23,42]]
Output: [[11],[3,2],[4,99,5],[16,8,3],[23,15],[42]]
</code>```

Scoring

This is code-golf, so the shortest answer in bytes wins.

答:
J, 6 bytes

```<code><@|./.
</code>```

This is an unnamed a monadic verb which takes a matrix and returns a list of antidiagonals:

```<code>   input =. i.3 4
   input
0 1  2  3
4 5  6  7
8 9 10 11

   <@|./. input
┌─┬───┬─────┬─────┬────┬──┐
│0│4 1│8 5 2│9 6 3│10 7│11│
└─┴───┴─────┴─────┴────┴──┘
</code>```

Test it here.

Explanation

/. is J's built-in to apply a function to each anti-diagonal. Unfortunately, these anti-diagonals are given in the opposite order of what we want here.

In <@|., we first apply |. which reverses the anti-diagonal and then < to box it (which is the only way to return a ragged array in J, since normal arrays are always rectangular, so the antidiagonals would be padded with zeroes).

Title: Diamondize a Matrix

Body:
J, 6 bytes

<@|./.

This is an unnamed a monadic verb which takes a matrix and returns a list of antidiagonals:

   input =. i.3 4
   input
0 1  2  3
4 5  6  7
8 9 10 11

   <@|./. input
┌─┬───┬─────┬─────┬────┬──┐
│0│4 1│8 5 2│9 6 3│10 7│11│
└─┴───┴─────┴─────┴────┴──┘

Test it here.

Explanation

/. is J's built-in to apply a function to each anti-diagonal. Unfortunately, these anti-diagonals are given in the opposite order of what we want here.

In <@|., we first apply |. which reverses the anti-diagonal and then < to box it (which is the only way to return a ragged array in J, since normal arrays are always rectangular, so the antidiagonals would be padded with zeroes).

ContentLicense: CC BY-SA 4.0


答:
Jelly, 7 bytes

```<code>ṚŒDṙZL$
</code>```

Try it online!

Explanation

```<code>Ṛ         Reverse the matrix vertically.
 ŒD       Get its diagonals. However these start from
          the main diagonal, not the corners.
    ZL$   Get the width of the input matrix.
   ṙ      Rotate the list of diagonals left by that many
          places to obtain the correct order.
</code>```

Title: Diamondize a Matrix

Body:
Given a matrix, output a representation of the matrix where the top left element is on top, the anti-diagonal is the central row and the bottom right element is at the bottom.

For example, consider the following matrix:

```<code>1 2 3
4 5 6
7 8 9
</code>```

The diamond version of this matrix is:

```<code>  1
 4 2
7 5 3
 8 6
  9
</code>```

Inputs and outputs

An input matrix will be given as a list of lists (or anything similar in your language of choice). The output shall be a list of lists as well.

The matrices will only contain positive integers.

The input matrix will not necessarily be square.

The input matrix will be at least 1×1.

Test Cases

```<code>Input:  [[1]]
Output: [[1]]

Input:  [[1,2],[3,4]]
Output: [[1],[3,2],[4]]

Input:  [[1,2,3],[4,5,6]]
Output: [[1],[4,2],[5,3],[6]]

Input:  [[11,2,5],[3,99,3],[4,8,15],[16,23,42]]
Output: [[11],[3,2],[4,99,5],[16,8,3],[23,15],[42]]
</code>```

Scoring

This is code-golf, so the shortest answer in bytes wins.

Body:
Jelly, 7 bytes

ṚŒDṙZL$

Try it online!

Explanation

Ṛ         Reverse the matrix vertically.
 ŒD       Get its diagonals. However these start from
          the main diagonal, not the corners.
    ZL$   Get the width of the input matrix.
   ṙ      Rotate the list of diagonals left by that many
          places to obtain the correct order.


问:
Snakify a String
A snakified string looks like this:

```<code>T AnE eOf ifi ing
h s x l A k e r
isI amp Sna dSt
</code>```

Your Task

Take a string s and a size n, then output the snakified string. The inputs ThisIsAnExampleOfaSnakifiedString and 3 would produce the example above.

Specifications

s will only contain ASCII characters between code points 33 and 126 inclusive (no spaces or newlines).

s will be between 1 and 100 characters long.

n is an integer representing the size of each output string segment. Each line of characters (up/down or left/right) that make up the curves in the "snake" is n characters long. See the test cases for examples.

n will be between 3 and 10 inclusive.

The output string always starts pointing downwards.

Trailing spaces on each line are allowed.

Trailing newlines at the end of the output are also allowed.

Leading spaces are not allowed.

code-golf means shortest code in bytes wins.

Test Cases

```<code>a 3

a

----------

Hello,World! 3

H Wor
e , l
llo d!

----------

ProgrammingPuzzlesAndCodeGolf 4

P  ngPu  Code
r  i  z  d  G
o  m  z  n  o
gram  lesA  lf

----------

IHopeYourProgramWorksForInputStringsWhichAre100CharactersLongBecauseThisTestCaseWillFailIfItDoesNot. 5

I   gramW   tStri   100Ch   gBeca   CaseW   DoesN
H   o   o   u   n   e   a   n   u   t   i   t   o
o   r   r   p   g   r   r   o   s   s   l   I   t
p   P   k   n   s   A   a   L   e   e   l   f   .
eYour   sForI   Which   cters   ThisT   FailI

----------

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~ 10

!        <=>?@ABCDE        `abcdefghi
"        ;        F        _        j
#        :        G        ^        k
$        9        H        ]        l
%        8        I        \        m
&        7        J        [        n
'        6        K        Z        o        ~
(        5        L        Y        p        }
)        4        M        X        q        |
*+,-./0123        NOPQRSTUVW        rstuvwxyz{
</code>```

答:
C, 138 bytes

```<code>char*h[]={"\e[B\e[D","","\e[A\e[D",""},t[999];i;main(n){system("clear");for(scanf("%s%d",t,&n),--n;t[i];++i)printf("%c%s",t[i],h[i/n%4]);}
</code>```

This uses ANSI escapes. Works in linux terminal.

Ungolfed:

```<code>char*h[]={"\e[B\e[D","","\e[A\e[D",""},
    /* cursor movement - h[0] moves the cursor one down and one left,
    h[2] moves the cursor one up and one left. */
t[999];i;
main(n){
    system("clear");
    for(scanf("%s%d",t,&n),--n;t[i];++i)
        printf("%c%s",t[i],h[i/n%4]);
}
</code>```

Title: Snakify a String

Body:
C, 138 bytes

char*h[]={"\e[B\e[D","","\e[A\e[D",""},t[999];i;main(n){system("clear");for(scanf("%s%d",t,&n),--n;t[i];++i)printf("%c%s",t[i],h[i/n%4]);}

This uses ANSI escapes. Works in linux terminal.

Ungolfed:

char*h[]={"\e[B\e[D","","\e[A\e[D",""},
    /* cursor movement - h[0] moves the cursor one down and one left,
    h[2] moves the cursor one up and one left. */
t[999];i;
main(n){
    system("clear");
    for(scanf("%s%d",t,&n),--n;t[i];++i)
        printf("%c%s",t[i],h[i/n%4]);
}


问:
Create your own Gym Exercise and follow along irl
Let's do an exercise together, shall we? Simply follow the movements of the stickman. What stickman you ask? The one we are about to create!

Make a program that allows for a string-input only containing integers, and outputs the following nine stick-figures with an interval in between:

```<code> @  \@/ |@_ _@| @/   \@ |@| \@\ /@/
/|\  |   |   |   |   |   |   |   |
/ \ / \ / \ / \ /|\ /|\ / \ / \ / \
 1   2   3   4   5   6   7   8   9
</code>```

This will result in an animation, which we can then follow along irl.

Here is an example of the output created when the input is "123245762":

Some rules:

This is tagged code-golf, so shortest code in bytes wins

The input is a string only containing [1-9] with a length potentially varying from 2 to 100 characters **

The animation must be an endless loop

The interval time must be 750 ms

No duplicated adjacent integers are allowed in the input - this also includes the first and last integers of the input (see the seventh and eight test cases) **

** How it fails in case of invalid input is completely up to you, but it should be clear that it fails. It could throw an exception; simply return at the start; reboot your PC; delete it's own source-/compile-code; hack into the military and shoots a rocket to the compiling GPS-location. Your call. EDIT: It's not allowed to simply display a (correct / half) animation as failing input. It must be clear that something is wrong with the input for the failing testcases below. Thanks to @user81655 for bringing this to my attention, hence the edit.

Subrules:

Please post the gif for one of the (non-failing) test cases below in your answer, which you can make very easy with the following program: screentogif.codeplex.com (Any other gif-creating program is also allowed of course.)

Bonus points if you actually do the exercise alongside the animation irl when it's finished. ;)

Correct testcases:

123245762

65

121234346565879879132418791576

Failing testcases:

2 // Too few characters

7282918274959292747383785189478174826894958127476192947512897571298593912374857471978269591928974518735891891723878 // Too much characters

1232405762 // Containing invalid input (the 0)

112212 // Fails because of the 11 and 22 present

1232457621 // Fails because of the starting and leading 1

答:
Pyth, 114 bytes

```<code>.V0j@cv.Z"xÚí» À0DW¹NM@+Üñ\">íÂA¸êÄÓw»`3±2þ&'NövfAé8é~0,p'ÆìÞúr_'¥,d!YÈBíéqs"3tv@zb.d.75
</code>```

Try it here

(No pausing, non-infinite)

Body:
Pyth, 114 bytes

.V0j@cv.Z"xÚí» À0DW¹NM@+Üñ\">íÂA¸êÄÓw»`3±2þ&'NövfAé8é~0,p'ÆìÞúr_'¥,d!YÈBíéqs"3tv@zb.d.75

Try it here

(No pausing, non-infinite)


问:
Simulate Cookie Clicker

Intro

Cookie Clicker is a popular browser game in which the player has to click a cookie to gain... Cookies!

The player can buy upgrades which will automatically farm cookies for you. Sometimes a golden cookie will enter the screen, which will provide a bonus upon clicking it.

Challenge

This challenge will not be about clicking a cookie. It will be about automating the cookie.

The player starts off with 0 cookies, and will gain cookies on every game tick. The goal of the challenge is to calculate the amount of ticks it takes to reach a certain amount of cookies.

Input

Your program or function will have 2 input parameters.

The first parameter is the goal, which is simply the amount of cookies you need to finish.

The second parameter is a list of upgrades the player can buy. The format of the list will be explained in further detail in the paragraph Upgrades - List format.

For sake of convenience, you may assume that the input is always in the correct format. You don't have to deal with bad input.

Output

Your program or function should print or return the amount of ticks it takes before the player has reached the goal amount of cookies.

Upgrades

Every game-tick will earn the player some cookies. These are based on the upgrades of the player. If a player has no upgrades, the cookie production will be 1 cookie per tick. Once the player reaches enough cookies to buy an upgrade, he will automatically purchase it (before the next tick happens).

When the player has enough cookies after a tick to buy multiple upgrades, he will attempt to buy all possible upgrades, starting with the cheapest.

List format

[[100, '+', 1],[1000, '*', 2],[10000, '*', 10]

This is an unordered list which contains 3 upgrades. The first parameter of an upgrade is the price of the upgrade. The second parameter is the effect of the upgrade and the third parameter is the effect amount.

The input list will contain no upgrades with the same price.

The only possible effects are + (addition) and * (multiplication).

[100, '+', 1] means that the cookie production will increase by 1 at a cost of 100 cookies.

Golden cookies

Every 1000 ticks, the player will find a golden cookie. When this happens, the player will receive a bonus of 10% (rounded down) of his cookie amount. If the player has 1000 cookies at tick 1000, he will have 1100 cookies before tick 1001 happens.

Order of events

After every tick a couple of actions can occur.

If the player has enough cookies, end the game

If the player has enough cookies, purchase upgrades

If the amount of ticks is a multiple of 1000; Golden cookie

Example input + output

```<code>(0    , [])                              -> 0
(100  , [])                              -> 100
(1    , [[10, '+', 1]])                  -> 1
(10   , [[10, '+', 1]])                  -> 10
(1000 , [[100, '*', 2], [10, '+', 1]])   -> 310
(10000, [[10, '+', 1], [100, '*', 2]])   -> 2263
</code>```

Detailed example

```<code>(1000 , [[100, '*', 2], [10, '+', 1]])   -> 310

ticks: 8, cookies: 8, production: 1
ticks: 9, cookies: 9, production: 1
ticks: 10, cookies: 0, production: 2
ticks: 11, cookies: 2, production: 2
ticks: 58, cookies: 96, production: 2
ticks: 59, cookies: 98, production: 2
ticks: 60, cookies: 0, production: 4
ticks: 61, cookies: 4, production: 4
ticks: 70, cookies: 40, production: 4
ticks: 80, cookies: 80, production: 4
ticks: 260, cookies: 800, production: 4
ticks: 280, cookies: 880, production: 4
ticks: 300, cookies: 960, production: 4
ticks: 309, cookies: 996, production: 4
ticks: 310, cookies: 1000, production: 4 -- done
</code>```

Disclaimer

I will not be held responsible for any addictions caused by this challenge.

Happy golfing!

Title: Simulate Cookie Clicker

Body:
Python2, 177 Bytes

def f(g,u):
 c=t=0;r=1
 while c<g:
    c+=r;t+=1
    if c>=g:break
    u=sorted(u)
    if u[0][0]<=c:
        h=u.pop(0);c-=h[0]
        exec('r'+h[1]+'=h[2]')
    if t%1000<1:c+=c//10
 return t


答:
JavaScript (Firefox), <s>86</s> 75 bytes

a=>a.concat(a[0]).slice(1).map((_,i)=>[for(v of a)if(n=v[i--])n].reverse())
```

Saved 11 bytes thanks to @Neil!

Works in Firefox 30+. Takes an array of arrays.

Title: Diamondize a Matrix

Body:
JavaScript (Firefox), <s>86</s> 75 bytes

a=>a.concat(a[0]).slice(1).map((_,i)=>[for(v of a)if(n=v[i--])n].reverse())

Saved 11 bytes thanks to @Neil!

Works in Firefox 30+. Takes an array of arrays.


答:
Mathematica, 252 bytes

```<code>i=0;Dynamic[Uncompress["1:eJxTTMoPSuNkYGAoZgESPpnFJcFCQIaCQ4yBoZF+TUwMmFaIiQELx8Q46IMEFGrgwoJA4RqHeCyi8Q41aKICQFEUA2qg5gIlHdCEIeaimyAMcQTEWWj26aO7DQDaqDEh"][[FromDigits[#~StringTake~{i=i~Mod~StringLength@#+1}]]],UpdateInterval->3/4,TrackedSymbols->{}]&
</code>```

Would be nice if someone could create a GIF. Run in a notebook.

Title: Create your own Gym Exercise and follow along irl

Body:
Let's do an exercise together, shall we? Simply follow the movements of the stickman. What stickman you ask? The one we are about to create!

Make a program that allows for a string-input only containing integers, and outputs the following nine stick-figures with an interval in between:

```<code> @  \@/ |@_ _@| @/   \@ |@| \@\ /@/
/|\  |   |   |   |   |   |   |   |
/ \ / \ / \ / \ /|\ /|\ / \ / \ / \
 1   2   3   4   5   6   7   8   9
</code>```

This will result in an animation, which we can then follow along irl.

Here is an example of the output created when the input is "123245762":

Some rules:

This is tagged code-golf, so shortest code in bytes wins

The input is a string only containing [1-9] with a length potentially varying from 2 to 100 characters **

The animation must be an endless loop

The interval time must be 750 ms

No duplicated adjacent integers are allowed in the input - this also includes the first and last integers of the input (see the seventh and eight test cases) **

** How it fails in case of invalid input is completely up to you, but it should be clear that it fails. It could throw an exception; simply return at the start; reboot your PC; delete it's own source-/compile-code; hack into the military and shoots a rocket to the compiling GPS-location. Your call. EDIT: It's not allowed to simply display a (correct / half) animation as failing input. It must be clear that something is wrong with the input for the failing testcases below. Thanks to @user81655 for bringing this to my attention, hence the edit.

Subrules:

Please post the gif for one of the (non-failing) test cases below in your answer, which you can make very easy with the following program: screentogif.codeplex.com (Any other gif-creating program is also allowed of course.)

Bonus points if you actually do the exercise alongside the animation irl when it's finished. ;)

Correct testcases:

123245762

65

121234346565879879132418791576

Failing testcases:

2 // Too few characters

7282918274959292747383785189478174826894958127476192947512897571298593912374857471978269591928974518735891891723878 // Too much characters

1232405762 // Containing invalid input (the 0)

112212 // Fails because of the 11 and 22 present

1232457621 // Fails because of the starting and leading 1

Body:
Mathematica, 252 bytes

i=0;Dynamic[Uncompress["1:eJxTTMoPSuNkYGAoZgESPpnFJcFCQIaCQ4yBoZF+TUwMmFaIiQELx8Q46IMEFGrgwoJA4RqHeCyi8Q41aKICQFEUA2qg5gIlHdCEIeaimyAMcQTEWWj26aO7DQDaqDEh"][[FromDigits[#~StringTake~{i=i~Mod~StringLength@#+1}]]],UpdateInterval->3/4,TrackedSymbols->{}]&

Would be nice if someone could create a GIF. Run in a notebook.


Body:
JavaScript (ES6), 210 bytes

s=>setInterval(_=>(c=console).clear(i=0)&c.log(`, @
/|\\
/ \\,\\@/
 |
/ \\,|@_
 |
/ \\,_@|
 |
/ \\,@/
 |
/|\\, \\@
 |
/|\\,|@|
 |
/ \\,\\@\\
 |
/ \\,/@/
 |
/ \\`.split`,`[s[i++%s.length]]),750)


问:
Where will the cat go? (orbital mechanics)
A nearly massless cat is dropped in space (don't worry, with a space suit and everything) at the point (x, y, z) with velocity (vx, vy, vz). There is an fixed, infinitely dense planet (with volume of 0) at the point (0, 0, 0) and it attracts objects at distance r with acceleration 1/r^2. According to Newtonian gravity, where does the object go after time t?

Nearly massless in this case means you are outputting the value of lim (mass --> 0) <position of cat>. The mass is affected by the planet's gravity, but the planet is not affected by the cat's gravity. In other words, the central body is fixed.

This is somewhat similar to Code Golf: What is the spaceship's fate? [floating point version], but this is different because it is measuring accuracy.

You may implement a solution based on a simulation, which must run in less than 3 seconds, OR you may implement a program that gives an exact value (must also run in less than 3 seconds). See scoring details below. If you implement a simulation, it does not have to be exact, but your score will be lower because of the inaccuracy.

Input: x y z vx vy vz t, not necessarily integers representing x, y, z coordinates, velocity in the x, y, and z directions and time respectively. It is guaranteed that the speed of the cat is strictly less than the escape velocity at that altitude. Input may be taken from anywhere, including parameters to a function. The program must run in less than three seconds on my laptop for t < 2^30, which means, if you are running a simulation, you must adjust your timestep accordingly. If you are planning on hitting the 3 second limit for every test case, make sure there is a tuneable parameter that can make it more accurate/less accurate for speed gains, so that I can make it run in three seconds on my computer.

Output: x y z, the position after time t.

Since the two-body problem can be solved perfectly, it is in theory possible to get a perfect, correct answer.

Scoring: For any test case, the error is defined to be the distance between your output and the "true" output. The true output is defined to be the one that the test case snippet generates. If the error is less than 10^(-8), the error is rounded down to zero. Your score is the average error on 100 (or more) random test cases. If you write a perfectly accurate answer, you should get a score of 0; lowest score wins, and ties will be broken by code length.

Test cases:

```<code>1 0 0 0 -1 0 1000000000 --> 0.83789 -0.54584 0
</code>```

In this case , the orbit is perfectly circular with period 2*pi, so after circling 159154943 times, the cat ends up at approximately (0.83789, -0.54584). This isn't a test case your code will be tested on; if you submit a perfectly accurate answer, however, you may want to test it on this.

The snippet below generates random additional test cases and will be used to judge submissions; let me know if there is a bug with this:

var generateTestCase = function() {
  var periapsis = Math.random() * 10 + 0.5;
  var circularVelocity = Math.sqrt(1 / periapsis);
  var escapeVelocity = Math.sqrt(2) * circularVelocity;
  var v = Math.random() * (escapeVelocity - circularVelocity) + circularVelocity;
  var a = 1 / (2 / periapsis - v * v);
  var ecc = (a - periapsis) * (2 / periapsis - v * v);
  var semiLatus = a * (1 - ecc * ecc);
  var angM = Math.sqrt(semiLatus);
  var period = 2 * Math.PI * a * Math.sqrt(a);
  var getCoordAfterTime = function(t) {
    var meanAnomaly = t / (a * Math.sqrt(a));

    function E(e, M, n) {
      var result = M;
      for (var k = 1; k <= n; ++k) result = M + e * Math.sin(result);
      return result;
    }
    var eAnomaly = E(ecc, meanAnomaly, 10000);
    var radius = a * (1 - ecc * Math.cos(eAnomaly));
    var theta = Math.acos((semiLatus / radius - 1) / ecc);
    if (t > period / 2) theta = -theta;
    return {
      r: radius,
      t: theta,
      x: radius * Math.cos(theta),
      y: radius * Math.sin(theta)
    };
  }
  var t1 = Math.random() * period;
  var p1 = getCoordAfterTime(t1);
  var v1 = Math.sqrt(2 / p1.r - 1 / a);
  var velocityAngle1 = Math.acos(angM / p1.r / v1);
  var horizontal = Math.atan(-p1.x / p1.y);
  if (t1 < period / 2) {
    velocityAngle1 *= -1;
    horizontal += Math.PI;
  }
  var vx = v1 * Math.cos(horizontal + velocityAngle1);
  var vy = v1 * Math.sin(horizontal + velocityAngle1);
  var t2 = Math.random() * period;
  var p2 = getCoordAfterTime(t2);
  var v2 = Math.sqrt(2 / p2.r - 1 / a);
  var dummyPeriods = Math.floor(Math.random() * 100) * period + period;
  var rotate = function(x, y, z, th1, th2, th3) {
    var x1 = x * Math.cos(th1) - y * Math.sin(th1);
    var y1 = x * Math.sin(th1) + y * Math.cos(th1);
    var z1 = z;
    var x2 = x1 * Math.cos(th2) + z1 * Math.sin(th2);
    var y2 = y1;
    var z2 = -x1 * Math.sin(th2) + z1 * Math.cos(th2);
    var x3 = x2;
    var y3 = y2 * Math.cos(th3) - z2 * Math.sin(th3);
    var z3 = y2 * Math.sin(th3) + z2 * Math.cos(th3);
    return [x3, y3, z3];
  }
  var th1 = Math.random() * Math.PI * 2,
    th2 = Math.random() * 2 * Math.PI,
    th3 = Math.random() * 2 * Math.PI;
  var input = rotate(p1.x, p1.y, 0, th1, th2, th3)
    .concat(rotate(vx, vy, 0, th1, th2, th3))
    .concat([t2 - t1 + dummyPeriods]);
  var output = rotate(p2.x, p2.y, 0, th1, th2, th3);
  return input.join(" ") + " --> " + output.join(" ");
}
var $ = function(x) {
  return document.querySelector(x);
}
$("#onetestcase").onclick = function() {
  $("#result1").innerHTML = generateTestCase();
}
$("#huntestcases").onclick = function() {
  var result = [];
  for (var i = 0; i < 100; ++i) result.push(generateTestCase());
  $("#result100").innerHTML = result.join("\n");
}```
<pre class="snippet-code-html lang-html prettyprint-override"><code><div>
  <button id="onetestcase">Generate</button>
  <pre id="result1">```

  <button id="huntestcases">Generate 100 test cases</button>

```
</div></code>```

答:
Javascript

This is just to get the ball rolling, since no one seems to be posting answers. Here is a very naive, simple way that can be improved a lot:

function simulate(x, y, z, vx, vy, vz, t) {
  var loops = 1884955; // tune this parameter
  var timestep = t / loops;
  for (var i = 0; i < t; i += timestep) {
    var distanceSq = x*x + y*y + z*z; // distance squared from origin
    var distance = Math.sqrt(distanceSq);
    var forceMag = 1/distanceSq; // get the force of gravity
    var forceX = -x / distance * forceMag;
    var forceY = -y / distance * forceMag;
    var forceZ = -z / distance * forceMag;
    vx += forceX * timestep;
    vy += forceY * timestep;
    vz += forceZ * timestep;
    x += vx * timestep;
    y += vy * timestep;
    z += vz * timestep;
  }
  return [x, y, z];
}
```
<p>Testing:</p>
<pre class="lang-js prettyprint-override"><code>simulate(1, 0, 0, 0, -1, 0, Math.PI*2) --> [0.9999999999889703, -0.0000033332840909716455, 0]
</code>```

Hey, that's pretty good. It has an error of about 3.333*10^(-6) which isn't enough for it to be rounded down... it's close.

Just for fun:

console.log(simulate(1, 0, 0, 0, -1, 0, 1000000000))
--> [-530516643639.4616, -1000000000.0066016, 0]
```
<p>Oh well; so it's not the best.</p>
<p>And on a random test case from the generator:</p>
<pre class="lang-js prettyprint-override"><code>simulate(-1.1740058273269156,8.413493259550673,0.41996042044140003,0.150014367067652,-0.09438816345868332,0.37294941703455975,7999.348650387233)
-->    [-4.528366392498373, 6.780385554803544, -9.547824236472668]
Actual:[-4.452695438880813, 6.932249293597744, -9.272924876103785]
</code>```

With an error of only approximately 0.32305!

This can be improved a lot by using something like Verlet integration or some fancy algorithm. In fact, those algorithms may even get perfect scores, despite being simulations.

Body:
Javascript

This is just to get the ball rolling, since no one seems to be posting answers. Here is a very naive, simple way that can be improved a lot:

function simulate(x, y, z, vx, vy, vz, t) {
  var loops = 1884955; // tune this parameter
  var timestep = t / loops;
  for (var i = 0; i < t; i += timestep) {
    var distanceSq = x*x + y*y + z*z; // distance squared from origin
    var distance = Math.sqrt(distanceSq);
    var forceMag = 1/distanceSq; // get the force of gravity
    var forceX = -x / distance * forceMag;
    var forceY = -y / distance * forceMag;
    var forceZ = -z / distance * forceMag;
    vx += forceX * timestep;
    vy += forceY * timestep;
    vz += forceZ * timestep;
    x += vx * timestep;
    y += vy * timestep;
    z += vz * timestep;
  }
  return [x, y, z];
}

Testing:

simulate(1, 0, 0, 0, -1, 0, Math.PI*2) --> [0.9999999999889703, -0.0000033332840909716455, 0]

Hey, that's pretty good. It has an error of about 3.333*10^(-6) which isn't enough for it to be rounded down... it's close.

Just for fun:

console.log(simulate(1, 0, 0, 0, -1, 0, 1000000000))
--> [-530516643639.4616, -1000000000.0066016, 0]

Oh well; so it's not the best.

And on a random test case from the generator:

simulate(-1.1740058273269156,8.413493259550673,0.41996042044140003,0.150014367067652,-0.09438816345868332,0.37294941703455975,7999.348650387233)
-->    [-4.528366392498373, 6.780385554803544, -9.547824236472668]
Actual:[-4.452695438880813, 6.932249293597744, -9.272924876103785]

With an error of only approximately 0.32305!

This can be improved a lot by using something like Verlet integration or some fancy algorithm. In fact, those algorithms may even get perfect scores, despite being simulations.


答:
Java 8, <s>663</s> <s>636</s> <s>634</s> <S>631</s> <s>596</s> <s>355</s> <s>354</s> 343 bytes

Just for the lols I tried to make the program in Java. <s>Admittedly I'm pretty bad at golfing and regexes, so it can most likely be golfed (a lot?) more. Nevertheless, here it is in Java 7.</s>

Now almost two years later and I almost halved the code in Java 8. Why did I ever made those rules about validating the input and requiring a full program, though... >.> I hate my past self now..

interface M{static void main(String[]a)throws Exception{for(;a[0].matches("[1-9]{2,100}")&!a[0].matches(".*(.)\\2.*|(.).*\\1");)for(int c:a[0].getBytes()){c-=48;System.out.printf("%s%n%s%n%s%n",c<2?" @ ":c<3?"\\@/":c<4?"|@_":c<5?"_@|":c<6?"@/ ":c<7?" \\@":c<8?"|@|":c<9?"\\@\\":"/@/",c<2?"/|\\":" | ",c%7>4?"/|\\":"/ \\");Thread.sleep(750);}}}
```
<p>-11 bytes thanks to <em>@ceilingcat</em>.</p>
<p><strong>Explanation:</strong></p>
<p>Try it online. (After it has timed out after 60 sec.)</p>
<pre class="lang-java prettyprint-override"><code>interface M{                     // Class
  static void main(String[]a)    //  Mandatory main-method
      throws Exception{          //    Required throws for the `Thread.sleep`
    for(;a[0].matches("[1-9]{2,100}")
                                 //   Validate 2-100 integers only containing 1-9
         &!a[0].matches(".*(.)\\2.*|(.).*\\1")
                                 //   Validate no adjacent duplicated char (with wrap-around)
        ;)                       //   And start an infinite loop if it's valid:
      for(int c:a[0].getBytes()){//    Inner loop over the characters of the input:
        c-=48;                   //     Convert the character-code to an integer
        System.out.printf("%s%n%s%n%s%n",
                                 //     Print:
          c<2?" @ ":c<3?"\\@/":c<4?"|@_":c<5?"_@|":c<6?"@/ ":c<7?" \\@":c<8?"|@|":c<9?"\\@\\":"/@/",
                                 //      The top part of the stick figure
          c<2?"/|\\":" | ",      //      The middle part of the stick figure
          c%7>4?"/|\\":"/ \\");  //      The bottom part of the stick figure
        Thread.sleep(750);}}}    //     Sleep 750 ms
</code>```

Gif:

(Note: Old gif, since it's clearly jdk1.8+ now.)

Body:
Java 8, <s>663</s> <s>636</s> <s>634</s> <S>631</s> <s>596</s> <s>355</s> <s>354</s> 343 bytes

Just for the lols I tried to make the program in Java. <s>Admittedly I'm pretty bad at golfing and regexes, so it can most likely be golfed (a lot?) more. Nevertheless, here it is in Java 7.</s>

Now almost two years later and I almost halved the code in Java 8. Why did I ever made those rules about validating the input and requiring a full program, though... >.> I hate my past self now..

interface M{static void main(String[]a)throws Exception{for(;a[0].matches("[1-9]{2,100}")&!a[0].matches(".*(.)\\2.*|(.).*\\1");)for(int c:a[0].getBytes()){c-=48;System.out.printf("%s%n%s%n%s%n",c<2?" @ ":c<3?"\\@/":c<4?"|@_":c<5?"_@|":c<6?"@/ ":c<7?" \\@":c<8?"|@|":c<9?"\\@\\":"/@/",c<2?"/|\\":" | ",c%7>4?"/|\\":"/ \\");Thread.sleep(750);}}}

-11 bytes thanks to @ceilingcat.

Explanation:

Try it online. (After it has timed out after 60 sec.)

interface M{                     // Class
  static void main(String[]a)    //  Mandatory main-method
      throws Exception{          //    Required throws for the `Thread.sleep`
    for(;a[0].matches("[1-9]{2,100}")
                                 //   Validate 2-100 integers only containing 1-9
         &!a[0].matches(".*(.)\\2.*|(.).*\\1")
                                 //   Validate no adjacent duplicated char (with wrap-around)
        ;)                       //   And start an infinite loop if it's valid:
      for(int c:a[0].getBytes()){//    Inner loop over the characters of the input:
        c-=48;                   //     Convert the character-code to an integer
        System.out.printf("%s%n%s%n%s%n",
                                 //     Print:
          c<2?" @ ":c<3?"\\@/":c<4?"|@_":c<5?"_@|":c<6?"@/ ":c<7?" \\@":c<8?"|@|":c<9?"\\@\\":"/@/",
                                 //      The top part of the stick figure
          c<2?"/|\\":" | ",      //      The middle part of the stick figure
          c%7>4?"/|\\":"/ \\");  //      The bottom part of the stick figure
        Thread.sleep(750);}}}    //     Sleep 750 ms

Gif:

(Note: Old gif, since it's clearly jdk1.8+ now.)


答:
Mathematica, <s>58</s> 56 bytes

```<code>a=Length;Reverse@#~Diagonal~b~Table~{b,1-a@#,a@#&@@#-1}&
</code>```

Anonymous function, takes nested arrays.

Title: Diamondize a Matrix

Body:
Mathematica, <s>58</s> 56 bytes

a=Length;Reverse@#~Diagonal~b~Table~{b,1-a@#,a@#&@@#-1}&

Anonymous function, takes nested arrays.


问:
Map string to Hilbert curve
Let's map some strings to 2d space, fractal style. Your task is to compute a Hilbert curve and lay a string along it.

Task

The task is to take the single-line input string, and lay it out along a Hilbert curve big enough to contain it, but no bigger. Try to make the byte count as low as possible; this is code-golf after all!

Conditions

Any gaps to be padded with whitespace, but padding is not required at the end of lines.

The start of the line should be in the top-left corner, and the end in the bottom-left.

You may create a program or function.

There may be some new test-cases appearing, so don't hardcode anything!

Bonuses

Note: Bonuses stack like this: -50% & -20% on 100B = -20% on 50B or -50% on 80B = 40B.

-50% If the input is a multi-line string, reverse the process to create the original input. Test cases for the bonus: just use the existing ones (including the bonus test cases!)

-20% If you strip all unnecessary whitespace from the output (e.g. at the end of a line).

-5% If you don't pollute the global namespace (you know what I mean!)

Test cases

```<code>abcdefghijklmn

adef
bchg
 nij
 mlk

The quick brown fox jumps over the lazy dog.

Thn f ju
 ewooxpm
qckr rs
ui btevo
    hlaz
    e  y
      do
      .g
</code>```

And for the whitespace-stripping bonus:

```<code>No  hitespac  her

Noher

hesc
itpa
</code>```

Leaderboard

/* Configuration */

var QUESTION_ID = 66958; // Obtain this from the url
// It will be like https://XYZ.stackexchange.com/questions/QUESTION_ID/... on any question page
var ANSWER_FILTER = "!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe";
var COMMENT_FILTER = "!)Q2B_A2kjfAiU78X(md6BoYk";
var OVERRIDE_USER = 43394; // This should be the user ID of the challenge author.

/* App */

var answers = [], answers_hash, answer_ids, answer_page = 1, more_answers = true, comment_page;

function answersUrl(index) {
  return "https://api.stackexchange.com/2.2/questions/" +  QUESTION_ID + "/answers?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + ANSWER_FILTER;
}

function commentUrl(index, answers) {
  return "https://api.stackexchange.com/2.2/answers/" + answers.join(';') + "/comments?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + COMMENT_FILTER;
}

function getAnswers() {
  jQuery.ajax({
    url: answersUrl(answer_page++),
    method: "get",
    dataType: "jsonp",
    crossDomain: true,
    success: function (data) {
      answers.push.apply(answers, data.items);
      answers_hash = [];
      answer_ids = [];
      data.items.forEach(function(a) {
        a.comments = [];
        var id = +a.share_link.match(/\d+/);
        answer_ids.push(id);
        answers_hash[id] = a;
      });
      if (!data.has_more) more_answers = false;
      comment_page = 1;
      getComments();
    }
  });
}

function getComments() {
  jQuery.ajax({
    url: commentUrl(comment_page++, answer_ids),
    method: "get",
    dataType: "jsonp",
    crossDomain: true,
    success: function (data) {
      data.items.forEach(function(c) {
        if (c.owner.user_id === OVERRIDE_USER)
          answers_hash[c.post_id].comments.push(c);
      });
      if (data.has_more) getComments();
      else if (more_answers) getAnswers();
      else process();
    }
  });
}

getAnswers();

var SCORE_REG = /<h\d>\s*([^\n,]*[^\s,]),.*?(\d+(?:[.]\d+)?)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/;

var OVERRIDE_REG = /^Override\s*header:\s*/i;

function getAuthorName(a) {
  return a.owner.display_name;
}

function process() {
  var valid = [];

  answers.forEach(function(a) {
    var body = a.body;
    a.comments.forEach(function(c) {
      if(OVERRIDE_REG.test(c.body))
        body = '
' + c.body.replace(OVERRIDE_REG, '') + '
';
    });

    var match = body.match(SCORE_REG);
    if (match)
      valid.push({
        user: getAuthorName(a),
        size: +match[2],
        language: match[1],
        link: a.share_link,
      });

  });

  valid.sort(function (a, b) {
    var aB = a.size,
        bB = b.size;
    return aB - bB
  });

  var languages = {};
  var place = 1;
  var lastSize = null;
  var lastPlace = 1;
  valid.forEach(function (a) {
    if (a.size != lastSize)
      lastPlace = place;
    lastSize = a.size;
    ++place;

    var answer = jQuery("#answer-template").html();
    answer = answer.replace("{{PLACE}}", lastPlace + ".")
                   .replace("{{NAME}}", a.user)
                   .replace("{{LANGUAGE}}", a.language)
                   .replace("{{SIZE}}", a.size)
                   .replace("{{LINK}}", a.link);
    answer = jQuery(answer);
    jQuery("#answers").append(answer);

    var lang = a.language;
    if (/ b.lang) return 1;
    if (a.lang < b.lang) return -1;
    return 0;
  });

  for (var i = 0; i < langs.length; ++i)
  {
    var language = jQuery("#language-template").html();
    var lang = langs[i];
    language = language.replace("{{LANGUAGE}}", lang.lang)
                       .replace("{{NAME}}", lang.user)
                       .replace("{{SIZE}}", lang.size)
                       .replace("{{LINK}}", lang.link);
    language = jQuery(language);
    jQuery("#languages").append(language);
  }

}```
<pre class="snippet-code-css lang-css prettyprint-override"><code>body { text-align: left !important}

#answer-list {
  padding: 10px;
  width: 290px;
  float: left;
}

#language-list {
  padding: 10px;
  width: 290px;
  float: left;
}

table thead {
  font-weight: bold;
}

table td {
  padding: 5px;
}</code>```

<link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b">

Leaderboard

	Author	Language	Size

Winners by Language

Language	User	Score

{{PLACE}}	{{NAME}}	{{LANGUAGE}}	{{SIZE}}	Link

{{LANGUAGE}}	{{NAME}}	{{SIZE}}	Link

```
</div>
</div>
</p>

<p>To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:</p>

```# Language Name, N bytes
```

<p>where <code>N</code> is the size of your submission. If you improve your score, you <em>can</em> keep old scores in the headline, by striking them through. For instance:</p>

```# Ruby, <s>104</s> <s>101</s> 96 bytes
```

<p>If there you want to include multiple numbers in your header (e.g. because your score is the sum of two files or you want to list interpreter flag penalties separately), make sure that the actual score is the <em>last</em> number in the header:</p>

```# Perl, 43 + 2 (-p flag) = 45 bytes
```

<p>You can also make the language name a link which will then show up in the leaderboard snippet:</p>

```# [><>](http://esolangs.org/wiki/Fish), 121 bytes
```

答:
JavaScript (ES6), 227 - 20% : 181.6 bytes

m=>{for(n=1<<((33-Math.clz32(m.length-1))/2),t='',y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}
```

<p><em>Trying to get the 5% bonus</em></p>

<pre class="lang-js prettyprint-override"><code>m=>{for(var n=1<<((33-Math.clz32(m.length-1))/2),t='',x,y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(var p,q,u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}
</code>```

241 * 0.8 *0.95 : 183.16 bigger

Less golfed

m=>
{
  // calc the size of the bounding square, clz32 is a bit shorter than ceil(log2()
  n = 1<<( (33-Math.clz32(m.length-1)) / 2);
  t = '';
  for(y = 0; y < n; y++)
  {
    for(x = 0 ; x < n; x++)
    {
      // for each position x,y inside the square
      // get the index postion in the hilbert curve
      // see https://en.wikipedia.org/wiki/Hilbert_curve (convert x,y to d)
      for(u=y, v=x, h=0, s=n; s >>= 1; )
      {
        h += s*s*(3 * !!(p = u & s) ^ !!(q = v & s));
        q || (p && (u = s+~u, v = s+~v),[u,v]=[v,u])
      }
      // add char at given index to output
      t += m[h]||' '; // blank if beyond the length of m
    }
    t += '\n'; // add newline add end line
  }
  return t.replace(/ +$/mg,'').trim() // to get the 20% bonus
}
```

<p><strong>Test</strong></p>

<p><div class="snippet" data-lang="js" data-hide="false">
<div class="snippet-code">
<pre class="snippet-code-js lang-js prettyprint-override"><code>F=m=>{for(n=1<<((33-Math.clz32(m.length-1))/2),t='',y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}

function Test() { O.textContent = F(I.value) }

Test()</code>```

#I { width: 90% }
#O { border: 1px solid #ccc}```
<pre class="snippet-code-html lang-html prettyprint-override"><code><input id=I oninput='Test()' value='The quick brown fox jumps over the lazy dog.'>
<pre id=O>``````

Title: Map string to Hilbert curve

Body:
JavaScript (ES6), 227 - 20% : 181.6 bytes

m=>{for(n=1<<((33-Math.clz32(m.length-1))/2),t='',y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}

Trying to get the 5% bonus

m=>{for(var n=1<<((33-Math.clz32(m.length-1))/2),t='',x,y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(var p,q,u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}

241 * 0.8 *0.95 : 183.16 bigger

Less golfed

m=>
{
  // calc the size of the bounding square, clz32 is a bit shorter than ceil(log2()
  n = 1<<( (33-Math.clz32(m.length-1)) / 2);
  t = '';
  for(y = 0; y < n; y++)
  {
    for(x = 0 ; x < n; x++)
    {
      // for each position x,y inside the square
      // get the index postion in the hilbert curve
      // see https://en.wikipedia.org/wiki/Hilbert_curve (convert x,y to d)
      for(u=y, v=x, h=0, s=n; s >>= 1; )
      {
        h += s*s*(3 * !!(p = u & s) ^ !!(q = v & s));
        q || (p && (u = s+~u, v = s+~v),[u,v]=[v,u])
      }
      // add char at given index to output
      t += m[h]||' '; // blank if beyond the length of m
    }
    t += '\n'; // add newline add end line
  }
  return t.replace(/ +$/mg,'').trim() // to get the 20% bonus
}

Test

F=m=>{for(n=1<<((33-Math.clz32(m.length-1))/2),t='',y=0;y<n;y++,t+=`
`)for(x=0;x<n;x++,t+=m[h]||' ')for(u=y,v=x,h=0,s=n;s>>=1;q||(p&&(u=s+~u,v=s+~v),[u,v]=[v,u]))h+=s*s*(3*!!(p=u&s)^!!(q=v&s));return t.replace(/ +$/mg,'').trim()}

function Test() { O.textContent = F(I.value) }

Test()

#I { width: 90% }
#O { border: 1px solid #ccc}

<input id=I oninput='Test()' value='The quick brown fox jumps over the lazy dog.'>
<pre id=O></pre>


答:
Reng v.3.3, 137 bytes, noncompeting

Try it here!

```<code>aií0#zl2,q!~1ø
:"]"eq!v:"}"eq!v:">"eq!v:")"eq!v)1z+#z
ve¤[2-2<       <       <     +1<
>]?v$$$zÀ0#z >ðq!vlqv¤l2%[1Ø
   \$2+)1z+#z/   ~n1/
</code>```

There's a bit more golfing to be done, but at least it works. I added a command ð to keep track of stacks after this challenge in order for this to be remotely possible/easily. I'll explain this in a bit, but it generally keeps track of all strings iterated over and looks for repeats; if there is a repeat, then the string is irreducible. Otherwise, the string will be reduced to the empty string/stack, and will output 1. Otherwise, no output will be produced.

Title: Are the brackets fully matched?

Body:
You must write a program or function that takes a string of brackets and outputs whether or not that string is fully matched. Your program should print a truthy or falsy value, and IO can be in any reasonable format.

Rules and definitions:

For the purpose of this challenge, a "bracket" is any of these characters: ()[]{}<>.

A pair of brackets is considered "matched" if the opening and closing brackets are in the right order and have no characters inside of them, such as

```<code>()
[]{}
</code>```

Or if every subelement inside of it is also matched.

```<code>[()()()()]
{<[]>}
(()())
</code>```

Subelements can also be nested several layers deep.

```<code>[(){<><>[()]}<>()]
<[{((()))}]>
</code>```

A string is considered "Fully matched" if and only if:

Every single character is a bracket,

Each pair of brackets has the correct opening and closing bracket and in the right order, and

Each bracket is matched.

You may assume the input will only contain printable ASCII.

Test IO

Here are some inputs that should return a truthy value:

```<code>()
[](){}<>
(((())))
({[<>]})
[{()<>()}[]]
[([]{})<{[()<()>]}()>{}]
</code>```

And here are some outputs that should return a falsy value:

```<code>(               Has no closing ')'
}{              Wrong order
(<)>            Each pair contains only half of a matched element
(()()foobar)    Contains invalid characters
[({}<>)>        The last bracket should be ']' instead of '>'
(((()))         Has 4 opening brackets, but only 3 closing brackets.
</code>```

As usual, this is code-golf, so standard loopholes apply, and shortest answer in bytes wins.

Body:
Reng v.3.3, 137 bytes, noncompeting

Try it here!

aií0#zl2,q!~1ø
:"]"eq!v:"}"eq!v:">"eq!v:")"eq!v)1z+#z
ve¤[2-2<       <       <     +1<
>]?v$$$zÀ0#z >ðq!vlqv¤l2%[1Ø
   \$2+)1z+#z/   ~n1/

There's a bit more golfing to be done, but at least it works. I added a command ð to keep track of stacks after this challenge in order for this to be remotely possible/easily. I'll explain this in a bit, but it generally keeps track of all strings iterated over and looks for repeats; if there is a repeat, then the string is irreducible. Otherwise, the string will be reduced to the empty string/stack, and will output 1. Otherwise, no output will be produced.


答:
Mathematica, 60 bytes

```<code>#&@@@#&/@GatherBy[Join@@MapIndexed[List,#,{2}],Tr@*Last]&
</code>```

where  is a Unicode character which Mathematica reads as the postfix \[Transpose] operator.

This is a bit longer than the other Mathematica solution but I figured I'd post it because it doesn't use the Diagonals built-in and uses a completely different approach.

Explanation

```<code>MapIndexed[List,#,{2}]
</code>```

This first transposes the matrix (such that the antidiagonals appear in the correct order if the matrix was flattened). Then we map List over the cells of the matrix together with the index, which turns each matrix element i into {i, {x, y}} where x and y are the coordinates of the element in the matrix.

```<code>Join@@...
</code>```

This flattens the outermost dimension, so that we now have a flat list of the matrix elements (with their coordinates) in column-major order.

```<code>GatherBy[..., Tr@*Last]
</code>```

This groups those elements by the sum of their coordinates. Note that antidiagonals are lines of constant x+y, so this does exactly the grouping we want. The order within each group is preserved. Now we just need to get rid of the coordinates again. This is done via the rather cryptic:

```<code>#&@@@#&/@...
</code>```

This maps the function #&@@@#& over each group, which itself applies #& to each element in the group, and # is simply the first argument, i.e. the original matrix element.

Title: Diamondize a Matrix

Body:
Mathematica, 60 bytes

#&@@@#&/@GatherBy[Join@@MapIndexed[List,#,{2}],Tr@*Last]&

where  is a Unicode character which Mathematica reads as the postfix \[Transpose] operator.

This is a bit longer than the other Mathematica solution but I figured I'd post it because it doesn't use the Diagonals built-in and uses a completely different approach.

Explanation

MapIndexed[List,#,{2}]

This first transposes the matrix (such that the antidiagonals appear in the correct order if the matrix was flattened). Then we map List over the cells of the matrix together with the index, which turns each matrix element i into {i, {x, y}} where x and y are the coordinates of the element in the matrix.

Join@@...

This flattens the outermost dimension, so that we now have a flat list of the matrix elements (with their coordinates) in column-major order.

GatherBy[..., Tr@*Last]

This groups those elements by the sum of their coordinates. Note that antidiagonals are lines of constant x+y, so this does exactly the grouping we want. The order within each group is preserved. Now we just need to get rid of the coordinates again. This is done via the rather cryptic:

#&@@@#&/@...

This maps the function #&@@@#& over each group, which itself applies #& to each element in the group, and # is simply the first argument, i.e. the original matrix element.


答:
CJam, 17 bytes

```<code>{eeSf.*::+W%zSf-}
</code>```

An unnamed block (function) which expects the matrix on the stack and replaces it with its antidiagonals.

Test it here.

This (found by Sp3000) works for the same byte count:

```<code>{_,,Sf*\.+W%zSf-}
</code>```

Explanation

This is best explained with an example. Consider the input:

```<code>[[0  1  2  3]
 [4  5  6  7]
 [8  9 10 11]]

ee    e# Enumerate matrix, turning each row [x ... z] into [i [x ... z]] where
      e# i is the vertical index from the top.

[[0 [0  1  2  3]]
 [1 [4  5  6  7]]
 [2 [8  9 10 11]]]

Sf.*  e# Replace each i with a string of i spaces.

[[""   [0  1  2  3]]
 [" "  [4  5  6  7]]
 ["  " [8  9 10 11]]]

::+   e# Prepend these strings to the rows.

[[0  1  2  3]
 ['  4  5  6  7]
 ['  '  8  9 10 11]]   e# Note that each '  corresponds to a space character.

W%    e# Reverse the rows.

[['  '  8  9 10 11]
 ['  4  5  6  7]
 [0  1  2  3]]

z     e# Zip/transpose.

[[ '  '  0]
 [ '  4  1]
 [ 8  5  2]
 [ 9  6  3]
 [10  7]
 [11]]

Sf-   e# Remove spaces from each row.

[[ 0]
 [ 4  1]
 [ 8  5  2]
 [ 9  6  3]
 [10  7]
 [11]]
</code>```

Title: Diamondize a Matrix

Body:
CJam, 17 bytes

{eeSf.*::+W%zSf-}

An unnamed block (function) which expects the matrix on the stack and replaces it with its antidiagonals.

Test it here.

This (found by Sp3000) works for the same byte count:

{_,,Sf*\.+W%zSf-}

Explanation

This is best explained with an example. Consider the input:

[[0  1  2  3]
 [4  5  6  7]
 [8  9 10 11]]

ee    e# Enumerate matrix, turning each row [x ... z] into [i [x ... z]] where
      e# i is the vertical index from the top.

[[0 [0  1  2  3]]
 [1 [4  5  6  7]]
 [2 [8  9 10 11]]]

Sf.*  e# Replace each i with a string of i spaces.

[[""   [0  1  2  3]]
 [" "  [4  5  6  7]]
 ["  " [8  9 10 11]]]

::+   e# Prepend these strings to the rows.

[[0  1  2  3]
 ['  4  5  6  7]
 ['  '  8  9 10 11]]   e# Note that each '  corresponds to a space character.

W%    e# Reverse the rows.

[['  '  8  9 10 11]
 ['  4  5  6  7]
 [0  1  2  3]]

z     e# Zip/transpose.

[[ '  '  0]
 [ '  4  1]
 [ 8  5  2]
 [ 9  6  3]
 [10  7]
 [11]]

Sf-   e# Remove spaces from each row.

[[ 0]
 [ 4  1]
 [ 8  5  2]
 [ 9  6  3]
 [10  7]
 [11]]


答:
Perl, <s>34</s> 33 bytes

Includes +2 for -lp

Run with input on STDIN:

```<code>./brackets.pl <<< "{<>()}"
</code>```

brackets.pl:

```<code>#!/usr/bin/perl -lp
s/\(\)|\[]|<>|{}//&&redo;$_=!$_
</code>```

Finds the first bracket pair without anything between them and removes it as long as there are any. Then checks if the final string is empty.

Title: Are the brackets fully matched?

Body:
Perl, <s>34</s> 33 bytes

Includes +2 for -lp

Run with input on STDIN:

./brackets.pl <<< "{<>()}"

brackets.pl:

#!/usr/bin/perl -lp
s/\(\)|\[]|<>|{}//&&redo;$_=!$_

Finds the first bracket pair without anything between them and removes it as long as there are any. Then checks if the final string is empty.