時隔3年，再做雙倍超立方數的題目，這次用Lisp

09年的時候，有一陣子大家都在做網易的“有道難題”里面的雙倍超立方數問題。我當時看到題目后也隨手用 Python 做了一下，但當時那個解法是最笨的窮舉法，所以除了練習了一次 Python 的基本語法之外，也沒有什么特別的收獲了。

最近又把 SICP 翻出來看到了流這一部分。在我發自內心的贊嘆那些用流描述的數列在形式上是何等優美的同時，我想實際測量一下應用了這么多復雜的流的函數操作之后，程序的執行效率如何。而雙倍超立方數正好是一個適合拿來練習的題目。

這次我用的開發環境是 DrRacket. 這個東西的前身其實叫做 DrScheme. 是 scheme 的一種開發環境，在 Windows 上比較容易安裝。而因為 SICP 是用 scheme 語言描述的，所以用這個環境來開發就是理所當然的事了。

雙倍超立方數的表述形式為：

c = a^3 + b^3

對給定的c，如果有且只有兩組正整數(a, b)使上式成立，那么c就稱為雙倍超立方數。

基於對稱性，可以設 a <= b.

下面我們來觀察一下從 1^3 + 1^3 開始的兩個自然數的立方和組成的數列是怎樣的：

容易觀察到，上述圖中的 (a, b) 組合是沒有重復的。也就是說，這代表了所有 (a, b) 的可能的組合情況。雙倍超立方數題目的要求是，比如我們需要求出不超過n（比如100萬）的所有雙倍超立方數的個數。如果我們能夠對上述流中的數字在僅進行一趟從小到大的遍歷，在遍歷過程中判斷出其中哪些立方和出現了2次，則可以有效的減少遍歷的次數，得到答案。因此問題的關鍵轉化為如何有效的按立方和從小到大的次序遍歷上述矩陣。

可以看到，矩陣中橫向、縱向、以及從左下到右上的每一條對角線方向上的數字的立方和都是單調遞增的。但是很難直觀的對整個矩陣找到這樣一條單調遞增的遍歷路徑。

進一步觀察計算結果，可以發現，沿對角線方向的數字序列，兩兩數字之間的間隔最小。可以把每一條對角線上的數字序列定義為一個流，然后所有對角線組成的流進行合並（merge）。合並的方式是，依次取每個流的頭部元素互相比較大小，按從小到大的次序組合進結果流中。因為按對角線方向元素間隔相對來說最小，在流進行合並的時候，能夠連續的從同一個流中獲取多個連續的元素，那么在合並時需要不斷切換從不同的流中獲取下一個元素的情況應該就比較少一些，如果實現得當，代碼的局部性 (locality) 也許會比較好，也許合並算法上也可以有進一步優化的可能。（當然這個只是我目前的推斷，還沒有證實）

按照這個思路，整個問題解法的步驟如下：

