算 π 的 Chudnovsky algorithm

Mastodon 上面看到 Rob Pike 提到計算圓周率 π 的 Chudnovsky algorithm:(原始連結在 https://hachyderm.io/@robpike/112794329524261491)

照維基百科的說明,應該是目前收斂速度最快的演算法群之一,從近期的世界紀錄可以常常看到這個演算法:

It was used in the world record calculations of 2.7 trillion digits of π in December 2009, 10 trillion digits in October 2011, 22.4 trillion digits in November 2016, 31.4 trillion digits in September 2018–January 2019, 50 trillion digits on January 29, 2020, 62.8 trillion digits on August 14, 2021, 100 trillion digits on March 21, 2022, and 105 trillion digits on March 14, 2024.

維基百科條目的後面也有提到在算這些世界紀錄時用到的最佳化手段,甚至連 Python 實作都直接列出來了,可以看到遞迴的部分全部是整數運算,而最後轉換的時候也才用到一個整數的根號計算。

用 TeX 寫數學公式時常見的錯誤與解法

Hacker News 上看到「The Art of LaTeX: Common mistakes and advice for typesetting proofs (fanpu.io)」,原文在「The Art of LaTeX: Common Mistakes, and Advice for Typesetting Beautiful, Delightful Proofs」這,在講用 TeX 寫數學公式時常見的錯誤 (以及正確的寫法)。

裡面給了蠻多還蠻實用的提醒,像是第一點提出來的這兩個差異:

以及:

上面是用最單純的 {( 寫的,要達到下面比較正確的效果,不應該用 \bigl 這種直接指定括弧大小的用法,而是用 \left( 這種告訴 TeX 排版引擎這邊的意義,然後算出需要拿多大的括弧出來用。

畢竟 TeX 是 Knuth 老人家弄出來的東西,你會看到很多設計是需要給 TeX 排版引擎正確的定義或是資訊,剩下的他就會幫你處理...

Python 上的 reals 套件 (需要 3.10+ 以上才能裝)

看到「A lightweight python3 library for arithmetic with real numbers.」這個有趣的 Python 延伸套件,可以用他進行高精度的實數運算...

一開始在 Python 3.9 環境裝,結果就跳出需要 3.10+ 的環境,想了一下,開了一個 Docker container 裝 pyenv 來測,測過以後覺得還蠻有趣的,看起來之後把預設環境變成 3.10+ 應該會裝起來用...

這個 reals 的重點在於保證顯示數字的正確性:

It allows you to compute approximations to an arbitrary degree of precision, and, contrary to most other libraries, guarantees that all digits it displays are correct.

目前支援的常數與操作有這些:

Constants: pi, e, phi
Functions related to powers: sqrt, exp, log
Operators: negation, addition, subtraction, multiplication, division, powers
Trigonometric functions: sin, sinh, csc, csch, cos, cosh, sec, sech, tan, tanh, cot, coth

用法的部份,先把 reals 拉進來:

>>> from reals import sqrt

然後用法算直覺:

>>> sqrt2 = sqrt(2)
>>> sqrt2
<reals._real.Real object at 0x10d182560 (approximate value: 1.41421)>
>>> sqrt2.evaluate(10)
'1.4142135624'
>>> '{:.10f}'.format(sqrt2)
'1.4142135624'
>>> sqrt2.to_decimal(10)
Decimal('1.4142135624')

不過作者有提到效能沒有處理到很好,所以應該是拿來快速做一些運算得到結果而已。

用 GPT-3 解讀程式碼

Hacker News 上看到的方法,Simon Willison 試著把程式碼餵進 GPT-3,然後問 GPT-3 程式碼的意思,看起來答的還不錯:「Using GPT-3 to explain how code works」,對應的討論 (包括 Simon Willison 的回應) 則可以在「Using GPT-3 to explain how code works (simonwillison.net)」這邊看到。

第一個範例裡面可以解讀 regular expression,雖然裡面對 (?xm) 的解讀是錯的,但我會說已經很強了...

第二個範例在解釋 Shadow DOM,看起來也解釋的很不錯...

第三個範例回來原來產生程式碼的例子,拿來生 SQL 指令。

後面的 bonus 題目居然是拿來解釋數學公式,他直接丟 TeX 文字進去要 GPT-3 解釋柯西不等式 (Cauchy–Schwarz inequality)。這樣我想到以前高微作業常常會有一堆證明題,好像可以丟進去要 GPT-3 給證明耶...

從三角函數 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 跟俄羅斯軍方的關係也是很知名,這些東西大概要到十來年後才會知道...