Rotating Doughnut Using Matrix Multiplication

My copy of Escher’s artwork got some attention and quite some people asked me, “what kind of programmer are you?”, well, this Youtube video explains it well.

Actually, Joma’s video basically depicted me when in real dev mode. Upon watching this video, I couldn’t wait to re-implement it in Python. However, the initial version is basically a trans-coding of the original C code and was pretty hard to read. So here is a version that actually uses Matrix multiplication in Python (In the mean time, I still have to fix the syntax highlighter…):


import time

import numpy as np

screen_size = 40
theta_spacing = 0.07
phi_spacing = 0.02

illumination = np.fromiter(".,-~:;=!*#$@", dtype="<U1")

A = 1
B = 1
R1 = 1
R2 = 2
K2 = 5
K1 = screen_size * K2 * 3 / (8 * (R1 + R2))


def render_frame(A: float, B: float) -> np.ndarray:
"""
Returns a frame of the spinning 3D donut.
Based on the pseudocode from: https://www.a1k0n.net/2011/07/20/donut-math.html
"""
cos_A = np.cos(A)
sin_A = np.sin(A)
cos_B = np.cos(B)
sin_B = np.sin(B)

output = np.full((screen_size, screen_size), " ") # (40, 40)
zbuffer = np.zeros((screen_size, screen_size)) # (40, 40)

cos_phi = np.cos(phi := np.arange(0, 2 * np.pi, phi_spacing)) # (315,)
sin_phi = np.sin(phi) # (315,)
cos_theta = np.cos(theta := np.arange(0, 2 * np.pi, theta_spacing)) # (90,)
sin_theta = np.sin(theta) # (90,)
circle_x = R2 + R1 * cos_theta # (90,)
circle_y = R1 * sin_theta # (90,)

x = (np.outer(cos_B * cos_phi + sin_A * sin_B * sin_phi, circle_x) - circle_y * cos_A * sin_B).T # (90, 315)
y = (np.outer(sin_B * cos_phi - sin_A * cos_B * sin_phi, circle_x) + circle_y * cos_A * cos_B).T # (90, 315)
z = ((K2 + cos_A * np.outer(sin_phi, circle_x)) + circle_y * sin_A).T # (90, 315)
ooz = np.reciprocal(z) # Calculates 1/z

xp = (screen_size / 2 + K1 * ooz * x).astype(int) # (90, 315)
yp = (screen_size / 2 - K1 * ooz * y).astype(int) # (90, 315)

L1 = (((np.outer(cos_phi, cos_theta) * sin_B) - cos_A * np.outer(sin_phi, cos_theta)) - sin_A * sin_theta) # (315, 90)
L2 = cos_B * (cos_A * sin_theta - np.outer(sin_phi, cos_theta * sin_A)) # (315, 90)
L = np.around(((L1 + L2) * 8)).astype(int).T # (90, 315)

mask_L = L >= 0 # (90, 315)
chars = illumination[L] # (90, 315)

for i in range(90):
mask = mask_L[i] & (ooz[i] > zbuffer[xp[i], yp[i]]) # (315,)

zbuffer[xp[i], yp[i]] = np.where(mask, ooz[i], zbuffer[xp[i], yp[i]])
output[xp[i], yp[i]] = np.where(mask, chars[i], output[xp[i], yp[i]])

return output


def pprint(array: np.ndarray) -> None:
"""Pretty print the frame."""
print(*[" ".join(row) for row in array], sep="\n")


if __name__ == "__main__":
for _ in range(screen_size * screen_size):
A += theta_spacing
B += phi_spacing

# clear screen using ANSI control code
print("\x1b[H")
time.sleep(0.05)
pprint(render_frame(A, B))

 

Chapter 02 The First Algorithm

摩西快速学习了算术和几何……
这些知识他是从埃及人那里学到的。埃及人研习数学和其它学问。

《形而上学》亚里士多德

算法 是用以完成某项计算任务的有限步骤。算法跟计算机编程有紧密的联系,以至于大多数知道这个词的人以为算法的使用是从计算机科学开始。但实际上,人们使用算法已经有几千年的历史。数学中充满了算法,有一些我们到今天还在使用。学生们学习的长加法就是一种算法。

