编程求最大最小值的函数（典型算法思想与应用2）

迭代法（Iteration）也称辗转法，是一种不断用变量的旧值递推新值的过程，跟迭代法相对应的是直接法（或者称为一次解法），即一次性解决问题。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令（或一定步骤）重复执行，在每次执行这组指令（或这些步骤）时，都从变量的原值推出它的一个新值。

// 自然数连加和 int adds1(int n) { int sum = 0; for(int i=1;i<=n; i ) sum =i; // 迭代法 return sum; } int adds2(int n) { return (1 n)*n/2; // 直接法 }

利用迭代算法解决问题，需要做好以下三个方面的工作：

① 确定迭代变量

在可以用迭代算法解决的问题中，至少存在一个可直接或间接地不断由旧值递推出新值的变量，这个变量就是迭代变量。

② 建立迭代关系式

所谓迭代关系式，指如何从变量的前一个值推出其下一个值的公式（或关系）。迭代关系式的建立是解决迭代问题的关键。

递推法也是一种利用旧值推新出新值的方法，但与迭代法不同的时，递推法往往会有一个初始条件，就如同分段函数一样，初始条件是递推表达式之一，而迭代法一般使用一个单一表达式。

③ 对迭代过程进行控制

在什么时候结束迭代过程？这是编写迭代程序必须考虑的问题。不能让迭代过程无休止地执行下去。迭代过程的控制通常可分为两种情况：一种是所需的迭代次数是个确定的值，可以计算出来；另一种是所需的迭代次数无法确定。对于前一种情况，可以构建一个固定次数的循环来实现对迭代过程的控制；对于后一种情况，需要进一步分析得出可用来结束迭代过程的条件。

迭代法又分为精确迭代和近似迭代。前者如辗转相除法，后者如“二分法”、"牛顿迭代法”来近似求函数值的问题。

1 二分法求函数值

对于函数y=f(x)，设x的取值区间为[a,b]满足f(a)*f(b)<0k，也就是是f(x)在x的取值区间有0值，通过不断地把函数f(x)的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫二分法。

用二分法求函数f(x)零点近似值的步骤如下：

1.1 确定区间[a，b]，验证f(a)·f(b)<0，给定精确度ξ。

1.2 求区间(a，b)的中点c。

1.3 计算f(c)。

(1) 若f(c)=0，则c就是函数的零点；

(2) 若f(a)·f(c)<0，则令b=c；

(3) 若f(c)·f(b)<0，则令a=c。

(4) 判断是否达到精确度ξ：即若|a-b|<ξ，则得到零点近似值a(或b)，否则重复2-4。

while (fabs(b-a) > e) { double c = (a b)/2.0; if (f(a)* f(c) < 0) b = c; else a = c; }

测试代码：

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <assert.h> double f(double x) // 测试函数f(x)=1 x-x*x*x { return 1 x-x*x*x; } double bisect(double a,double b) { double e = 1e-5; if (fabs(f(a)) <= e) printf("solution: %lg\n", a); else if (fabs(f(b)) <= e) printf("solution: %lg\n", b); else if (f(a)*f(b) > 0) printf("f(%lg)*f(%lg) > 0 ! need <= 0 !\n", a, b); else { while (fabs(b-a) > e) { double c = (a b)/2.0; if (f(a)* f(c) < 0) b = c; else a = c; } return (a b)/2.0; } } int main() { printf("solution: %lg\n", bisect(1,2)); // solution: 1.32472 getchar(); return 0; }

二分法求平方根：

double sqrt1(double x) // 二分法 { double EPSINON = 1e-5; double low = 0.0; double high = x; double mid = (low high) / 2; while ((high - low) > EPSINON){ if (mid*mid > x) high = mid; else low = mid; mid = (high low) / 2; } return mid; }

2 牛顿迭代法求平方根

令f(x)=x^2-n，如图所示：

取x0，如果x0不是解，做一个经过(x0,f(x0))这个点的切线，与x轴的交点为x1。

同理，如果x1不是解，做一个经过(x1,f(x1))这个点的切线，与x轴的交点为x2。

依次类推。

以这样的方式得到的 xi 会无限趋近于 f(x)=0 的解。

(f(x)-f(xi))/(x-xi)=f’(x)，f’(x)是斜率也是f(x)的导函数，即f’(x)=2x。

