引言

很少有程序员需要自己写开平方函数，因为基本所有的语言都自带开平方的函数。更不用说1999年发布的SSE指令集，从CPU指令的级别上支持了开平方运算。当然这依然阻挡不住程序员们探索各种开平方方法的热情。其中最著名的当属John Carmack在开发Quake III时写的那段用精度换速度的Magic Code。如果你更关心如何提高开方函数的速度，你可以参考这篇比较各种方法速度文章，其中有些比单纯地使用CPU指令sqrtss开平方更快。

本文的关注点并不是速度，你甚至可以在本文中看到或许是业界最慢的一些方法。本文专注从软件工程领域的各各视角考量几种风格迥异的开平方函数，通过它们展示如何让代码的可重用度更高、使用上的灵活性更强、复杂度更低，以便于理解及维护。你总是可以在需要的时候，比较容易地提高一个具有高可维护性、设计良好的系统的性能。而一个系统一但由于代码质量上的问题，进入了一种难于理解、难于维护、难以变动、不可重用、无法预测运行结果的状态，如果再想提高其性能什么的，基本上就只能重写了。

导读

本文代码使用的语言是F#。但是你不必担心自己是否学过这个语言，当英语起来应该不会有什么问题。 本文中的示例与求解思路主要来自《Structure and Interpretation of Computer Programs 2nd Edition》，书中使用的是Scheme语言；类似地，那本书是某大学的计算机专门的入门教材，甚至不要求读者有编程经验。但是会需要微积分的基本知识，否则你应该只能搞懂第一个函数。

如果你很熟悉F#，那太好了，本文中的代码都是我一边看F#编程参考一边写的。如果有什么更好的实现方式，也希望你能指出来。

如果你根本不是计算机专业的，本文也值得一读，它或许能解释为什么你的那些计算机专业的同学们既不会修电脑修手机，也无法帮你找回被你遗忘的某网站的密码。因为下面这些才是他们所能解决的问题。

问题引出

求平方根运算的方法有很多，我们初中的时候就学习过手算开平方，还发过一个小册子，上面有各种数据表，用来直接查询一些数的平方根等。但是那些方面都不能方便才用于计算机：手算开平方规则较多，数据结构也复杂，代码量会很多；查表法需要过多的数据存储空间，也不切实际。

非计算机专业的人，很容易错估计算机领域问题的复杂度，一方面是因为人们总是喜欢这样做：先把问题在自己熟悉的领域重新定义，然后用自己在这个领域所熟悉的知识解决这个新问题，并声称原来的那个问题应该也具有同等的复杂度，并可以类似地解决掉。逻辑不好的人会很吃这套。

事实上，很多计算机专业人士在这方面也不惶多让。通常而言，在计算机领域，针对一个问题，找到一个可行的方案通常都是很简单的，并不需要什么高深的知识和精巧的构思。比如，如果你并没有高等数学的知识，仅仅用初、高中数学知识，你或许可以自己想出折半查找法求x的根，在[0,x)区间内以二分法搜索解，并得出像下面这样的代码：

open LanguagePrimitives

let inline abs x = if x > GenericZero then x else -x
let sqrt_v0 x = 
  let rec search f a b = 
    let close_enough target tolerance value =
      (abs (target - value)) < tolerance
    let avg = (a + b) / 2.0
    if a - b |> abs |> close_enough 0.0 0.001 then
      avg
    else
      match f avg with
      | v when v > 0.0 -> search f a avg
      | v when v < 0.0 -> search f avg b
      | v -> avg
  let half_interval_method f a b =
    let valueA = f a
    let valueB = f b
    match (valueA, valueB) with
    | (x, y) when x < 0.0 && y > 0.0 -> search f a b
    | (x, y) when x > 0.0 && y < 0.0 -> search f b a
    | (x, y) ->
      invalidArg "a" "f a must be of different sign with f b."
  half_interval_method (fun v -> v * v - x) 0.0 x

代码1

这个代码可行，结果也正确，只是有些小问题，比如：代码量多、速度慢、精度低。通常人们不把代码量当问题，因为它不能当作任何标准的度量；如果同时又对性能也没有什么过高要求的话，那么这段代码堪称完美。只是如果哪天需要改进其精度，会发现性能会随精度要求的提高而急剧下降。如果又想同时改进其性能，基本就只能整个重写一遍。这种半调子代码，几乎是所有软件项目的潘多拉魔盒，一但被打开，就会变得无法收拾。

问题的根源