尽管使用算法的历史很长,算法这个概念并不是一开始就存在,它是后来发明出来的。尽管我们不知道是谁首先发明了这个概念,却明确知道有些算法至迟在4000多年前的古埃及就存在了。

 

古埃及文明的核心是尼罗河,因为其农业有赖于尼罗河泛滥带来的肥沃泥土。问题是,尼罗河一发水,用来划分权属的界标也都被冲掉了。埃及人惯用绳子测量长度,因此发展出一套方法,可以根据书面记录重建界标。一群制定的祭司专司此事,他们学过这些数学技术,因此成为“拉绳子的人(rope-strecher)”。希腊人后来会把他们叫做几何学家(geometers),意思是“土地丈量者”。

很不幸我们并没有多少古埃及数学知识的文字记录,当时的数学文件仅有两篇留存至今。我们关心的一个,叫做Rhind Mathematical Papyrus,名字来自于19世纪在埃及买到它的苏格兰收藏者。这篇文献成于约公元前1650年,抄写者名叫Ahmes。文献中有一系列算术和几何问题,还有一些辅助计算的表格。其中包含一个快速乘法技术和一个快速除法技术,是最早的有记录的算法。我们首先来看一看这个快速乘法算法,(我们很快会看到)它至今仍然是重要的计算技术。

2.1 埃及乘法

跟其它古文明一样,埃及人使用的数字系统中没有数位概念,也没有零。因此,乘法非常复杂,只有少数专家能够掌握。(想象一下只使用罗马数字进行两位数乘法)

首先,什么是乘法?粗略地说,乘法就是“把一个数跟它自身反复相加”。严格一些,我们可以首先把乘法分为两种情况:乘以1,或者乘以一个大于1的数。

乘数为1的乘法的定义是:

$$1a=a$$                                   (2.1)

然后,我们考虑如何计算比我们已经计算出来的再多一倍的情形。有些读者可以意识到这是一个递推,我们会在后面正式使用这个技术。

$$(n+1)a=na+a$$                 (2.2)

计算na的厂家做法是把n个a加起来。但是,在n很大的时候,这可能非常繁琐,因为需要进行n-1次加法。使用C++,这个算法是这样的:

int multiply0(int n, int a){
    if (n==1) return a;
    return multiply0(n-1, a)+a;
}

这样两行代码对应于(2.1)和(2.2)。a和n都必须是正数,古埃及当时的约定也是如此。

Ahmes所描述的算法,这个算法在古希腊被称为“埃及乘法”,而很多现代作者则称之为“俄罗斯农民算法”,建立在如下洞察上:

$$ 4a=((a+a)+a)+a=(a+a)+(a+a) $$

这个优化利用了加法的结合律:

$$a+(b+c)=(a+b)+c $$

这样,我们只需要计算一次a+a,减少了计算次数。

核心思想在于,不断地把n取半,同时把a加倍,使用数的2的次方倍相加来计算乘法结果。在当时,人们不会使用变量a和n来描述算法,而是由作者给出示例,然后说“其它数字同理”。Ahmes也是这么做的,他使用下面的表格来演示这种算法,n=41,a=59:

1        ♦         59
2                  118
4                  236
8        ♦        472
16                944
32      ♦      1888

每个条目的最左边是2的幂,右边是前一条目的数字的倍数(数字加倍相对容易计算)。有标记的值对应41的二进制表示中数值为1的位。这个表格的意思是说:

$$41*59=(1*59)+(8*59)+(32*59)$$

其中右边的每一个乘积可以由对59加倍响应的次数而得到。

这个算法需要检查n是奇数还是偶数,所以我们可以推断古埃及人知道它们的区别,尽管我们没有直接的证据,但是古希腊人,他们自称从古埃及人那里学到数学,显然知道奇偶的区别。下面是他们的定义方法用现代数学标记表达出来:

$$n=n/2+n/2$$    说明n是偶数
$$n=(n-1)/2+(n-1)/2+1$$  说明n是奇数

