Category Archives: TTTech

T_T Tech…

[Programming Pearls] Column 1 – uninitialized space

Source:
Problem 1.9

Problem Description:
Reserve a set of memory without initialization, meaning the data in there is random. Name the memory as ‘data’.
Design a mechanism to track and tell whether given index in data is initialized or not.

Algorithm 1:
1. Use a set data structure (naming it to initialized_set) to store initialized indices. Initially it is empty.
2. Write value to data[i]: need to initialized_set.add(i)
3. Access value at data[i]: need to check initialized_set.has(i)
Flaw: ops on initialized_set are O(logN).

Algorithm 2:
Purpose: to reduce op cost on initialized_set from O(logN) to O(N).
1. Use array initialized_arr to store initialized indices. Initially it has length of len(data).
2. Use array track_positions: track_positions[i] stores index in initialized_arr, i.e. initialized_arr[track_positions[i]] should equal to i if data[i] is initialized. Initially it has length of len(data).
3. Use a int variable ‘track’ to indicate (assert) initialized_set[0:track] is valid set.
4. Write value to data[i]: need to set initialized_arr[track] = i; track_positions[i] = track; track += 1.
5. Access value at data[i]: need to check initialized_arr[track_positions[i]] == i && track > track_positions[i]

[Programming Pearls] Column 1 – bitvec

持续到2月,坚持至少一周2更, on algorithms…

Use case of a bitvec:
v = bitvec(capacity)
v.set(45322)
v.set(1000234)

if v.has(453222): print(‘ok it’s there’)

Data structure:
bytearray

Key concepts on storage:
- 1M bytes == 8M bits, meaning bitvec on 1MB memory will be able to handle 8M numbers, with the max number to be 8M – 1 (assuming we need to store 0 in the bitvec)
- 10M numbers requests 10/8M bytes – Assuming ‘M’ is always 2^20, not necessarily 1,000,000.

Key implementation problems:
In the case of v.set(num),
- What is the corresponding position (naming it as pos) for num in the underlying bytearray?
pos = num / 8 (or num >> 3)
- What is the corresponding mask index (naming it as r, remainder) in bytearray[pos]?
r = num % 8 (or num & 0×7)

The implementation of bit vector:

class prog_1_2_bitvec:
    """
    >>> v = prog_1_2_bitvec(6)
    >>> v.set(2)
    >>> v.set(5)
    >>> print('%x' % v.bytes[0]) # should be 100100
    24
    >>> v.has(2)
    True
    >>> v.has(5)
    True
    >>> v.unset(2)
    >>> v.has(2)
    False
    >>> v.unset(5)
    >>> v.has(5)
    False
    >>> v.set(3)
    >>> v.has(3)
    True
    """

    def __init__(self, max):
        if max <= 0: raise Exception('invalid max - length <= 0')
        self.bytes = bytearray(int(max / 8 + 1))
        self.max = max

    def set(self, num):
        self._valid(num)
        pos, mask = self._divide(num)
        self.bytes[pos] |= mask

    def has(self, num):
        self._valid(num)
        pos, mask = self._divide(num)
        return self.bytes[pos] & mask == mask

    def unset(self, num):
        self._valid(num)
        pos, mask = self._divide(num)
        self.bytes[pos] &= ~mask

    def count(self):
        # lie: TODO: best way to count the number of 1
        # will cover it later
        pass

    def _divide(self, num):
        d = num >> 3
        r = num & 0x7
        return (d, 1 << r)

    def _valid(self, num):
        if num < 0: raise Exception('invalid num - negative number')
        if num >= self.max: raise Exception('invalid num - num > %d' % self.max)

The implementation of multi-pass bitvec sort:

def prog_1_3_bitsort_multipass(l, max, k):
    """
    TODO: margin conditions...

    >>> max = 1 << 8
    >>> l = prog_1_4_gen_list(max, max >> 1)
    >>> print(l[0:20])  # doctest:+ELLIPSIS
    [...]
    >>> l = prog_1_3_bitsort_multipass(l, max, 2)
    >>> print(l[0:20])  # doctest:+ELLIPSIS
    [...]
    >>> test_verify_sorted(l)
    True
    """
    def sort_singlepass(l, max, i, k):
        # assumptions:
        # 1. let slice = max / k, the func is to handle numbers from range i * slice ~ (i + 1) * slice
        # 2. The mapping relationship is n => n - i * slice
        # exception cases (needs assertion):
        # 1. max / k has remainder
        slice = int(max / k)
        upper = (i + 1) * slice
        lower = i * slice
        v = prog_1_2_bitvec(slice)
        for n in l:
            if n < lower or n >= upper: continue
            mapped_n = n - lower
            v.set(mapped_n)
        sortedlist = [n + i * slice for n in range(slice) if v.has(n)]
        return sortedlist

    sortedlist = []
    for i in range(0, k):
        sortedlist += sort_singlepass(l, max, i, k)
    return sortedlist

Hints:
- Overall: define use case – UT, method signature, etc
- Write logic: write mainstraem logic with notes (in comments) or assertions for exception / corner cases.
- Organize mind: cases!!

Tests (TODO):
- Areas: UT – like UT for sort_singlepass, basic E2E (positive), negative like to pass numbers less than 0 or exceeding data volume, code coverage (previous test cases may have covered most), performance

Bresenham算法

上回说到, 在看一本书《Windows游戏编程大师技巧》 (Tricks of Windows Game Programming Gurus). 这次继续书里的内容: 直线光栅化的Bresenham算法. 书上讲的比较含糊, 没有讲算法的推导过程, 更没讲算法是怎么想出来的. 所以我们只好自己动手, 丰衣足食…

直线光栅化

直线光栅化是指用像素点来模拟直线. 比如下图中用蓝色的像素点来模拟红色的直线. 图中坐标系是显示器上的坐标系: x轴向右, y轴向下.

bresenham

设deltaX = endX – startX, deltaY = endY – startY. 那么斜率为k = deltaY / deltaX. 我们先考虑简单的情况: 当 0 < k < 1即直线更贴近x轴. 在这种情况下deltaY < deltaX, 所以在光栅化的过程中, 在y轴上描的点比在x轴上描点少. 那么就有一个很直观的光栅化算法:

line_bresenham(startX, startY, endX, endY)
{
    deltaX = endX - startX;
    deltaY = endY - startY;
    k = deltaY / deltaX;
 
    <span style="color: #0000ff">for</span> (x = startX, y = startY; x &lt;= endX; ++x)
    {
        <span style="color: #0000ff">if</span> (满足一定条件)
        {
            ++y;
        }
        drawPixel(x, y);
    }
}

基于斜率 / 距离的两个简单直线光栅化算法

好了,貌似很简单, 就剩一个问题: “满足一定条件”是什么? 可以用斜率判断, 也可以用上图中直线与光栅线交点 (红点) 光栅点 (蓝点) 的距离来判断. 继续用伪代码说话:

<span style="color: #008000">// 算法1: 用斜率判断</span>
<span style="color: #0000ff">void</span> line_bresenham_k(startX, startY, endX, endY)
{
    deltaX = endX - startX;
    deltaY = endY - startY;
    k = deltaY / deltaX;
 
    <span style="color: #0000ff">for</span> (x = startX, y = startY; x &lt;= endX; ++x)
    {
        <span style="color: #0000ff">if</span> (x - startX != 0)
        {
            <span style="color: #008000">// 计算当前斜率</span>
            currentK = (y - startY) / (x - startX);
 
            <span style="color: #008000">// 如果当前斜率 &lt; k, 则增加y坐标</span>
            <span style="color: #0000ff">if</span> (currentK &lt; k)
            {
                ++y
            }
        }
        drawPixel(x, y);
    }
}
 
