Я написал 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 = время [мс]):
Мое понимание производительности наихудшего случая здесь заключается в том, что объем работы, выполняемой в главном цикле while, пропорционален N .Основной цикл требует log N (C) (где C> 1 - это постоянное отношение начального диапазона поиска к запрашиваемой точности) для выполнения итераций.Это приводит к следующему уравнению:
График выглядит очень похоже на фиолетовую кривую, которую я использовал выше для приблизительного соответствия данных:
Теперь яУ меня есть (надеюсь, правильные) выводы и вопросы:
- Во-первых, действительно ли мой подход верен вообще?
- Действительно ли описанная мной функция описывает отношение между числом операций и N ?
- Мне кажется, это самый интересный вопрос - мне интересно, что вызывает столь существенную разницу между N = 2 и всеми остальными?Я последовательно получаю этот шаблон для всех моих тестовых данных.
Более того:
- Эта функция имеет минимум в e≈2.718 ... , чтоближе к 3, чем к 2. Также, f (3) .Учитывая, что уравнение, которое я придумал, является правильным, я думаю, это будет означать, что метод трисекции на самом деле может быть более эффективным, чем бисекция.Это кажется немного противоречивым, но измерения, кажется, подтверждают это.Это правильно?
- Если так, то почему метод трисекции кажется настолько непопулярным по сравнению с бисекцией?
Спасибо