從三角函數 cosine 的實做問題學一些週邊知識...

前幾天在 Hacker News 上看到「Implementing Cosine in C from Scratch (2020) (austinhenley.com)」這篇 2020 的文章,原文是「Implementing cosine in C from scratch」,裡面內在講自己刻三角函數的 cosine 所遇到的一些嘗試。

cosine 是很基本的函數,所以可以使用的地方很多。另外一方面,也因為他不是那麼直覺就可以實做出來,在現代的實做裡面其實藏了超多細節...

不過真的有趣的是在翻 Hacker News 上的討論時陸陸續續翻其他的資料看到的知識。

第一個看到的是 Intel 對於 FPU-based 指令集內的 FSIN 因為 π 的精度不夠而導致誤差超大 (尤其是在 0 點附近的時候):「Intel Underestimates Error Bounds by 1.3 quintillion」,然後 AMD 是「相容」到底,所以一樣慘:「Accuracy of FSIN and other x87 trigonometric instructions on AMD processors」。

這個就是有印象,但是太久沒有提到就會忘記...

第二個是 musl libc 裡的 cosine 實做 (看註解應該是從 FreeBSD 的 libc 移植過來的?):「__cos.c\math\src」與「cos.c\math\src」(話說 cgit 在 html 內 title 的內容對路徑的表達方式頗有趣,居然是反過來放...)。

拆開的部份是先將範圍限制在 [-\pi/4, \pi/4] 後 (這個部份看起來是透過 __rem_pio2.c 處理),再丟進公式實際運算。

另外帶出來第三個知識,查資料的時候翻到 binary64 (這也是 C 語言裡面的 double) 與 binary128 的差異:

而大家很常拿來惡搞的 double double 則是利用兩個 double 存放,形式是 v = head + tail,利用不同的 exponent 表示來不同部份的值,以提高經度:

A common software technique to implement nearly quadruple precision using pairs of double-precision values is sometimes called double-double arithmetic.

不過這樣的精確度只能到 106 bits,雖然跟 binary128 能達到的 113 bits 相比低了一些,但在大多數的情況下也還算夠用:

Using pairs of IEEE double-precision values with 53-bit significands, double-double arithmetic provides operations on numbers with significands of at least[4] 2 × 53 = 106 bits (...), only slightly less precise than the 113-bit significand of IEEE binary128 quadruple precision.

Python 裡使用超過 Double Precision 的運算

Hacker News 上爬到的,是一篇 2019 的文章:「When Double Precision Is Not Enough (adambaskerville.github.io)」,原文在「T>T: When Double Precision is Not Enough」。

作者是拿矩陣 (matrix) 的運算當例子,遇到了 double precision 造成的計算誤差問題:

There is no error with the program; this discrepancy is caused by a loss of numerical accuracy in the eigenvalue calculation due to the limitation of hardware double precision (16-digit).

解法是用 mpmath 增加精度,算是一種暴力解,到要注意計算會慢很多:

Note that this library is incredibly slow for large matrices, so is best avoided for most applications.

另外在 Hacker News 的討論串裡面看到個有趣的東西:「Herbie: Automatically Improving Floating Point Accuracy」這個專案,你把公式丟進去,Herbie 會試著提供等價公式來維持精度,像是 \sqrt{x+1} - \sqrt{x} = 1 / ( \sqrt{x+1} + \sqrt{x} ) 這種東西。

半自動化幫你改善...

在電腦上計算二進制圓周率 π 的公式

Hacker News 的最新列表上面看到的演算法,用來計算圓周率 π 的公式:「Bellard's formula」。

\pi = \frac1{2^6} \sum_{n=0}^\infty \frac{(-1)^n}{2^{10n}} \, ( -\frac{2^5}{4n+1} - \frac1{4n+3} + \frac{2^8}{10n+1} - \frac{2^6}{10n+3} - \frac{2^2}{10n+5} - \frac{2^2}{10n+7} + \frac1{10n+9} )

這個公式特別的地方在於 \frac{1}{2^{10n}} 的部份,這使得整個公式可以很迅速的算出某個二進制位置上的值 (要做一些條件判斷,相較於十進制轉二進制的方式快多了),而且可以馬上想到一些平行運算的方式...

另外一個讓我注意到的點是,這居然是 Fabrice Bellard 丟出來的公式,真的是到處看到他的蹤跡耶...

第 48 個梅森質數確認

在「Great Internet Mersenne Prime Search - PrimeNet」這邊看到 GIMPS 確認了 M_{57885161} (也就是 2^{57885161}-1) 是第 48 個梅森質數:

All tests smaller than the 48th Mersenne Prime, M(57885161), have been verified

第 47 個梅森質數是 M_{43112609},但一直都沒有確認完到 M_{57885161} 中間是否還有梅森質數,這次宣佈全部都算完,而且雙重確認過,所以可以確認第 48 個梅森質數。

下一個要確認的是第 49 個梅森質數,目前已知的下一個是 M_{74207281},也就是要看到 M_{74207281} 中間是否還有其他的梅森質數...

在 Hacker News 上看到選擇公理

沒想到會在 Hacker News 的首頁上看到這麼硬核的主題,選擇公理 (Axiom of choice,通常縮寫成 AC):「What is the Axiom of Choice?」,對應的討論在「What is the Axiom of Choice? (jaydaigle.net)」。

出自「xkcd: Set Theory

應該是大一教集合論的時候學到的,算是一個非常重要的公設,雖然的確有些數學系統是可以假定 AC 不成立,但用起來會不太好用,主要是因為「對於集合 S,取出任意一個元素」這類用法太常出現,在沒有 AC 的情況下這件事情就不一定能操作了...

我們目前常用的數學一般是建立在 Zermelo-Fraenkel Set Theory (ZF) 這個公理系統加上 AC,簡寫變成 ZFC。而 AC 在集合論常常會被拿出來說明,主要還是因為在歷史上花了不少力氣才證明 ZF 與 AC 的相對協調性 (ZF 與 AC 不衝突),以及 ZF 與 AC 獨立性 (ZF 無法推導出 AC)。

有了 AC 後就會再解釋連續統假設 (Continuum hypothesis,簡稱 CH),也就是 \mathbb{N}2^{\mathbb{N}} 之間存不存在一個集合 S 使得 |\mathbb{N}| < |S| < |2^{\mathbb{N}}|

然後再打臉一次,說明 ZFC 與 CH 的協調性 (ZFC 與 CH 不衝突),與獨立性 (ZFC 無法推導出 CH)。

當時學的時候的確是很頭痛,不過現在回頭看倒是覺得很有趣:在數學上你可以證明「某個敘述無法被證明」,這點應該是以前沒想過的...

Kaspersky Password Manager 的漏洞

Hacker News Daily 上看到「Kaspersky Password Manager: All your passwords are belong to us」這篇,講 Kaspersky Password Manager (KPM) 嚴重的安全漏洞,另外在 Hacker News 上的討論「Kaspersky Password Manager: All your passwords are belong to us (ledger.com)」也有提到一些有趣的東西。

標題的 All your passwords are belong to us 是出自「All your base are belong to us」這個梗的變形。

這包安全問題主要的原因是因為 KPM 沒有使用 CSPRNG,而且也沒有正確 seed,所以極為容易被猜出密碼本身。

KPM 的 Web 版使用了 Math.random(),在各家瀏覽器主要是用 xorshift128+ 實做 Math.random(),作者沒有針對這塊再花時間研究,但很明顯的 Math.random() 不是個 CSPRNG:

The underlying PRNG used by Chrome, Firefox and Safari for Math.random() is xorshift128+. It is very fast, but not suitable to generate cryptographic material. The security consequences in KPM has not been studied, but we advised Kaspersky to replace it with window.crypto.getRandomValues(), as recommended by the Mozilla documentation page previously mentioned.

Note: Math.random() does not provide cryptographically secure random numbers. Do not use them for anything related to security. Use the Web Crypto API instead, and more precisely the window.crypto.getRandomValues() method.

而桌機版則是用了 MT19937,理論上取得 624 bytes 的輸出後就可以重建整個 PRNG 的內部狀態 (於是就可以預測後續的 output),但這代表你要知道其他網站的密碼,這點其實有點困難。

但作者發現 KPM 在產生 MT19937 的 seed 只跟時間有關,超級容易被預測:

So the seed used to generate every password is the current system time, in seconds. It means every instance of Kaspersky Password Manager in the world will generate the exact same password at a given second.

於是可以直接暴力解出所有的可能性:

The consequences are obviously bad: every password could be bruteforced. For example, there are 315619200 seconds between 2010 and 2021, so KPM could generate at most 315619200 passwords for a given charset. Bruteforcing them takes a few minutes.

Hacker News 上有不少陰謀論的討論,像是:

Getting some DUAL_EC prng vibes.

Insert Kaspersky owned by Russia intelligence conspiracy here...