<span style="color: #008000">// 算法2: 用距离判断. 计算直线与光栅线交点y坐标我们需要用到</span>
<span style="color: #008000">// 直线方程 y = k (x - startX) + startY</span>
line_bresenham_dist(startX, startY, endX, endY)
{
    deltaX = endX - startX;
    deltaY = endY - startY;
    k = deltaY / deltaX;
 
    <span style="color: #0000ff">for</span> (x = startX, y = startY; x &lt;= endX; ++x)
    {
        <span style="color: #008000">// 计算直线与光栅线交点的y坐标, 以及与光栅点的距离</span>
        ptY = k * (x - startX) + startY;
        dist = ptY - y;
 
        <span style="color: #008000">// 如果距离 &gt; 0.5或者 &lt; -0.5, 说明我们需要增加y以</span>
        <span style="color: #008000">// 将距离的绝对值控制在0.5之类</span>
        <span style="color: #0000ff">if</span> (dist &gt; 0.5 || dist &lt; -0.5)
        {
            ++y;
        }
        drawPixel(x, y);
    }
}

消灭浮点数!

以上都是很直观的算法, 下面不直观的来了 – 上面的算法都需要在循环体内执行乘法, 准确的说, 是进行浮点数的乘法. 我们怎么能减少这些浮点数的乘法开销呢? 以基于距离的算法2为例: 首先, k是一个浮点数, 0.5也是浮点数. 我们可以通过将这些表达式都乘以2 * deltaX (整数) 来解决浮点数的问题. 伪代码:

<span style="color: #008000">// 算法3: 在算法2的基础上消灭浮点数!</span>
line_bresenham_dist(startX, startY, endX, endY)
{
    deltaX = endX - startX;
    deltaY = endY - startY;
 
    <span style="color: #0000ff">for</span> (x = startX, y = startY; x &lt;= endX; ++x)
    {
        <span style="color: #008000">// 计算直线与光栅线交点的y坐标, 以及与光栅点的距离</span>
        ptY1 = deltaY * (x - startX) + startY * deltaX;
        dist1 = ptY1 - y * deltaX;
        dist1 = dist1 &lt;&lt; 1; <span style="color: #008000">// dist1 = dist1 * 2</span>
 
        <span style="color: #008000">// 如果距离 &gt; 0.5或者 &lt; -0.5, 说明我们需要增加y以</span>
        <span style="color: #008000">// 将距离的绝对值控制在0.5之类</span>
        <span style="color: #0000ff">if</span> (dist1 &gt; deltaX || dist &lt; -deltaX)
        {
            ++y;
        }
        drawPixel(x, y);
    }
}

消灭乘法!

圆满解决浮点数运算问题! 不过…乘法运算还在. 消灭乘法问题的办法比较不直观, 让我们想一想: 还有什么办法能简化运算. 直线方程已经不能再简化, 所以唯一的突破口就是能不能利用递推 / 用上一次循环的计算结果推导下一次循环的计算结果.

首先我们来看看在算法2的基础上 (因为算法2计算红点蓝点之间的距离, 比较直观), 怎么通过第n – 1次循环计算出的dist值 (设为d1) 来推导出第n次循环的dist值 (设为d2). 先回顾一下: dist = 直线与光栅线交点的y坐标 – 相应光栅点的y坐标. 我们从几何上直观地考虑: 在第n次循环中, 我们先根据上一次循环所计算出来的d1, 暂时令d2 = d1 + k, 因为我们要保证-0.5 < d2 < 0.5, 而d1 + k满足d1 + k > –0.5, 所以我们只需要考虑当d1 + k > 0.5时, 我们需要将光栅点y坐标增加1, 并且将d2减去1. 显然, 设y1是第n – 1次循环中光栅点的y坐标, y2是第n次循环中光栅点的y坐标. 我们有
1) d2 = d1 + k –  (y2 – y1)
2) 当d1 + k > 0.5时y2 = y1 + 1, 否则y2 = y1
我们已经能根据上面的两个关系式写出算法, 不过为了消除乘法和浮点数, 我们将这两个关系式两端同时乘以2 * deltaX, 并且设e = 2 * deltaX * d, 则我们有
3) e2 =  e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)
4) 当e1 + 2 * deltaY > deltaX时y2 = y1 + 1, 否则y2 = y1
终于, 没有了乘法 (2 * deltaY在循环体外计算且被简化为左移一位的运算), 没有了浮点数, 根据关系式3) 和 4), 写出算法:

<span style="color: #008000">// 算法4: 在算法2, 3的基础上利用递推消灭乘法和浮点数!</span>
line_bresenham(startX, startY, endX, endY)
{
    deltaX = endX - startX;
    deltaY = endY - startY;
    e = 0;
    deltaX2 = deltaX &lt;&lt; 1;
    deltaY2 = deltaY &lt;&lt; 1;
 
    drawPixel(startX, startY);
 
    <span style="color: #0000ff">for</span> (x = startX + 1, y = startY; x &lt;= endX; ++x)
    {
        <span style="color: #008000">// 关系式3) e2 =  e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)</span>
        <span style="color: #008000">// 关系式4) 当e2 + 2 * deltaY &gt; deltaX时y2 = y1 + 1, 否则y2 = y1</span>
        e += deltaY2;
        <span style="color: #0000ff">if</span> (e &gt; deltaX)
        {
            e -= deltaX2;
            ++y;
        }
        drawPixel(x, y);
    }
}

消灭浮点数! 代数推导

上面递推关系的推导过程是从图形上”直观”地分析得来的, 但是不严密. 我们能不能形式化地证明关系式1), 2), 3), 4)呢? 因为关系式3), 4)和1), 2)能互相推导, 我们只证明3), 4)如下:

在算法3的基础上设第n – 1次循环计算出的dist1值为e1, 对应的y值为y1, 第n次循环计算出的dist1值为e2, 对应的y值为y2. 根据算法3,
dist1 = 2 * deltaY * (x – startX) + 2 * startY * deltaX – 2 * y * deltaX, 则
e2 – e1
= 2 * deltaY * (x – startX) + 2 * startY * deltaX – 2 * y2 * deltaX – [2 * deltaY * (x – 1 – startX) + 2 * startY * deltaX – 2 * y1 * deltaX ]
=  – 2 * y2 * deltaX + 2 * deltaY + 2 * y1 * deltaX
= 2 * deltaY – 2 * deltaX * (y2 – y1)
所以e2 = e1 + deltaY – deltaX * (y2 – y1). 所以我们有关系式
1) e2 = e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)
2) –deltaX <  e1 < deltaX
3) –deltaX < e2 < deltaX
4)  y2 – y1 = 0 或者 1 
我们根据e1 + 2 * deltaY的取值范围进行讨论. 首先, 因为不等式6), 我们有
2 * deltaY – deltaX < e1 + 2 * deltaY < 2 * deltaY + deltaX

情况1: 如果2 * deltaY – deltaX < e1 + 2 * deltaY < deltaX, 则
2 * deltaY – deltaX – 2 * deltaX * (y2 – y1) < e2 < deltaX– 2 * deltaX * (y2 – y1) 
证: 若y2 – y1 = 1, 则 2 * deltaY – deltaX – 2 * deltaX < e2 < deltaX – 2 * deltaX = -deltaX, 所以y2 – y1 = 1不成立. 即情况1中y2 = y1.

情况2:  如果 deltaX < e1 + 2 * deltaY < 2 * deltaY + deltaX, 则
deltaX – 2 * deltaX * (y2 – y1) < e2 < 2 * deltaY + deltaX – 2 * deltaX * (y2 – y1)
反证: 若y2 – y1 = 0, 则 deltaX < e2 < 2 * deltaY + deltaX 所以y2 – y1 = 0不成立. 即情况2中y2 = y1 + 1.

打了这么多字, 累…以上就是当0 < k < 1的情况, 剩余几种情况 (k > 1, –1 < k < 0, k < –1. 不要挑剔我不用”>=”这种符号…) 都可以通过简单的x, y交换以及正负交换来搞定.