Производительность алгоритма поиска корня n-раздела - PullRequest
2 голосов
/ 20 сентября 2019

Я написал n-секционный алгоритм для нахождения корней функции.Принцип работы точно такой же, как и при делении пополам, только диапазон вместо этого делится на N равных частей.Это мой код C:

/*
    Compile with: clang -O3 -o nsect nsect.c -Wall -DCOUNT=5000000 -DNSECT=? -funroll-loops
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>

#define N 6

#ifndef COUNT
#error COUNT not defined!
#endif

#ifndef NSECT
#define NSECT 2
#endif

float polynomial[N];
float horner( const float poly[N], float x )
{
    float val = poly[N-1];
    for ( int i = N - 2; i >= 0; i-- )
        val = poly[i] + x * val;
    return val;
}

float f( float x )
{
    return horner( polynomial, x );
}

float nsect( float a, float b, const float eps_x )
{
    float fa = f( a );
    float fb = f( b );
    if ( fa == 0 ) return a;
    else if ( fb == 0 ) return b;
    else if ( fa * fb > 0 ) return 0;

    float x[NSECT];
    float fx[NSECT];

    while ( b - a > eps_x )
    {   
        x[0] = a;
        if ( ( fx[0] = f( x[0] ) ) == 0 ) return x[0];

        int found = 0;
        for ( int i = 0; i < NSECT - 1; i++ )
        {
            x[i + 1] = a + ( b - a ) * (float)( i + 1 ) / NSECT;
            if ( ( fx[i + 1] = f( x[i + 1] ) ) == 0 )
                return x[i + 1];
            else if ( fx[i] * fx[i + 1] < 0 )
            {
                a = x[i];
                b = x[i + 1];
                found = 1;
                break;
            }
        }

        if ( !found )
            a = x[NSECT - 1]; 

    }

    return ( a + b ) * 0.5f;
}

int main( int argc, char **argv )
{
    struct timeval t0, t1;
    float *polys = malloc( COUNT * sizeof( float ) * N );
    float *p = polys;
    for ( int i = 0; i < COUNT * N; i++ )
        scanf( "%f", p++ ); 

    float xsum = 0; // So the code isn't optimized when we don't print the roots
    p = polys;
    gettimeofday( &t0, NULL );
    for ( int i = 0; i < COUNT; i++, p += N )
    {
        memcpy( polynomial, p, N * sizeof( float ) );
        float x = nsect( -100, 100, 1e-3 );
        xsum += x;

        #ifdef PRINT_ROOTS
            fprintf( stderr, "%f\n", x );
        #endif
    }
    gettimeofday( &t1, NULL );

    fprintf( stderr, "xsum: %f\n", xsum );
    printf( "%f ms\n", ( ( t1.tv_sec - t0.tv_sec ) * 1e6 + ( t1.tv_usec - t0.tv_usec ) ) * 1e-3 );
    free( polys );
}

РЕДАКТИРОВАТЬ: Это команда, которую я использовал для компиляции кода: clang -O3 -o nsect nsect.c -Wall -DCOUNT=5000000 -DNSECT=? -funroll-loops.Я запустил все на i7-8700k.

Я решил проверить производительность алгоритма для разных значений N .Тест состоит из измерения времени, необходимого для нахождения любого корня в диапазоне (-100; 100) для каждого из 5 000 000 полиномов степени 5. Полиномы генерируются случайным образом и имеют действительные коэффициенты в диапазоне от -5 до 5. Значения полинома рассчитываются с использованиемМетод Хорнера.

Это результаты, которые я получил при выполнении кода 10 раз для каждого N (x = N , y = время [мс]): image

Мое понимание производительности наихудшего случая здесь заключается в том, что объем работы, выполняемой в главном цикле while, пропорционален N .Основной цикл требует log N (C) (где C> 1 - это постоянное отношение начального диапазона поиска к запрашиваемой точности) для выполнения итераций.Это приводит к следующему уравнению:

Equation

График выглядит очень похоже на фиолетовую кривую, которую я использовал выше для приблизительного соответствия данных: image

Теперь яУ меня есть (надеюсь, правильные) выводы и вопросы:

  • Во-первых, действительно ли мой подход верен вообще?
  • Действительно ли описанная мной функция описывает отношение между числом операций и N ?
  • Мне кажется, это самый интересный вопрос - мне интересно, что вызывает столь существенную разницу между N = 2 и всеми остальными?Я последовательно получаю этот шаблон для всех моих тестовых данных.

Более того:

  • Эта функция имеет минимум в e≈2.718 ... , чтоближе к 3, чем к 2. Также, f (3) .Учитывая, что уравнение, которое я придумал, является правильным, я думаю, это будет означать, что метод трисекции на самом деле может быть более эффективным, чем бисекция.Это кажется немного противоречивым, но измерения, кажется, подтверждают это.Это правильно?
  • Если так, то почему метод трисекции кажется настолько непопулярным по сравнению с бисекцией?

Спасибо

1 Ответ

3 голосов
/ 23 сентября 2019

Я комментирую общий вопрос, а не ваш конкретный код, который я не до конца понимаю.

Предполагая, что существует единственный корень, который, как известно, находится в интервале длины L и желаемогоТочность ε, вам понадобятся этапы деления log (L / ε) / log (N).Каждый этап подразделения требует N-1 оценки функции (не N), чтобы указать, какой подинтервал среди N содержит корень.

Следовательно, без учета издержек общая стоимость пропорциональна (N-1) / log(Н).Значения этого отношения, начиная с N = 2:

1.44, 1.82, 2.16, 2.49, 2.79...

и выше.

Следовательно, теоретический оптимум при N = 2.Вот почему трисекция не используется.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...