写出高质量的代码，是每个软件工程师的职责，一个软件工程师的水平也可以通过其代码质量而得到体现。所以其实没有人会刻意地去写出像上面那样的半调子代码。出现这样的代码的原因，更多是非主观因素。比如相关知识储备不足、比如项目本身风格约束、比如工程进度的压迫等等。他们专注于问题本身，并尽自己所能写出了自己最满意的代码，却可能意识不到它其实是有多糟糕。这也是软件工程和土木工程的最显著不同：土木危楼，一目了然；软件危机，至死方现。

除了对计算机本身和软件工程本身的研究，多数软件工程项目都不仅仅需要软件方面的知识，更需要把软件工程知识与相关领域知识结合起来才能发挥其价值，这也是软件开发工作相比其它工作最特殊的地方之一。所以人们常说搞软件的要一直学习，这是没错的，因为在软件业，任何已经解决过的问题都不是问题，一方面是要学习如何解决这些问题，另一方面，如果想创造出些新的价值出来的话，势必要做一些别人没有想过、没有做过的事情。这些是题外话。

问题分析

求平方根函数，就是这种跨领域问题，解决起来总有两方面的基本工作要做：

1. 熟悉领域知识。
2. 找到这个问题在计算机领域的合适的解决方式。

开平方是一个数学问题，在着手设计并编写代码之前，从数学分析的角度深入而全面地了解一下这个问题是必不可少的。在这方面的工作都没有做充分的情况下，就会产生出类似上面的代码。

问题的解决大体如是：对问题本身有一个大局观，才能从中找到更合适的解决方案，同时也可以避免陷入自己的先入为主的知识体系的牛角尖而无法自拔。在这个阶段，有这样几个核心问题要回答：

1. 问题所在领域的知识有哪些？
2. 哪些知识可以更好地解决这个问题？
3. 问题的环境（上下文）是怎么样的（问题没有普适的解法，所以需要针对具体环境）？

这些问题或许有些抽象而空洞，我们下面就开始一步步地求解开平方问题。最后这个函数写成什么样都不重要，重要的是我们如何走到那一步的，以及为什么要走那一步。

问题的求解

基本需求

回到问题上来。我们的目标是希望有一个函数，可以用来计算一个数的平方根，并需要满足以下的要求：

• 时间复杂度在特定精度下为常量并稳定。
• 可以灵活地调整所需要的精度。
• 其他程序员能看得懂。

折半查找法并不符合前两个基本要求。所以它不是一个好的方法。

古巴比伦算法

我们只要简单地搜索一下就可以知道常见的开平方解法有哪些。而且不难发现有些方法在计算机系统上实现起来比折半法还要简单。比如公元前2000年，美索不达米亚人所使用的古巴比伦算法。而且不需要什么高深的数学知识，只是我们初高中没有学习过，不广为人知而已。

其计算过程很简单，如果我们猜g是x的一个根，那么$$\frac{g+\frac{x}{g}}{2}$$会是一个更好的根。用代码实现起来也很简单：

let square x = x * x
let sqrt_v1 x = 
  let rec sqrt_it guess =
    if (abs (square guess - x) / x) < 0.001 then
      guess
    else
      sqrt_it ((guess + x / guess) / 2.0)
  sqrt_it 1.0

代码2

代码2中我们把1.0当作所有数的平方根的初始猜测，然后迭代上面的算法直到误差小于0.001。这个函数基本满足前两个需求。只是没有背景知识的话，有些不太好看懂。对于一些对代码质量没有要求的环境，这个函数也已经完美了。

但是我们希望代码好读一些。而有心的人也会问，巴比伦人的算法的原理是什么。这些问题都还没有解决。

函数抽象

一个非常简单的提高代码可读性的手段是给你读不懂的部分起一个名字。

open LanguagePrimitives

let average x y = (x + y) / (GenericOne + GenericOne)
let close_enough target tolerance root = 
    (abs (target - root) / target) < tolerance 
let avg_damping f = fun x -> average x (f x) 
let sqrt_v2 x = 
  let getNext = avg_damping (fun v -> x / v)
  let rec sqrt_it guess =
    if square guess |> close_enough x 0.001 then
      guess
    else
      sqrt_it (getNext guess)
  sqrt_it 1.0

代码3

与前一个版本相比，这个函数其像仅仅是提取出了几个独立的小函数，给他们起了更明白的名字而已。有没有独立的函数，有没有名字都不是重点。每个语句都可以被提取成一个单独的函数，但是那样做一点意义都没有。重点是：哪些需要被提取出来，哪些不需要。引导你做出那些决定的策略才是重点。很多人不屑于做这种小函数的提取，尤其是在这些函数只有自己用的时候。但是有些有价值的东西，只有你这样做了之后才会发现。