我们也将使用这个特点:

odd(n) 意味着 half(n)=half(n-1)

下面是埃及乘法的C++实现:

int multiply1(int n, int a){
    if (n==1) return a;
    int result=multiply1(half(n), a+a);
    if (odd(n)) result=result +a;
    return result;
}

odd(x)很容易实现,测试x的最低位就可以了,half(x)则可以通过对x进行一次右移实现:

bool odd(int n) { return n&0x1;}
int half(int n) { return n>>1; }

multiply1需要进行多少次加法运算?每次调用这个函数,我们需要在a+a的地方进行一次加法,因为我们是不断对n取半,我们需要调用这个函数logn次。在某些时候,我们还需要在result+a的地方进行加法运算,因此加法运算的总次数为:

$$\#(n)=logn + v(n)-1$$

其中v(n)是n的二进制表达中1的个数(population count或者pop count)。由此我们把一个O(n)算法优化成了O(logn)算法。

这个算法已经最优了吗?不见得。比如,如果我们要乘以15,使用之前的公式,我们需要

$$\#(15)=3+4-1=6$$

进行6次加法。但是我们使用下面的办法,就只需要5次加法运算:

int multiply_by_15(int a){
    int b=(a+a)+a;    //b == 3*a
    int c+b+b;        //c == 6*a
    return (c+c)+b;   //12*a + 3*a
}

这样一系列加法运算又称为加法链。这里我们发现了15的最优加法链,尽管如此,Ahmes的算法在多数情况下已经足够好了。

练习2.1

为n<100的数找到最优加法链

有读者可能已经注意到其中有些乘法运算可以更快,只要我们在第一个参数大于第二个参数的情况下交换一下参数顺序就可以了(比如,3*15会比15*3更容易计算)。的确如此,而且古埃及人也知道这一点。但是我们不会在这里进行这样的优化,因为到第七章我们会发现,我们最终要把算法推广到不同的数据类型。在这种情况下,参数的顺序是不能交换的。

2.2 改进算法

函数multiply1的加法运算次数还可以,但是它还要进行logn次递归调用。鉴于函数调用的开销很高,我们希望能够优化程序以避免这个开销。

我们将要利用的一个原则是:It is often easier to do more work rather than less. 特别地,我们准备计算

r+na

其中r是加法计算的中间结果,也就是说,我们要实现一个乘法累计函数,而不仅仅是乘法函数。这种思维方式不仅仅适用于程序编写,也适用于硬件设计和数学——比如在数学中,时常有一般性判断比特定判断更容易证明的情形。

这里是我们的乘法累计函数:

int mult_acc0{int r, int n, int a) {
    if (n==1) return r+a;
    if (odd(n)){
        return mult_acc0(r+a, half(n), a+a);
    }else{
        return mult_acc0(r, half(n), a+a);
    }
}

这里有一个不变形式:r+na=r0+n0a0,其中r0、n0和a0是r、n和a这些变量的初始值。

我们可以通过简化递归获得进一步优化。注意在上面的两行递归调用代码只有第一个参数不同,这样对于奇数情况和偶数情况分别进行递归,不如我们在进行递归前根据情况修改第一个参数的值:

int mult_acc1(int r, int n, int a) {
    if (n==1) return r+a;
    if (odd(n)) r=r+a;
    return mult_acc1(r, half(n), a+a);
}

现在我们的函数是尾递归函数——递归仅仅发生在返回值里面——我们马上会利用这一特点。

我们注意到两点:

  • n=1的情况很少发生;
  • 如果n是偶数,则完全没有必要判断它还是不是1.

因此,我们可以首先判断n的奇偶性,这样把n与1进行比较的次数降低一半:

int mult_acc2(int r, int n, int a) {
    if (odd(n)) {
        r=r+a;
        if (n==1) return r;
    }
    return mult_acc2(r, half(n), a+a);
}

有些程序员认为编译器可以为我们做这样的优化,其实不太可能。编译器不能把一种算法优化成另外一种算法。

到目前为止都很好,不过我们的目标是避免递归以节省函数调用的开销。要达到目的,如果我们的函数是完全尾递归就好了。

