n 段求根算法的性能
Performance of n-section root finding algorithm
我写了一个 n 截面算法来求函数的根。工作原理与二分法完全一样,只是将运行ge分成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 值的算法性能。该测试包括测量在 运行ge (-100;100) 中为 5,000,000 个 5 次多项式中的每一个找到任何根所需的时间。多项式是 运行domly 生成的并且具有实系数 运行从 -5 到 5。多项式值是使用 Horner 方法计算的。
这些是我从 运行 代码中得到的结果,每个 N (x=N, y=时间 [毫秒]):
我对这里最坏情况性能的理解是,主 while 循环中要完成的工作量与 N 成正比。主循环需要 logN(C) (其中 C > 1 是一个常数 - 比率初始搜索 运行ge 到要求的精度)迭代完成。这会产生以下等式:
该图看起来与我在上面用来大致匹配数据的紫色曲线非常相似:
现在,我有一些(希望是正确的)结论和问题:
- 首先,我的方法是否有效?
- 我想出的函数真的描述了运算次数和N之间的关系吗?
- 我认为这是最有趣的一个 - 我想知道是什么导致 N=2 和所有其他人之间存在如此显着的差异?对于我的所有测试数据,我始终得到这种模式。
此外:
- 该函数在 e≈2.718... 中有最小值,它更接近 3 而不是 2。此外,f(3) < f( 2)成立。鉴于我得出的方程式是正确的,我认为这意味着三等分法实际上可能比二等分法更有效。这似乎有点违反直觉,但测量结果似乎承认了这一点。这样对吗?
- 如果是这样,为什么与二分法相比,三分法似乎如此不受欢迎?
谢谢
我评论的是一般问题,而不是你的特定代码,我不完全理解。
假设已知有一个单根位于长度为L的区间内,所需精度为ε,则需要log(L/ε)/log(N)个细分阶段。每个细分阶段都需要对函数进行 N-1 次评估(不是 N),以判断 N 中哪个子区间包含根。
因此,忽略开销,总成本与 (N-1)/log(N) 成正比。这个比率的值是,从 N=2:
1.44, 1.82, 2.16, 2.49, 2.79...
及更高。
因此,理论上的最优值是 N=2。这就是不使用三等分的原因。
我写了一个 n 截面算法来求函数的根。工作原理与二分法完全一样,只是将运行ge分成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 值的算法性能。该测试包括测量在 运行ge (-100;100) 中为 5,000,000 个 5 次多项式中的每一个找到任何根所需的时间。多项式是 运行domly 生成的并且具有实系数 运行从 -5 到 5。多项式值是使用 Horner 方法计算的。
这些是我从 运行 代码中得到的结果,每个 N (x=N, y=时间 [毫秒]):
我对这里最坏情况性能的理解是,主 while 循环中要完成的工作量与 N 成正比。主循环需要 logN(C) (其中 C > 1 是一个常数 - 比率初始搜索 运行ge 到要求的精度)迭代完成。这会产生以下等式:
该图看起来与我在上面用来大致匹配数据的紫色曲线非常相似:
现在,我有一些(希望是正确的)结论和问题:
- 首先,我的方法是否有效?
- 我想出的函数真的描述了运算次数和N之间的关系吗?
- 我认为这是最有趣的一个 - 我想知道是什么导致 N=2 和所有其他人之间存在如此显着的差异?对于我的所有测试数据,我始终得到这种模式。
此外:
- 该函数在 e≈2.718... 中有最小值,它更接近 3 而不是 2。此外,f(3) < f( 2)成立。鉴于我得出的方程式是正确的,我认为这意味着三等分法实际上可能比二等分法更有效。这似乎有点违反直觉,但测量结果似乎承认了这一点。这样对吗?
- 如果是这样,为什么与二分法相比,三分法似乎如此不受欢迎?
谢谢
我评论的是一般问题,而不是你的特定代码,我不完全理解。
假设已知有一个单根位于长度为L的区间内,所需精度为ε,则需要log(L/ε)/log(N)个细分阶段。每个细分阶段都需要对函数进行 N-1 次评估(不是 N),以判断 N 中哪个子区间包含根。
因此,忽略开销,总成本与 (N-1)/log(N) 成正比。这个比率的值是,从 N=2:
1.44, 1.82, 2.16, 2.49, 2.79...
及更高。
因此,理论上的最优值是 N=2。这就是不使用三等分的原因。