魔法书1.1补充:平方根算法的论证

SICP 1.1 节「程序设计的基本元素」给出了一个计算平方根的算法,原书并没有给出证明,大概是因为不太适合书的主题,我在此写这个证明只是为了补充书上所缺,满足自己的强迫症,如有错误请指出!

证明

数学表达如下:

对于任意实数 [tex]b>0[/tex],数列 [tex]a_n[/tex] 的递推公式是 [tex]a_{n+1}=\frac{a_{n}+b/a_{n}}{2}[/tex],且[tex]a_1=1[/tex],那么当 n 充分大时,[tex]a_n[/tex] 可以无限接近 [tex]\sqrt{b}[/tex]。下面给出证明。

由 [tex]a_{n+1}-a_{n}=\frac{b-a_n^2}{2a_n}[/tex] 可知,当 [tex]a_n<\sqrt{b}[/tex] 时,数列 [tex]a_n[/tex] 是严格递增的;当 [tex]a_n>\sqrt{b}[/tex] 时,数列 [tex]a_n[/tex] 是严格递减的。容易知道,当 b > 1 时,数列 [tex]a_n[/tex] 严格递增且有上界;当 b = 1 时,无疑 [tex]a_n=1[/tex] 恒成立;当 b < 1 时,数列 [tex]a_n[/tex] 严格递减且有下界。因此确信数列 [tex]a_n[/tex] 收敛且存在极限,记作 [tex]a[/tex]。

递推式两边取极限得:[tex]a=\frac{a+\frac{b}{a}}{2}[/tex]

整理即得:[tex]a=\sqrt{b}[/tex]

也就是,当 n 充分大时,[tex]a_n[/tex] 可以无限接近 [tex]\sqrt{b}[/tex]。

注记

我一开始认为 [tex]a_1[/tex] 是可以选取任意正数的,后来计算一番发现,存在特定的组合 [tex]a_1[/tex] 和 [tex]b[/tex],使得数列不收敛。如当 [tex]{a_1}^2=\frac{\sqrt{2}+1}{3}b[/tex] 时,[tex]a_3 = a_1[/tex],此时所有奇数项相等,所有偶数项相等,数列围绕 [tex]\sqrt{b}[/tex] 摆动而不收敛。尽管在计算机上存在浮点数舍入误差,不一定能产生预期的效果,但是只有当[tex]a_1=1[/tex]时才能在数学上严格证明出对于所有正数 [tex]b[/tex] 都一定能求出算数平方根的近似值。

实际上迭代算平方根的方法很多,比如说还可以由 [tex]c^2-1=b-1[/tex] 得到 [tex]c_{n+1}=\frac{b-1}{c_n+1}+1[/tex],可以证明该数列的极限也是 [tex]\sqrt{b}[/tex](选取恰当的 [tex]c_1[/tex])。在此不再赘述。

在 Tcl 里面创建一个 Gtk+ 窗口…

This may sound ridiculous, but that's actually happening...

package require critcl
foreach item [exec pkg-config --cflags gtk+-3.0] {
	critcl::cheaders $item
}
foreach item [exec pkg-config --libs gtk+-3.0] {
	critcl::clibraries $item
}
critcl::ccode {
#include <gtk/gtk.h>
}
critcl::cproc gtk {} void {
	gtk_init(NULL, NULL);
	GtkWidget *win = gtk_window_new(GTK_WINDOW_TOPLEVEL);
	gtk_widget_show(win);
	gtk_main();
}
gtk

"apt-get install critcl libgtk3.0-dev", and be sure to kill it using C-c.

Critcl is amazing, isn't it?

递推式的迭代表达

考虑 Fibonacci 数列:f(0)=0, f(1)=1, f(n)=f(n-1)+f(n-2) (n=2, 3, 4, ...)

可以很轻松地写出其递归形式表达(这里使用 Racket):

(define (fib n)
  (if (< n 2)
      n
      (+ (fib (- n 1)) (fib (- n 2)))))

很显然,这个函数是无法被尾递归优化的。但是我们也可以很简单地将它转化为线性迭代过程:

(define (fib-iter a b n)
  (if (= n 0)
      a
      (fib-iter b (+ a b) (- n 1))))
(define (fib n)
  (fib-iter 0 1 n))

这样做的原理是什么呢?实际上,我们只不过是把每次运算的结果丢弃掉,然后换上下一次需要的数据而已。比如说,我们要计算 f(n+1),那么就需要知道 f(n-1) 和 f(n)。循环的终止条件是 n=0,然后返回 a。如何验证这样做是正确的?可以列出如下表:

i=n     a=f(0)   b=f(1)
i=n-1   a=f(1)   b=f(2)
i=n-2   a=f(2)   b=f(3)
...
i=1     a=f(n-1) b=f(n)
i=0     a=f(n)   b=f(n+1)

我们得到了我们想要的结果 f(n)!循环终止条件不用 n=1 的理由是:如果一开始传入的 n 就是 0,那么循环将永远无法终止。而使用 n=1 所得到的效率优化微乎其微、可以忽略不计。

可以用完全相同的思路解决 SICP 练习 1.11:写出计算 f(n)=n (n<3), f(n)=f(n-1)+2f(n-2)+3f(n-3)(n>=3) 的递归和迭代的代码。

递归非常简单且自然:

(define (foo n)
  (cond ((< n 3) n)
        (else (+ (f (- n 1)) (* 2 (f (- n 2))) (* 3 (f (- n 3)))))))

使用相同的表格:

i=n    a=f(0)   b=f(1)   c=f(2)
i=n-1  a=f(1)   b=f(2)   c=f(3)
i=n-2  a=f(2)   b=f(3)   c=f(4)
...
i=0    a=f(n)   b=f(n+1) c=f(n+2)

可以帮助我们写出正确的迭代代码:

(define (foo-iter a b c count)
  (if (= count 0)
      c
      (foo-iter b c (+ (* 3 a) (* 2 b) c) (- count 1))))
(define (foo n)
  (foo-iter 0 1 2 n))

可以使用循环不变式证明该代码是正确的。

将以上结论推广,可以得到:

定义在 N 上的函数 f(n),对于 n>=N,有 f(n)=A1f(n-1)+A2f(n-2)+A3f(n-3)+...,其中 A1, A2, ... 都是常数,那么对于任意 n<N,f(n) 都是需要另行定义的。其递归伪代码必为如下格式:

(define (f n)
  (cond ((n=1) ...)
        ((n=1) ...)
        (...)
        ((n=N) ...)
        (else (+ (* A1 (f (- n 1))
                 (* A2 (f (- n 2))
                 ...))))

其迭代代码一定需要 N 个变量,必为如下格式:

(define (f-iter a b c d e ... count)
  (if (= count 0)
      a
      (f-iter b c d e ... (+ ...) (- count 1))))
(define (f n)
  (f-iter f0 f1 f2 ... fN n))

当 N 比较大时,可能 f(n+1), f(n+2), f(n+3), ..., f(n+N-1) 的计算成本就无法忽略了,这时我们可以在 f 里面想想办法。

yield

今天看 C# 发现这样一种用法:

class IntManager
{
	int[] arr = { 0, 1, 2, 3 };
	public IEnumerator<int> GetEnumerator()
	{
		for (int i = 0; i<4; i++)
		{
			yield return arr[i];
		}
	}
}

是的……对我这样孤陋寡闻的人来说,这就是这么神奇……居然可以这样写出迭代器!

而且根据一些资料来看,这个迭代器是随用随取的,也就是说如果资源不在内存上,那就可以需要多少取多少,不会一次性全部取回来,非常方便。

以前我只知道 yield 用在 coroutine 上。不过把这两个一结合,我就发现了神奇的现象……Tcl 的 yield 其实也不过是这样而已。

proc y {} { 
    set i 0
    while 1 {
        yield [incr i]
    }
}
coroutine coro y ;; -> 1,这里会调用一次 y 所以是 1
coro ;; -> 2
coro ;; -> 3
coro ;; -> 4

没错,这还是一个……嗯,我们不如来研究一些好玩的。

proc fibo {} {
    ;# 第一次 yield 之前可以尝试做一些初始化工作
    yield ;# 然后用 yield 返回结果,这里自然是不需要的
    set a 0
    set b 1
    while 1 {
        ;# a,b = b,a+b
        set t $b
        set b [expr $a+$t]
        set a $t
        yield $a
    }
}
coroutine next fibo
next ;; -> 1
next ;; -> 1
next ;; -> 2
-> 就是 1,1,2,3,5,8,13,21,34,55,...

实际上用这个求这个数列的效率不高,因为求过一次数据就丢弃了,不如以前写过的一个缓存版本快。但是这个版本不会爆栈(相比递归版本),因为根本只有一个栈帧和循环版的效率理论上来说差不多,没测试。

CSAPP 2.1:字节序、大端与小端

终于比较认真地读了一下 CSAPP 的一点点,顺便做点记录,也算是以往经验的小小总结吧。诸如 [DCL13-C] 的标注,是 CERT C 安全编码标准的编号(看来我写的代码还是不错的嘛)。

大端、小端是机器中两种表示数据的方法(字节序),大端法是把最高有效数字排在低地址,小端法是把最高有效数字排在高地址。举例来说,0x12345678 这个数字,在小端机器中会被存储为 78 56 34 12(低->高),相反在大端机器中是 12 34 56 78(低->高)。

以程序验证:

#include <stdio.h>
#include <stdlib.h>

void show_bytes(const void *bytes, size_t len)
{
	const unsigned char *bytes_; // [DCL00-C], [INT07-C]
	size_t i;

	bytes_ = bytes;
	for (i = 0; i < len; i++)
	{
		printf(" %.2X", bytes_[i]);
	}
	printf("\n");
}

int main(void) {
	int val = 0x12345678;
	show_bytes(&val, sizeof(val));
	return EXIT_SUCCESS;
}

运行该程序,输出 78 56 34 12,说明我的机器就是小端的。

注意,在这里我稍微修改了一下 CSAPP 上的例程。有如下修改:

  1. 将 byte_pointer 修改为了 const void *:任何指针都能被隐式转换为 void * 而没有警告,方便使用。另外加上 const 确保不会误修改传入数据。[DCL13-C]
  2. 将 len 的类型修改为了 size_t:对 len 来说有符号是毫无意义的,len 不可能小于 0。实际上,在运行中如果传入负数将会溢出(指针都无符号)。[INT01-C]
  3. 这里的计数器 i 并不是通常的 int 类型,而是与和它比较的 len 一致,关于这一点可以看 CS:APP (中文版)的第 53 页关于 FreeBSD getpeername 漏洞的讨论。[INT02-C]

那么如何判断大小端呢?在运行时可以用类似上面的方法这样判断:

bool bigendian()
{
	static int val = 0x12345678;
	if (*(((char *)&val)+3) == 0x78) // warning: magic number 3
	{
		return true;
	}
	return false;
}

或者用 union:

union
{
	int val;
	char ch[4];
} _bigendianstub;

bool bigendian() // requires C99-compliant compilers; stdbool.h must be included
{
	_bigendianstub.val = 0x12345678;
	if (_bigendianstub.ch[3] == 0x78)
	{
		return true;
	}
	return false;
}

一般来说,在一个机器中只能是大端或小端,而操作系统也必然与之对应,于是可以直接在编译时确定:

const char ch[4] = { 0x00, 0xff, 0xff, 0xff }; // [DCL00-C]
#define LITTLEENDIAN ((*(int*)ch)<0) // [PRE02-C]
#define BIGENDIAN (!LITTLEENDIAN) // [PRE02-C]

这里比较流氓,直接把数字写到数据段里面去了,然后利用了 int 最高位表示符号的特点解决问题。相当于直接把 0x00FFFFFF 当作 int,显然我们发现,在大端机器中这个数字表示 0x00FFFFFF,在小端机器中这个数字表示 -256。

当然,以上不过是示例而已,在真实环境中,我们有编译器提供的宏来确定。gcc 的相关宏是 __BYTE_ORDER。

说了这么多,大小端究竟有什么影响呢?主要影响就是与其他机器交互了。个人电脑基本都是小端,而网络字节序则是大端,意味着我们要在发送数据时进行转换。操作系统提供了 ntoh(network to host)、hton (host to network)系列函数可以实现这个目的。但是在大端机器上,这两个函数显然无法得到小端结果!因此我们自己实现一套,以下是完整的测试代码:

#include <stdio.h>
#include <stdlib.h>

const char ch[4] = { 0x00, 0xff, 0xff, 0xff };
#define LITTLEENDIAN ((*(int*)ch)<0)
#define BIGENDIAN (!LITTLEENDIAN)

void show_bytes(const void *bytes, size_t len)
{
	const unsigned char *bytes_;
	size_t i;

	bytes_ = bytes;
	for (i = 0; i < len; i++)
	{
		printf(" %.2X", bytes_[i]);
	}
	printf("\n");
}

int swap_byteorder(int orig)
{
	int ret;
	char *ptrf, *ptrl;

	ret = orig;
	ptrf = (char*)&ret;
	ptrl = ((char*)&ret)+3; // warning: magic number 3
	for (; ptrf < ptrl; ptrf++, ptrl--)
	{
		int t;
		t = *ptrf;
		*ptrf = *ptrl;
		*ptrl = t;
	}
	return ret;
}

inline int tobe(int orig) // requires C99-compliant compilers
{
	if(BIGENDIAN) return orig;
	return swap_byteorder(orig);
}

inline int tole(int orig) // requires C99-compliant compilers
{
	if(LITTLEENDIAN) return orig;
	return swap_byteorder(orig);
}

int main(void) 
{
	int orig, le, be;

	orig = 0x12345678;
	le = tole(orig);
	be = tobe(orig);
	printf("original (%s):", LITTLEENDIAN?"little endian":"big endian");
	show_bytes(&orig, sizeof(int));
	printf("little endian:");
	show_bytes(&le, sizeof(int));
	printf("big endian:");
	show_bytes(&be, sizeof(int));
	return EXIT_SUCCESS;
}

只实现了 int 的大小端转换,实际上其余类型乃至任意类型的大小端转换都很容易实现,故从略。

提醒:字符串(特指 char[])没有大小端转换,每个 char 就是一个字节,从低地址端开始排列(想想「字节序」这个名称)。实际上我们上面的一切都是依赖这个的。另外还依赖了「64 位、32 位系统 int 都是 4 字节」。如果需要很好的可移植性,这些还不够。以上关于字节序的讨论有一个直接推论:跨机器交换文件应当尽量使用文本文件,二进制文件必须设计良好、经过大量测试才能使用。

(翻到后面才发现好多内容都是作业……)

C++11 的 decltype 一种用法

#include <string>
#include <iostream>
using namespace std;

class clsA
{
public:
    int value()
    {
        return 1;
    }
};

class clsB
{
public:
    string value()
    {
        return "yes/no";
    }
};

template<class Type>
auto foo(Type x) -> decltype(x.value())
{
    return x.value();
}

int main(void)
{
    clsA oa;
    clsB ob;
    cout << foo(oa) << foo(ob) << endl;
    return 0;
}

(lib)Tcl 现代方法:变量

之前有人问到如何从 C 端访问 Tcl 变量,此处做一总结。庶几不失「现代」之名。有关 Tcl_Obj 的维护者资料击此

于是终于回到了主线。《Tcl 现代方法》本来是作为记录 Tcl 现代用法的系列博文,一开始是介绍 C 库的,但是后来 Tcl 脚本的内容反倒更多,结果现在纯粹谈起来 C 库又会导致歧义,因此在标题前加上 lib 标记,以示内容是 C 库而非 Tcl 脚本。

继续阅读

Tcl 现代方法:你知道几种写类的方法?

本文是从底向上逐步构建其 TclOO 大框架的,如果你更喜欢自顶向上的,可以参考别的教程。

另外,本文只是入门级别,读完本文后再去看文档或者其他教程应该会容易些,毕竟是汉语写的。这里我推荐 Ashok 写的《Tcl Programming for Windows》中的 TclOO 一章。

除去协程,面向对象编程也是现代编程语言不可或缺的一部分。Tcl 在 8.6 之前都没有内建支持面向对象编程,但这不意味着 Tcl 没法 OO——恰恰相反,Tcl 有一大堆 OO 库,我随口就能叫出几个出名的:[incr Tcl]、XOTcl、Snit。这些库各有千秋,有的模仿 C++,有的是原型继承,还有的是模仿 Common Lisp Object System,千姿百态。可能和大家所想不同,这样的「大好」局势并未给 Tcl 带来许多好处:依然有人批评 Tcl 没有 OO 支持;有的人虽然知道有,但批评选择太多让用户无所适从,或者为自己所爱的库口诛笔伐;有的人虽然不在乎别人怎么选择,但很苦恼其他扩展包依赖了自己根本不需要的 OO 库……在这种情况下,Tcl Core Team 站了出来,提出了 TIP 257,由 Donal K. Fellows 主要实现。这套面向对象机制起名叫做 TclOO,并随 Tcl 本身提供。TclOO 的设计十分灵活,可以以此为基础重新制作 [incr Tcl] 等库,让大家其乐融融地生活在一起。

思考一下,开源系统的桌面环境是不是也是这样呢?

继续阅读

Tcl 现代方法:协程(coroutine)

协程已经是现代语言的必备要素了。Tcl 的协程仿效的是 Lua,在 8.6 版本中引入。(感谢 NRE 引擎!)

 

继续阅读

Tcl 现代方法:TclHttpd 入门

TclHttpd 是 Tcl 界人尽皆知的有名软件,功能强大,易于使用。在其原作者不再维护之后,Cliff 和 Sean 于最近接手维护,使之支持最新的 8.6。

继续阅读