定义2.1 一个完全尾递归函数是一个所有递归调用的形式参数都跟函数本身一致的函数。

同样地,我们可以通过在递归调用之前给变量重新赋值来达到这个要求:

int mult_acc3 (int r, int n, int a) {
    if (odd(n)) {
        r=r+a;
        if (n==1) return r;
    }
    n=half(n);
    a=a+a;
    return mult_acc3(r,n,a);
}

到了这里,把它转化成一个循环就很简单了,只需要使用一个while(true)结构替换掉尾递归就可以了:

int mult_acc4(int r, int n, int a) {
    while (true) {
        if (odd(n)) {
            r=r+a;
            if (n==1) return r;
        }
        n=half(n);
        a=a+a;
    }
}

使用这个全新的乘法累加函数,我们可以写一个新版的乘法函数:

int multiply2(int n, int a) {
    if (n==1) return a;
    return mult_acc4(a, n-1, a);
}

注意我们在调用mul_acc4的时候预先把累加结果置为a,这样就减少了一次循环。

效果很不错,除了在n是2的幂的时候。由于我们第一次调用就把n减去1,这样mult_acc4函数被调用的时候第二个参数是一个二进制表达是全1的数,这对于我们的算法是最坏情况。因此,我们要在n是偶数的情况下做一些预处理,以避免这种情形。我们把n取半(同时把n加倍),直到n变成奇数。

int multiply3(int n, int a) {
    while (!odd(n)) {
        a=a+a;
        n=half(n);
    }
    if (n==1) return a;
    return mult_acc4(a, n-1, a);
}

但这样的话,mul_acc4就没有必要再检查n=1的情况了,因为n是奇数,n-1是偶数。因此,我们可以在调用的时候直接把n-1取半,把a再加倍,这样就是我们的最终版:

int multiply4(int n, int a) {
    while (!odd(n)) {
        a=a+a;
        n=half(n);
    }
    if (n==1) return n;
    return mult_acc4(a, half(n-1), a+a);
}

重写代码

在不断变换我们的乘法算法的过程中,可以看到重写代码是很重要的。没有谁一上来就能写出很好的代码,很多时候需要反复尝试才能找到最高效或者最通用的方法。程序员不应该有一遍就过的想法。

在重写过程中你可能会想,一条指令的多少不见得能造成多大的区别。但是也可能你的代码在后来很长时间内被重用很多次(实际上,临时处理措施经常变成代码里面存活时间最长的部分),而且,本来可以省掉的操作现在不起眼,留下来也可能会被修改成非常开销非常高的代码。

努力追求效率的另外一个好处是,在此过程中你不得不逼迫自己深入理解问题本身。与此同时,理解深入必然带来更有效率的实现,这是一个良性循环。

2.3 本章的思考

初等代数的学生都会学习如何通过变换表达式的形式来简化它们。在我们实现埃及乘法的不同实现中,我们也经历了类似的变换过程,重新安排代码语句,使之更清晰更高效。每个程序员都应该养成尝试对代码进行变换的习惯,直到得到最佳形式。

本章我们看到了数学在古埃及是怎样出现的,如何留下了第一个有记录的算法。我们要到靠后面的内容才会回到这个算法并对其进行推广。接下来,我们要向后跳跃1000年,领略一些古希腊人的数学发现。

通过截获gettext调用来修改WordPress的部分翻译

这是上一篇文章的中文版。

WordPress网站的页面底端通常有一句话,英文原文是“Proudly Powered by WordPress”,链接到WordPress网站。这句话当前的官方翻译是,“自豪地采用WordPress”。我个人认为,这个翻译是蹩脚的字面翻译,中国人不大会用这种语言致谢。更好的翻译应该是“低调地使用WordPress”。

但是要修改这一句话还不太容易,比较了多种方法之后,我采用了如下方法:

  1. 创建一个child theme。
  2. 在child theme的functions.php中加入如下代码:

[code language=”php”]
function change_attribute_line( $translated_text, $text, $domain ) {
        switch ( $text ) {
            case ‘Proudly powered by %s’ :
               $translated_text = __( ‘低调地使用%s’);
           break;
}
return $translated_text;
}
add_filter(‘gettext’, ‘change_attribute_line’,20,3);
[/code]
就可以了。

代码基本没有需要解释的地方。我在此处搜索$test而不是$translated_text仅仅是因为我觉得匹配英文字符串可能更快一些。

值得注意的是,此处的字符串中仍然包含有占位符”%s”。这就是说,gettext和它的filter是在占位符”%s”被替换之前调用的。

如果我搜索替换后的字符串“Proudly powered by WordPress”,是搜索不到的。我刚开始不明白(不仅WordPress文档没有指出这一点,网上好像也没有人遇到这个问题)所以在此花了很长时间。希望看到帖子的同学们可以节省一些时间。

Replacing WordPress translation by hooking to gettext

In a typical WordPress website, you see this line in the page footer:

“Proudly powered by WordPress”

This line links to WordPress website. The official Chinese translation of this line is:

“自豪地采用WordPress”

I’m using WordPress in several of my websites and I found it’s annoying: It’s a word to word translation and it’s simply not what we Chinese would say when we want to attribute something to someone. After some thinking, I think the appropriate translation should be:

低调地使用WordPress

“采用”emphasis the choosing of WordPress, “使用”emphasis the fact we are using WordPress right now.

The message is clear and typical Chinese: I’m a humble webmaster. All the glories go to WordPress.

So I set out to change the translation. It turned out to be rather complicated, but here’s how I finally accomplished it:

  1. Create a child theme of your current theme, because you’ll have to rewrite its logic and you don’t want your effort to be overwritten by an upgrade of that theme. The procedure of creating a child theme can be found here.
  2. In your function.php file, add the following lines:

[code language=”php”]
function change_attribute_line( $translated_text, $text, $domain ) {
        switch ( $text ) {
            case ‘Proudly powered by %s’ :
               $translated_text = __( ‘低调地采用%s’);
           break;
}
return $translated_text;
}
add_filter(‘gettext’, ‘change_attribute_line’,20,3);
[/code]

That’s it.

Notice that I search for $text instead of $translated_text. You can also search for $translated_text. That would be:
[code language=”php”]
        switch ( $translated_text ) {
            case ‘自豪地采用%s’ :
               $translated_text = __( ‘低调地使用%s’);
           break;
}
return $translated_text;
[/code]
I search for $text because I assume (I’m not sure) search for non-unicode string will be faster.

It’s worth noting that both strings have placeholders in it. Initially I was trying to search the whole string, “Proudly powered by WordPress” or “自豪地采用WordPress” and failed. The fact that the above code works tells us, gettext() and its filters are called before the placeholders get replaced. I spent a lot of time figuring this out. It’s not mentioned in the document, nor did anyone post this on the web. This is actually the key reason why I write this post. Someone from WordPress should update the document.

RSA illustration with not-so-small numbers – part 2

Let’s have a closer look at the encryption. During the communication, what’s been exposed are:

Alice’s public key (n=2627, e=13) , and the encrypted message.

For anyone who’s entered the world of modern cryptography from the old age, it’s tempting to try to decrypt the encrypted message using the encrypting key, the public key.

For these people, I have the below chart that shows the mapping between the plain text and the encrypted data:

encryption_mapping

x-axis is the plain-text data (sorted from 1 to 2627) and y-axis is the encrypted data(from 0 to 2626). I did the calculation using this line of script:

~$ for i in `seq 1 2627`; do echo "$i^13 %2627" | bc; done > /tmp/encryption.mapping

Below is part of this chart zoomed-in:

encryption_mapping_part

So you know the encrypted data, let’s say 2144, and you know the public key (n=2627, e=13). How do you find the number x such that x^13 % 2627 = 2144.

You cannot unless you compute everyone possible 1<x<2627 and then find the correct one. That’s brutal force. This is one of the basic assumption behind the security of RSA: There’s no efficient way to find x. This is called the discrete logarithm problem.

In real world scenarios, the 2 prime numbers will be so large that brutal force is simple impractical.

