Przybliżenie specjalnego przypadku funkcji Riemanna Thety

27

Wyzwanie polega na napisaniu szybkiego kodu, który może wykonać trudną obliczeniowo nieskończoną sumę.

Wkład

nPrzez nmatrycę Pz pozycji całkowitych, które są mniejsze niż 100wartości bezwzględnej. Podczas testowania z przyjemnością dostarczam dane wejściowe do Twojego kodu w dowolnym rozsądnym formacie, jakiego chce Twój kod. Domyślnie będzie to jeden wiersz na wiersz matrycy, oddzielone spacją i zapewnione na standardowym wejściu.

Pbędzie pozytywnie określony, co oznacza, że ​​zawsze będzie symetryczny. Poza tym tak naprawdę nie musisz wiedzieć, co oznacza pozytywny zdecydowanie, aby odpowiedzieć na wyzwanie. Oznacza to jednak, że faktycznie będzie odpowiedź na sumę zdefiniowaną poniżej.

Musisz jednak wiedzieć, czym jest produkt macierz-wektor .

Wydajność

Twój kod powinien obliczyć nieskończoną sumę:

wprowadź opis zdjęcia tutaj

w granicach plus lub minus 0,0001 poprawnej odpowiedzi. Oto Zzestaw liczb całkowitych, podobnie jak Z^nwszystkie możliwe wektory z nelementami liczb całkowitych i ejest to słynna stała matematyczna, która w przybliżeniu wynosi 2,71828. Zauważ, że wartość wykładnika jest po prostu liczbą. Poniżej wyraźny przykład.

Jak to się ma do funkcji Riemanna Thety?

W notacji tego artykułu o aproksymacji funkcji Riemanna Thety próbujemy obliczyć wprowadź opis zdjęcia tutaj. Nasz problem jest szczególnym przypadkiem z co najmniej dwóch powodów.

  • Ustawiamy parametr początkowy wywoływany zw powiązanym papierze na 0.
  • Matrycę tworzymy Pw taki sposób, aby minimalny rozmiar wartości własnej wynosił 1. (Zobacz poniżej, jak tworzona jest macierz).

Przykłady

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

Output: 1.07551411208

Bardziej szczegółowo, zobaczmy tylko jeden termin w sumie dla tego P. Weźmy na przykład tylko jeden termin w sumie:

wprowadź opis zdjęcia tutaj

a x^T P x = 30. Zauważ, że e^(-30)to jest około 10^(-14)i dlatego jest mało prawdopodobne, aby było ważne dla uzyskania prawidłowej odpowiedzi do danej tolerancji. Przypomnij sobie, że nieskończona suma faktycznie wykorzysta każdy możliwy wektor długości 4, gdzie elementy są liczbami całkowitymi. Właśnie wybrałem jeden, aby podać wyraźny przykład.

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

Wynik

Przetestuję twój kod na losowo wybranych macierzach P o rosnącym rozmiarze.

Twój wynik jest po prostu największy, ndla którego otrzymuję poprawną odpowiedź w mniej niż 30 sekund, gdy uśrednia się ponad 5 przebiegów z losowo wybranymi matrycami Ptego rozmiaru.

A co z krawatem?

W przypadku remisu zwycięzcą zostanie ten, którego kod działa najszybciej, uśredniony w ciągu 5 przebiegów. Jeśli czasy te są równe, zwycięzcą jest pierwsza odpowiedź.

Jak powstanie losowe wejście?

  1. Niech M będzie losową macierzą m na n macierzą m <= n i pozycjami -1 lub 1. W Pythonie / numpy M = np.random.choice([0,1], size = (m,n))*2-1. W praktyce będzie ustawić mna około n/2.
  2. Niech P będzie macierzą tożsamości + M ^ T M. W Pythonie / numpy P =np.identity(n)+np.dot(M.T,M). Mamy teraz gwarancję, że Pjest pozytywnie określony, a wpisy są w odpowiednim zakresie.

Zauważ, że oznacza to, że wszystkie wartości własne P wynoszą co najmniej 1, co sprawia, że ​​problem jest potencjalnie łatwiejszy niż ogólny problem aproksymacji funkcji Riemanna Thety.

Języki i biblioteki

Możesz użyć dowolnego języka lub biblioteki, którą lubisz. Jednak do celów punktacji uruchomię Twój kod na moim komputerze, więc podaj jasne instrukcje, jak uruchomić go na Ubuntu.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja Ubuntu na 8-rdzeniowym procesorze AMD FX-8350 8 GB. Oznacza to również, że muszę być w stanie uruchomić Twój kod.