比如上面的代码，经过整理之后，你就会发整个代码段中，与开平方运行有关的只有加粗的那个函数。这是一个非常强的信号，这个sqrt_v2函数的求解过程应该有着某种通用性。这个感觉没有错，那通用的部分就是求解方程的不动点。合理的函数提取就是过程抽象的手段，而提取之后又可以发现新的模式可以进一步抽象。过程抽象本身可以推进过程抽象。

于是那个通用的部分就又可以被提取出来了。于是开平方的代码本身就只剩一行了。它的含义就是：函数$$y=\frac{y+\frac{x}{y}}{2}$$的不动点，就是x的根。

let rec find_fix_point f guess = 
  let next = f guess
  if close_enough next guess 0.001 then
    next
  else
    find_fix_point f next

let sqrt_v3 x = 
  find_fix_point (fun v -> average v  (x / v)) 1.0

代码4

比较代码2、3、4，这三个版本的代码都源自古巴比伦算法。只是代码4具有更高的：

• 可重用性：find_fix_point, avg_damping, close_enough都是通用方法。
• 稳定性：基本可以在3次迭代内得到的一个很近似的解。
• 可读性：代码量少，含义明确。

剩下的两个问题

根据我们之前提出的需求，代码4还存在两个主要的问题。

• 可维护性比代码1低：因为依赖更多的领域知识，如果不知道古巴比伦算法，基本看不懂这个代码是怎么写出来的。
• 其精度不可灵活控制：其内部使用了0.001来判断一个解是不是够好。但是显然对于计算$$\sqrt{0.0001}$$就不能正常工作了。

这两个问题是在我们选择古巴比伦法之前没有遇见到的。在软件开发的过程时常会出现类似这种状况：有些问题只有做出来甚至上线之后才会出现。

这两个问题不是数学问题，而是纯粹编程实现上的问题。这也正是我们在分析领域问题时所需要注意的，如何把领域问题与计算机系统结合起来，并找到可能存在的问题及适合的解决方案的地方。

类似的，常见的计算机系统在数学领域的约束还有：

1. 数值表示的精度。
2. 数值表示的范围。
3. 数值运算的截断误差。
4. 不同数值运算间的速度差异。
5. 数值求值策略对数学公式的敏感性。（等价的数学公示可能会有不同的运算结果。有些是数学本身的问题，有些有计算机系统约束的问题。）

我们下面继续分别解决这两个问题。但是不要忘了基于上面的话，它告诉我们：我们为解决这个两问题而引入的技术，很可能会带来新的问题。这是不可控的。只有做了才能发现。这是软件工程复杂性的另一个体现。《Principles of Computer System Design》一书对此进行了更详尽透彻的讨论。在此就不展开了。

牛顿法

当把一个函数化简到代码4的样子，我们不难又可以发现：如果我们不是要开平方，是要开立方。唯一需要改的地方很可能就是那个不动点方程。但是古巴比伦人没有给出开立方的迭代求解方法。即使有立方的，那任意次方的呢？更进一步地，有没有一个不动点方程可以求解任意一元多项式方程呢？牛顿说有：

假设a是对方程 f(x)=0的一个估算解，且f(x)可导，令：

