從三角函數 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} 中間是否還有其他的梅森質數...

各種反直覺的項目

前陣子看到「The most counterintuitive facts in all of mathematics, computer science, and physics」這篇,講各種反直覺的項目,有空的時候拉一兩個來讀還蠻有趣的...

第一個提到的是 Homomorphic encryption。在密碼學的保護概念中,密文是不能被操作的,但我們可以透過重新設計密碼學系統,讓密文可以被運算:

1. It is possible to compute over encrypted data without access to the secret key: https://en.wikipedia.org/wiki/Homomorphic_encryption

2009 年的時候第一個 FHE (fully homomorphic encryption) 演算法被提出來後就被大量討論,這代表你可以把資料加密後丟上雲端,利用雲端的大量運算資源運算完後,再把資料拉回地端解密,得到運算後的結果。這算是 Lattice-based cryptography 的一個很重大的應用。

其他的主題也蠻有趣的,先丟著找機會翻...

在 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)。

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

證明圓周率 π 是無理數

前陣子在 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 (不得不說大家都很愛這味...)

Miles 換算 KM 的方式

Twitter 上看到很有趣的方式:

這邊可以這樣算是因為 1.609 跟黃金比率很接近,而 Fibonacci number 的也有黃金比率的特性,所以可以直接拿來用...