Then to decrypt the message, one would need the private key. The private key is the modular inverse of phi(n). However, in order to get phi(n), he has to know the factors that form n. And factoring large number is mathematically hard. That is the other assumption behind the security of RSA: There’s no efficient way to factor a large number.

As you will see in other places, these 2 assumptions are the corner stones of modern cryptography.

RSA illustration with not-so-small numbers

Modern cryptography is difficult to understand without illustrations. One of the reason is, modern cryptography involves very large numbers that easily exceed the capacity of a standard calculator, let alone human comprehension. There are some illustrations out there using small numbers. The problem is, the numbers are too small to be convincing. So I’d like to try some no-so-small numbers here. Most of the necessary calculations can be done with GNU bc, so you can try yourself on just any GNU Linux distribution.

Let’s say Bob wants to send the below number to Alice (and make sure only Alice can decrypt the message):

520

Here’s what Alice will do first:

  1. Pick up two distinct prime numbers. The numbers should be sufficiently large so that brutal force is difficult. Here we choose p=37 and q=71.
  2. Calculating n=pq=37*71=2627.
  3. Calculating the n‘s totient function: phi(n)=(p-1)*(q-1)=2520.
  4. Pick a number e between 1 and phi(n) that is co-prime with phi(n). Here we choose 13.
  5. Find number d so that e*d mod (phi(n)) =1. Here we choose 1357. This step cannot be done with bc. Intead, you can try this online calculator. Just put “modinv(13,2520)” in the text field and then press “go” you’ll get the result.

Now Alice has a public key (n=2627, e=13) and a private key (n=2627, d=1357). She can simply distribute her public key to everyone, including Bob.

Now for Bob to encrypt the message 520 to Alice, he has to encrypt the message using Alice’s public key:

520^13 % 2627 = 2235

Now Alice received this number 2235 from Bob. In order to decrypt this message, she do the following calculation(using her private key):

2235^1357 % 2627 = 520

Actually, here Bob can encrypt just any number that is less than or equal to n=2627 in this way.

Bob:

1^13 % 2627 = 1

Alice:

1^1357 % 2627 = 1

Bob:

2^13 % 2627 = 311

Alice:

311^1357 % 2627 = 2

Bob:

3^13 % 2627 = 2361

Alice:

2361^1357 % 2627 = 3

Bob:

4^13 % 2627 = 2149

Alice:

2149^1357 % 2627=4

Bob:

137^13 % 2627 = 2431

Alice:

2431^1357 % 2627 = 137

If his message is large, then he has to split his message into chunks that are smaller than n and encrypt them one by one.

Note that this only illustrates how Bob can send secrete messages to Alice. If Alice wants to send secrete messages to Bob then she has to have Bob do the same first:

  1. Pick up 2 sufficiently large prime numbers;
  2. Get the product of these 2 prime numbers – This is part of the keys;
  3. Get the totient of this product;
  4. Pick a number that is co-prime with this totient but smaller – This combined with the product is the public key;
  5. Find the number that is the multiply modular inverse of this number – This combined with the product is the private key;

Then Bob sends his public key to Alice and Alice can encrypt the messages using Bob’s public key. Upon receiving the messages, Bob can decrypt the messages using his private key.

An ugly implementation

Recently, I’ve been working on an CPropertySheet based application. It has multiple property page. Along with every property page, there’s a working thread behind the scene. I use a timer and a critical section kernel object to syncronize the working thread and the GUI update.

I tested the code with just one page, it works fine. But when I add the second page to the property sheet, the code ends up with an assert error on the SetTimer function call with the OnInitDialog funciton of the PropertySheet.

I was confused by this odd – it takes me a whole day to figure out what happened to the second property page. Finally I realized that the page would not be created till the user clicked the tab and activate the hiding page. So in OnInitDialog function, all propertypages’ m_hwnd is NULL – except for the first page, it has been activated automatically – that lead to the ASSERT error.

To solve this problem, I use a loop to activate the pages within the OnInitDialog function of propertysheet and then SetTimer. It works, but what an ugly implementation! If anyone has a neat approach, pls let me know.