1. 得到一個按大小順序排列的兩個數立方和的流（通過多個對角線的數字流 merge 得到）
在合並的同時，對每個和數出現的次序進行 count. 同時在結果流里面保留每個和出現的次數；
2. 對步驟1得到的流，進行過濾操作 (stream-filter)，僅保留上述出現次數恰好為 2 的那些結果；
3. 對步驟2得到的流計算長度 (stream-length），結果就是 <= N 的雙倍超立方數的個數。

測試下來發現，N=100萬的時候，執行速度還是非常快的，跟 Python 版解法的速度差不多。可以看到，利用流，以及作用於流的各種函數（filter, map, merge 等等操作）作為一種高級的數據結構抽象，可以從另一個角度理解問題。而因為流是惰性求值的，這個特性也使得運算的時間和空間開銷等方面都處於一個可接受的范圍。

下面是我調試通過的代碼。其中還有很大的優化空間，比如很簡單的就可以對同一個數字反復求立方的操作進行記憶，避免重復運算的時間開銷。

代碼的上半部分都是用來操作流的一些通用的幫助函數，有些是 SICP 書里的，另外有幾個是我自己寫的。雙倍超立方數問題的相關實現可以直接看代碼的后半部。

執行這個代碼可以正確的打印出100萬以下雙倍超立方數的數目：43.

下面代碼中關於流的一些輔助函數的解釋，可以參考我上一篇隨筆：Racket, SICP stream learning

#lang racket
;***********************************************
; 雙倍超立方數問題求解。
; 
; 木野狐(Neil Chen)@博客園
;
; 開發環境：DrRacket v5.2.1
; 2012-05-15
;***********************************************

;=========== 准備工作（流操作相關函數）開始 =====
(define (stream-car stream) (car stream))

(define (stream-cdr stream) (force (cdr stream)))

(define-syntax cons-stream
  (syntax-rules ()
    [(cons-stream x y) (cons x (delay y))]))

(define the-empty-stream '())

(define (stream-null? stream)
  (null? stream))

(define (stream-filter pred stream)
  (cond ((stream-null? stream) the-empty-stream)
        ((pred (stream-car stream))
         (cons-stream (stream-car stream)
                      (stream-filter pred (stream-cdr stream))))
        (else (stream-filter pred (stream-cdr stream)))))

(define (stream-ref s n)
  (if (stream-null? s) the-empty-stream
      (if (= n 0)
          (stream-car s)
          (stream-ref (stream-cdr s) (- n 1)))))

(define (stream-map proc . argstreams)
  (if (stream-null? (car argstreams))
      the-empty-stream
      (cons-stream (apply proc (map stream-car argstreams))
                   (apply stream-map 
                          (cons proc (map stream-cdr argstreams))))))

(define (stream-enumerate-interval low high)
  (if (> low high)
      the-empty-stream
      (cons-stream low
                   (stream-enumerate-interval (+ low 1) high))))

(define (enumerate-interval low high)
  (if (> low high)
      '()
      (cons low
            (enumerate-interval (+ low 1) high))))

(define (stream-length stream)
  (define (iter s acc)
    (if (stream-null? s) acc
        (iter (stream-cdr s) (+ 1 acc))))
  (iter stream 0))
;=========== 准備工作（流操作相關函數）結束 ======

(define (integers-starting-from n)
  (cons-stream n (integers-starting-from (+ n 1))))

(define integers
  (integers-starting-from 1))

(define (cube n)
  (* n n n))

(define cube-numbers
  (stream-map cube integers))

(define (merge-with-count s1 s2)
  (cond ((stream-null? s1) s2)
        ((stream-null? s2) s1)
        (else
         (let ((s1car (stream-car s1))
               (s2car (stream-car s2)))
           (let ((e1 (if (list? s1car) (car s1car) s1car))
                 (c1 (if (list? s1car) (cadr s1car) 1))
                 (e2 (if (list? s2car) (car s2car) s2car))
                 (c2 (if (list? s2car) (cadr s2car) 1)))                         
             (cond ((< e1 e2) (cons-stream (list e1 c1) 
                                           (merge-with-count (stream-cdr s1) s2)))
                   ((> e1 e2) (cons-stream (list e2 c2)
                                           (merge-with-count s1 (stream-cdr s2))))
                   (else (cons-stream (list e1 (+ c1 c2)) (merge-with-count (stream-cdr s1) (stream-cdr s2))))))))))

(define (stream-merge-with-count . argstreams)
  (if (= 0 (length argstreams))
      the-empty-stream
      (merge-with-count (car argstreams)
                        (apply stream-merge-with-count (cdr argstreams)))))

; 得到第 n 條對角線上的數字的流
(define (diag n)
  (stream-map (lambda (i)   
                (+ (cube i)
                   (cube n)))
              (stream-enumerate-interval 1 n)))

(define (new-double-cubic-numbers-below n)
  (let ((max (floor (expt (- n 1) (/ 1 3)))))
    (stream-map car
                (stream-filter 
                 (lambda (s)
                   (and (list? s)
                        (= 2 (cadr s))
                        (not (> (car s) n))))
                 (apply stream-merge-with-count 
                        (map diag 
                             (enumerate-interval 1 max)))))))

(stream-length (new-double-cubic-numbers-below 1000000))