另外 Kaspersky 跟俄羅斯軍方的關係也是很知名,這些東西大概要到十來年後才會知道...

證明圓周率 π 是無理數

前陣子在 Michael Penn 教授的 YouTube 頻道上看到的證明:

這個證明方式在其他人的 YouTube 頻道也有做過,像是 MindYourDecisions 在 2019/3/14 (Pi Day) 發表的影片:

這兩篇使用的方法是 Ivan M. Niven 在 1946 年提出來的方法 (發表在 1947 年 AMS 的期刊,第 53 期的 Bulletin of the American Mathematical Society 上):「A simple proof that π is irrational」。

先不討論證明方法本身,我把這份 PDF 印出來慢慢看,當年用打字機打的字型讀起來超有味道的... (這份 PDF 看起來是掃描檔)

從 PDF 可以看出證明超級短,只有一頁,但畢竟這是大師丟出來的證明,裡面其實省了很多步驟。這接步驟不難推導,但是是屬於考試作答時不能省略的部份,如果加上去的話應該會兩到三頁 (還是非常短)。

整個證明的過程很巧妙的設計了兩個函數搭配反證法。裡面有用到微積分,但只用到最基本的微積分,其中微分的部份用到的公式就這三個,其中是三角函數的微分公式:

\frac{d}{dx} sin(x) = cos(x)\newline\newline\frac{d}{dx} cos(x) = -sin(x)

然後是兩個函數相乘後的微分公式:

\frac{d}{dx} f(x)g(x) = f'(x)g(x) + f(x)g'(x)

積分的部份只用了定義的部份 (主要是在微分的部份證完)。

然後證明是圍繞在「整數」上面做文章,夾擠出一個不可能的情況,因此得到矛盾,證明了 π 是無理數。

整個證明可以在高三範圍 (如果現在高三還有教微積分的話) 或是大一教過微積分之後的程度證完,作者講的 "Simple" 主要是出自這邊...

四個機率問題

在「Some Useful Probability Facts for Systems Programming」這邊看到的四個機率問題,而且都有接近解:

  • 每一次有 \frac{1}{N} 成功的機率,跑 N 次後最少成功一次的機率。
  • N 個球隨機丟到 N 個籃子後有空籃子的機率。
  • N 個數字在隨機排序後有數字不改變位置的機率。
  • N 種 coupon 平均隨機出現,要抽中一套 N 種都有的 coupon 需要抽幾次的期望值。

會想要寫這篇是因為發現最後一題就是「大人買」抽卡問題的簡化版本 (這邊簡化成機率相同,一般的情境下應該是不同的)。

可以看到歐拉常數 e 與自然對數 \ln{} (i.e. 以 e 為底的對數) 滿天飛 XDDD

用程式解數學邏輯問題...

Hacker News Daily 上看到的數學邏輯問題:「“Which answer in this list is the correct answer to this question?”」。

問題是這樣:

Which answer in this list is the correct answer to this question?

  • All of the below.
  • None of the below.
  • All of the above.
  • One of the above.
  • None of the above.
  • None of the above.

accepted 的那個是推演的答案,但最高分的那個是寫程式窮舉 XDDD (不得不說大家都很愛這味...)

二戰時德國坦克製造速度的估算問題

看到「The German Tank Problem」這篇在講二戰很有名的統計應用。這個主題在中文的維基百科寫得還蠻完整的,讀起來應該會更快一些:「德國坦克問題」:

在統計學理論的估計中,用不放回抽樣來估計離散型均勻分布最大值問題中著名的德國坦克問題(英語:German tank problem),它因在第二次世界大戰中用於估計德國坦克數量而得名。

如同上面所說的,這個方法是因為估算的準確度極高而知名:

對坦克車輪的分析產生了對使用中的車輪模具數量的估計。在與英國車輪製造商討論過後,他們估計了這麼多的模具可以生產多少車輪,進而是每個月可生產的坦克數量。對兩輛坦克(每輛32個車輪,總計64個車輪)車輪的分析的結果是1944年2月的生產數量估計在270左右,大大超出此前預期。

德國戰後公布的記錄顯示,1944年2月一個月的生產量是276輛。統計方法結果的精確度是常規情報收集方法所遠遠不能達到的,而「德國坦克問題」這個詞也成為了這種統計分析問題的標誌。

而且之後被拿來推敲經典的 Commodore 64 的數量也還蠻準的:

該公式在非軍事中也有使用,如估計Commodore 64計算機的總數,其結果(1.25億)與官方數字相當匹配。