化简得：f(xi)=f(x)-f’(x)(x-xi)，令f(xi)=0得：

(x²-n)-2x(x-xi)=0

持续化简得：

x² - n - 2x² 2xxi=0

2xxi=x² n

2xi=x n/x

xi=(x n/x)/2 // 迭代表达式

double sqrt2(double x) // 牛顿迭代法 { if (x == 0) { return 0; } double last = 0.0; double res = 1.0; while (res != last) // 判断前后两个解 xi 和 xi-1 是否无限接近 { last = res; res = (res x / res) / 2; // 迭代表达式 } return res; }

3 欧几里德辗转相除法

最经典的迭代算法是欧几里德辗转相除法，用于计算两个整数a,b的最大公约数。枚举法虽然也可以求出最大公约数，但效率明显较差。

辗转相除法的计算原理：

gcd(a,b) = gcd(b,a % b)

a可以表示成a = kb r，则r = a % b。假设d是a,b的一个公约数，则有 a%d==0，b%d==0，而r = a - kb，因此r%d==0（因为a%b==0，kb%d==0） ，因此d是(b,a % b)的公约数。

同理，假设d 是(b,a % b)的公约数，则 b%d==0，r%d==0 ，但是a = kb r ，因此d也是(a,b)的公约数。

因此(a,b)和(b,a % b)的公约数是一样的，其最大公约数也必然相等。

欧几里德算法就是根据这个原理来做的，欧几里德算法又叫辗转相除法，它是一个反复迭代执行，直到余数等于0停止：

#include <iostream> using namespace std; int gcdIter(int a,int b) { if (a<=0 || b<=0) return 0; int temp; while (b > 0) { temp = a % b; // 迭代关系式 a = b; b = temp; } return a; } int gcd(int a,int b) { return b?gcd(b,a%b):a; } // 一次gcd()调用有a=b,b=a%b int main() { cout<<gcdIter(48,128)<<endl;//16 cout<<gcd(48,128)<<endl;//16 cin.get(); return 0; }

从上面的程序我们可以看到a，b是迭代变量，迭代关系是temp = a % b；根据迭代关系我们可以由旧值推出新值，然后循环执a = b; b = temp；直到迭代过程结束（余数为0）。

需要注意的是，一些问题可以用迭代法解决，也可以用递归法解决，如上面的辗转相除法。其迭代表达式隐含在函数递归调用的参数表达式中，通过隐式栈存储，其返回值也是如此。

附迭代法求平方根代码：

#include <iostream> #include <cmath> using namespace std; double sqrt1(double x) // 二分法 { double EPSINON = 1e-5; double low = 0.0; double high = x; double mid = (low high) / 2; while ((high - low) > EPSINON){ if (mid*mid > x) high = mid; else low = mid; mid = (high low) / 2; } return mid; } double sqrt2(double x) // 牛顿迭代法 { if (x == 0) { return 0; } double last = 0.0; double res = 1.0; while (res != last) // 判断前后两个解 xi 和 xi-1 是否无限接近 { last = res; res = (res x / res) / 2; // 迭代表达式 } return res; } double sqrt3( double x ) { float y; float delta; float maxError; if ( x <= 0 ) return 0; y = x / 2; // initial guess maxError = x * 0.00001; // refine do { delta = ( y * y ) - x; y -= delta / ( 2 * y ); } while ( delta > maxError || delta < -maxError ); return y; } double sqrt4(double x) // 引自John Carmack 的Quake引擎中的源代码 { if(x == 0) return 0; float result = x; float xhalf = 0.5f*result; // float后加f转换成double类型 int i = *(int*)&result; i = 0x5f3759df - (i>>1); // 通过预先猜测的一个神奇数字 0x5f3759df, 来将迭代次数降到极限的一次 // (另外通过整形数的移位来进一步加速) result = *(float*)&i; result = result*(1.5f-xhalf*result*result); result = result*(1.5f-xhalf*result*result); return 1.0f/result; } int main() { cout<<sqrt(2)<<endl; //1.41421 cout<<sqrt1(2)<<endl; //1.41421 cout<<sqrt2(2)<<endl; //1.41421 cout<<sqrt3(2)<<endl; //1.41421 cout<<sqrt4(2)<<endl; //1.41421 cin.get(); return 0; }

－End－

,