$$b=a-\frac{f(a)}{f'(a)}$$

那么b是个比a更好的近似解，其中f'(a)是 f(x) 的导函数 f'(x)在a点的值。

很少有编程语言能在语言层面支持求导运算。但是任何语言都可以根据导数的如下定义：

$$Dg(x)=\frac{g(x + dx) – g(x)}{dx}$$

实现出一个求取函数f导数值的函数：

let deriv f dx = fun x -> (f (x + dx) - f x) / dx

并得到一个使用牛顿法的开平方函数：

let newton_method f = 
  let Df = deriv f 0.001
  find_fix_point (fun x -> x - f x / Df x) 1.0

let sqrt_v4 x =
  newton_method (fun y -> y * y - x)

代码5

我们看到这个版本书使用的函数$$y=y * y – x$$就是解为x的方程本身。那么即使你不明白牛顿法的工作方式是什么，你大概也可以猜到求立方根的方法可以写成这样：

let cube_root x =
  newton_method (fun y -> y * y * y - x)

代码6

看上去，基本所有的可导函数的解都可以用牛顿法解。但是事情其实不会这么简单。但是对牛顿法的深入讨论不在本文的讨论范围内。如果有兴趣可以自行参看微积分相关书籍。

另外需要注意的是，在代码5中的Df函数中，又引入了一个新的精度假设值0.001。我们提高了代码的可读性和可扩展性，但是进一步损害了其精度（4次迭代之后在小数点第9位会体现出来）。之后会介绍可以解决这个问题的一个方案。

无限延迟求值流

对于精度问题，有一个看似十分简单且可行但是实际没有太大价值的解决方案：把精度当参数之一。它只是把同样的问题抛给了调用者而已。而对调用者而言，我相信他们更希望拿到一个解之后才能真正确定这个解够不够好，他们也无法预知对所有定义域都合理的精度是什么。当然，他们可以简单地选择一个绝对足够好的精度，比如7位。但是这样做的后果就在是一些并不需要这么好的精度的情况下产生了额外的无用运算。

考虑到开平方运算是迭代进行的，一个思路就是把迭代过程中的每个值都返回给调用者，并在调用者需要下一个更好的解的时候才进行下一步的计算，直到调用者满意为至。

当一个思路产生之后，无论你脑子里是否能立马闪现出具体实现方案，一个好的习惯就是在先计算机领域看看，有没有合适的技术与理念可以用来实现这个思路的，并和你自己的方案进行比较、甄选。然后再着手实现。思路的寻找往往是最难也是最有价值的部分。如果你总是想到什么就做什么，可能所有问题倒是也都可以解决，但是在这个过程中很难会有思路上、水平上的提升。把一年的工作经验用十年，就是这个意思。

了解到计算机系统中可以用来实现上面所述思路的概念至少有：

1. 流(Stream – a view of system structure)：具有的延时求值特性。
2. 按名调用(Call by name – an evaluation strategy)
3. 协程(Coroutine – a language feature)

Stream会更适合实现我们所需要的功能。它不依赖于语言，而且应用广泛。LINQ， Reactive Extensions都是基于Stream View设计的。

参考《SICP》上的流实现，不难构建出如下的Stream类和用之实现的Stream化的Sqrt函数：

type Stream<'a>(empty:bool, head:'a, tail:unit -> Stream<'a>) = 
  new(head:'a, tail:unit -> Stream<'a>) = Stream(false, head, tail)
  member this.Head = match empty with 
                      | false -> head
                      | true -> invalidOp "An empty stream has no head."
  member this.Tail = tail()
  member this.IsEmpty = empty
  member this.iterate action = 
    if not empty then 
      action this.Head; this.Tail.iterate action
  member this.iterateWhen action condition = 
    if not empty then
      action this.Head;
      if condition this.Head then
        this.Tail.iterateWhen action condition
  static member Empty = 
    Stream<'a>(true, Unchecked.defaultof<'a>, 
               fun unit -> invalidOp "An empty stream has no tail.")
  static member Map func (s:Stream<'a>) = 
    if not s.IsEmpty then 
      Stream(func(s.Head), fun () -> Stream.Map func s.Tail)
    else
      Stream.Empty
  static member Infinit producer head = 
    let rec stream = 
      fun () -> Stream(head, fun () -> stream() |> Stream.Map producer)
    stream()

let sqrt_stream x =
  let inline avg_damping f = fun x -> average x (f x)
  let getNext = avg_damping (fun v -> x * 1.0 / v)
  Stream.Infinit getNext 1.0

代码7

在这个版本的开平方函数中，没有任何精度设定，也就没有实现上的精确度损失。但是限于我目前的F#水平，我相信这个Stream的实现并不简洁，欢迎大家给出改进建议。（这遍文章里的F#代码大概占我写过的全部F#的50%，选F#只是因为我觉得照搬书上的Scheme代码和用我熟悉的语言解决这么简单的问题都很无聊。）

符号求导

在上面的牛顿法一节中，我们为了求解函数的导数而引入了一个新的精度假设值并导致了一定的精度损失。其根本原因是我们其实根本没有求解出导函数，而是利用了计算机强大的运算能力直接计算出指定点的导数值而已。

如果我们能把导函数方程先求解出来，并用这个导函数计算导数，就不会再有那个精度误差了。也就是说。我们希望能把方程$$f(x)=ax^2+bx+c$$当参数，并求解出$$f(x)=2ax+b$$这个导函数。

高级编程语言基本都自带表达式解析库。F#也不例外。所以符号化求导实现起来也并不复杂，至少不用自己写词法分析器。你当然可以自己写一个彰显实力，只是不太好读的代码多了也就不好管理了。但是需要避免的状况是：你自己写的一个，仅仅是因为你不知道语言里自带了个库，或是不屑于用。这种情况很常见，比如一个例子是:

let abs a = if a > 0 then a else -a

你看不出来问题在哪？那说的就是你。

在F#中实现符号求导，我们只需要把所需要用的如下求导法则：

• $$\frac{dc}{dx}=0$$
• $$\frac{dx}{dx}=1$$
• $$\frac{d(u+v)}{dx}=\frac{du}{dx}+\frac{dv}{dx}$$
• $$\frac{d(uv)}{dx}=u(\frac{dv}{dx})+v(\frac{du}{dx})$$

代码化成一组模式匹配：

let inline sum     e1 e2 = <@@ (%%e1: float) + (%%e2: float) @@>
let inline sub     e1 e2 = <@@ (%%e1: float) - (%%e2: float) @@>
let inline product e1 e2 = <@@ (%%e1: float) * (%%e2: float) @@>

let rec deriv (body:Expr) (p:Var) =
    match body with 
    | Var(x) ->
        if x = p then <@@ 1.0 @@> else <@@ 0.0 @@>
    | ValueWithName(x) ->
        <@@ 0.0 @@>
    | SpecificCall <@ (+) @> (_, _, [lArg;rArg]) ->
        sum (deriv lArg p) (deriv rArg p)
    | SpecificCall <@ (-) @> (_, _, [lArg;rArg]) ->
        sub (deriv lArg p) (deriv rArg p)
    | SpecificCall <@ (*) @> (_, _, [lArg;rArg]) -> 
        let pl = product lArg (deriv rArg p)
        let pr = product rArg (deriv lArg p)
        sum pl pr
    | _ -> invalidArg "body" ("Invalid body: " + body.ToString())

let derivE (exp:Expr) = 
    let (var, body) =
        match exp with 
        | Lambda(var, body) -> (var, body)
        | _ -> invalidArg "exp" ("Invalid exp: " + exp.ToString())
    let derivBody = deriv body var
    Expr.Lambda(var, derivBody)

let sqrt_v5 (y: float) = 
 let exp = <@ fun x -> x * x - y @>
 let f  = exp |>           EvaluateQuotation :?> ( float -> float )
 let Df = exp |> derivE |> EvaluateQuotation :?> ( float -> float )
 find_fix_point (fun x -> x - f x / Df x) 1.0

代码8

现实的项目中没有人会用符号求导技术来写个开平方函数。因为它比前面的版本要慢上1000倍。但是符号求导所带来的灵活性使得这个框架可以用来解决比开平方复杂得多的问题。

这个版本的代码没有使用流，如果需要，可以把它也实现成一个流。在软件工程领域，不同技术方案基本都可以合并而生成一个新的方式。而在非软件工程领域，这种合并常常会由于物理法则的约束而无法达成。这是软件工程复杂性的又一个体现。

综述

本文的主要意图是阐释，充实自身业务领域知识与计算机领域知识对于编写高质量代码的必要性，使用合适的方案解决一个问题能带来多方面的好处。这似乎是一句废话，但是其实这只是一个认识陷阱：人们，尤其是经验越丰富的人，一但认为自己理解了什么，或是利用起了他们的既往经验，他们就常常停止了探索，停止了学习，停止了思考。这对于软件开发是非常有害的。

但是另一方面，如上的发现问题、解决问题、发现新问题的循环是没有尽头的。具有钻研精神的程序员们总能在代码中发现问题并用各种方式尝试解决它，同时引入新的问题。专业的软件工程师能够根据最终目标和环境的不同知道应该往哪个方向走，并且到站了也知道停下来。在一个仅仅需要开平方的地方引入符号求导机制同样有害。

本文中的例子和解决问题的过程纯粹为阐释上述意图而刻意设定，请勿作编程指导，如果你能从中发现自己坠入的认识陷阱就再好不过了。

参考书目

Jerome H.Saltzer, M. Frans Kaashoek. Priciples of Computer System Design （《计算机系统设计原理》）

Harold Abelson, Gerald Jay Sussman, Julie Sussman. Structure and Interpretation of Computer Programs 2nd Edition （《计算机程序的构造和解释》）

Tomas Petricek, Jon Skeet. Real-world Functional Programming – With Examples in F# and C# （《C#与F#编程实践》）

Adrian Banner. The Calculus Lifesaver – All the tools you need to excel at Calculus （普林斯顿微积分读本）

Gerald M.Weinberg. The Psychology of Computer Programming （《程序开发心理学》）

谷强强 版权所有