Wiodące odpowiedzi

  • n = 47w C ++ przez Ton Hospel
  • n = 8w Python autorstwa Maltysen
Glorfindel
źródło
Warto wspomnieć, że dodatnia określona macierz jest z definicji symetryczna.
2012rcampion
@ 2012rcampion Dzięki. Dodany.
Ok, może to głupie pytanie, ale ja patrzył na to przez wieki i nie mogę dowiedzieć się, w jaki sposób dostał xod [-1,0,2,1]. Czy możesz to rozwinąć? (Wskazówka: nie jestem guru matematyki)
wnnmaw
@wnnmaw Przepraszamy za zamieszanie. W tym przypadku suma ma jeden składnik na każdy możliwy wektor x długości 4. [-1,0,2,1] to tylko jeden, który wybrałem losowo, aby wyraźnie pokazać, jaki byłby termin w takim przypadku.
1
@Lembik Sposób, w jaki generujesz macierze SPD, oznacza, że ​​żadna pojedyncza wartość nie ma nigdy wartości bezwzględnej poniżej 1. Czy możemy skorzystać z tej wiedzy?
flawr

Odpowiedzi:

15

C ++

Nigdy więcej naiwnego podejścia. Oceniaj tylko wewnątrz elipsoidy.

Korzysta z bibliotek armadillo, ntl, gsl i pthread. Zainstaluj za pomocą

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

Skompiluj program używając czegoś takiego:

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

W niektórych systemach może być konieczne dodanie -lgslcblaspóźniej -lgsl.

Uruchom z rozmiarem macierzy, po której następują elementy na STDIN:

./infinity < matrix.txt

matrix.txt:

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

Lub wypróbować dokładność 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;
}
Ton Hospel
źródło
Moim zdaniem jest to imponujące i znacznie lepsze niż naiwne podejście. Czekam na dokumentację :)
1
@TonHospel Czy możesz nam powiedzieć coś więcej o tym, jak wymyśliłeś granice?
flawr
2
Korzystam z Arch Linux i potrzebowałem -lgslcblasflagi do skompilowania. Nawiasem mówiąc, niesamowita odpowiedź!
Rhyzomatic 12.04.16
2

Python 3

12 sekund n = 8 na moim komputerze, rdzeń Ubuntu 4.

Naprawdę naiwny, nie mam pojęcia, co robię.

from itertools import product
from math import e

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

N = 2

n = [1]

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

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

Będzie to zwiększało zakres wykorzystywanych zasobów Z, dopóki nie uzyska wystarczająco dobrej odpowiedzi. Napisałem własne mnożenie macierzy, powinienem używać numpy.

Maltysen
źródło
Dzięki ! Czy możesz pokazać niektóre wyniki i czasy na swoim komputerze?
Twój kod działa pypy, co jest świetne i szybkie. Niestety [[6,0, -1,0, -3,0, 1,0, 3,0, -1,0, -3,0, 1,0, 3,0], [-1,0, 6,0, -1,0, -5,0, 1,0, 5,0, -1,0, -5,0, 1,0 ], [-3,0, -1,0, 6,0, 1,0, -5,0, -1,0, 5,0, 1,0, -5,0], [1,0, -5,0, 1,0, 6,0, -1,0, -5,0, 1,0, 5,0, -1,0] , [3,0, 1,0, -5,0, -1,0, 6,0, 1,0, -5,0, -1,0, 5,0], [-1,0, 5,0, -1,0, -5,0, 1,0, 6,0, -1,0, -5,0, 1,0], [-3,0, -1,0, 5,0, 1,0, -5,0, -1,0, 6,0, 1,0, -5,0], [1,0, -5,0, 1,0, 5,0, -1,0, -5,0, 1,0, 6,0, -1,0], [ 3,0, 1,0, -5,0, -1,0, 5,0, 1,0, -5,0, -1,0, 6,0]] daje tylko złą odpowiedź.
8.1443647932-8.14381938863 = 0,00054540457> 0,0001.
3
@Maltysen Twój program sprawdza tylko, czy ostatni termin jest krótszy niż podana dokładność. Ale błąd, który popełniasz, jest zdecydowanie większy, ponieważ musisz wziąć pod uwagę sumę wszystkich innych terminów dotyczących błędu!
flawr