趣题:理想模型下的排序算法(上)

    当我们研究复杂度时,我们往往会将现实机器进行理想化。例如,我们说冒泡排序是O(n^2)的,这其实是不准确的。这个论断假设整数之间的比较运算是O(1)的,而事实上它们是O(log(min(|a|,|b|)))的。多数时候我们都认为这种机器模型的理想化是合理的,毕竟这让问题简化了不少,并且也能反映出算法的本质。但大家有想过吗,这个“大整数随便算”的假设其实是一个超级大漏洞,我们可以利用理想模型中的这一漏洞来作弊,获得时间复杂度更低的算法。上个月,Michael Brand在他的UyHiP里就提出了这样一个问题:假设计算机对任意大整数的赋值、四则运算、取余运算、比较运算、位运算(包括左移右移)的复杂度都是常数级别,你能否设计出一个O(n)的排序算法来?

    我非常喜欢这个题目。月初的时候我就提交了一个正确的算法。我们将用左移和加法运算把整数序列编码成一个超大整数,然后利用排序网络进行并行排序。这个算法比较复杂,你可以按照下面的思路一步一步得到这个算法。

1. 如何用位运算来取绝对值

2. 给出两个正整数a, b,不用比较运算和判断语句如何把小数赋给a,大数赋给b?
    提示:和加差除以2等于大数,和减差除以2等于小数

3. 如何利用位运算把整数序列编码成一个超大整数?
    例如把(二进制数)11, 1011, 1110, 1编码为一个数00011 01011 01110 00001

4. 如何用位运算给超大整数中的所有数同时取绝对值?

5. 给出两个超大整数a, b,不用比较运算和判断语句如何把对应位置上的小数赋给a的对应位置,大数赋给b的对应位置? 例如把
      a = 000010 000111 000100 001001
      b = 000001 001011 000011 011111
    变成
      a = 000001 000111 000011 001001
      b = 000010 001011 000100 011111

6. 如何实现奇偶移项排序

    最后,由于奇偶移项排序只有O(n)层,因此整个算法是O(n)的。

    但是,这个算法太繁琐了,不具有美观性。虽然这个算法是我自己想出来的,但我仍然很不满意。待我看了这个月Michael Brand发布的答案后,我一拍大腿,哎呀,还有一个如此简单巧妙的算法我没想到!相比之下,我的算法太复杂了,原因就在于我还没有充分挖掘到“大整数的常数级运算”的潜力。这个理想模型的假设太强大了。打开思路,放宽思维,大胆想象,从更大的尺度上来思考,我们可以得到一个简单得出奇的线性排序算法来。

Read more…

趣题:用位运算生成下一个含有k个1的二进制数

    在所有8-bit的整数中,含有k个数字“1”的二进制数一共有C(8,k)个。给出其中的一个二进制数,你如何利用位运算快速找到下一个恰有k个“1”的数?例如,如果给你二进制数01011100,那么下一个(含4个“1”的)数就是01100011。在继续阅读下去之前,建议你仔细思考一下。你或许会想看看我很早以前写的一篇介绍位运算的文章。这是一道很好的题目,很多位运算技巧在这里都有体现。

    在草稿纸上随便举几个例子,规律很容易看出来。由于“1”的个数是固定的,为了让这个二进制数更大,我们必须把第一个出现在“1”左边的“0”改成“1”;同时,为了让这个二进制数尽可能小,我们必须把它右边那些“1”重新排到最低位去。
    更具体地说,下一个二进制数可以通过以下步骤得到:找到右起第一个单个的或连续的数字“1”,把它们全改成“0”,同时把它们左边的那个“0”改为“1”。此时,“1”的个数可能减少了,我们只需把还差的“1”摆在最右边就行了。举个例子,01011100的右起第一个“1”在第三位,把它和左边紧挨着的“1”一并变为“0”,并把再左边那个“0”变为“1”,于是我们得到01100000。我们还差两个“1”,把这两个“1”补在最低位得到01100011即可。现在我们的任务是,想出一个用位运算来实现这些步骤的办法。
    我们已经熟知,用x & -x可以提取最右边的那个“1”。当意识到可以利用加法来消除连续的“1”时,我们很快得到了第一步操作的位运算实现:把x & -x加到x上,利用二进制加法的进位把“..01111..”变成“..10000..”。现在,我们需要计算出刚才的操作中一共“跳过”了多少个“1”,换句话说现在的x的右起第一个“1”和原来的x的右起第一个“1”差了多少位。关键就在这里!我们可以用除法来完成这一步,例如100000除以100就相当于把被除数右移2位,得到的结果即可以表示两个数中的“1”差了多少位。在最低位产生指定数量的“1”需要用到另一个技巧:减1操作可以把右边连续的“0”都变成“1”,即把…10000变成…01111。我们得到了该问题的第一个算法:

b = x & -x;
t = x + b;
c = t & -t;
m = (c/b >> 1) - 1;
r = t | m; //最终结果

    我们对上述算法做一个简单的说明:

操作              | 样例     |  说明
——————+———-+—————————-
x                 | 01011100 |  原数
b = x & -x        | 00000100 |  提取x的右起第一个“1”
t = x + b         | 01100000 |  把x的右起第一个位于某个“1”左边的“0”变成“1”,并把它右边的那些“1”都变为“0”
c = t & -t        | 00100000 |  提取t的右起第一个“1”
c / b             | 00001000 |  右移c中的那个“1”,其结果中最低位连续的“0”的个数正好是c和b中的“1”相差的距离
m = (c/b >> 1) – 1| 00000011 |  在最低位产生数字“1”,其个数比上述的“距离”少1
r = t | m         | 01100011 |  最终结果

    除去赋值,我们一共用了9个运算符。有可能用更少的运算么?

Read more…

旧闻一则:神秘的0x5f3759df 不可思议的Quake III源码

    Quake III公开源码后,有人在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍:
float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 ); // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}

    code/common/cm_trace.c中也出现了这样一段解释sqrt(x)的函数,与上面的代码唯一不同的就是这个函数返回的是number*y:
/*
================
SquareRootFloat
================
*/
float SquareRootFloat(float number) {
    long i;
    float x, y;
    const float f = 1.5F;

    x = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( f - ( x * y * y ) );
    y  = y * ( f - ( x * y * y ) );
    return number * y;
}

    这样的代码速度肯定飞快,我就不用多说了;但算法的原理是什么呢?其实说穿了也不是很神,程序首先猜测了一个接近1/sqrt(number)的值,然后两次使用牛顿迭代法进行迭代。根号a的倒数实际上就是方程1/x^2 – a = 0的一个正实根,它的导数是-2/x^3。运用牛顿迭代公式x' = x – f(x)/f'(x),式子化简为x' = x * (1.5 – 0.5a * x^2)。迭代几次后,x的值将趋于1/sqrt(a)。
    但这段代码真正牛B的是那个神秘的0x5f3759df,因为0x5f3759df – (i >> 1)出人意料地接近根号y的倒数。人们都不知道这个神秘的常数是怎么来的,只能把它当作神来膜拜。这个富有传奇色彩的常数到底咋回事,很少有人说得清楚。这篇论文比较科学地解释了这个常数。

Tupper自我指涉公式:图象里竟然包含式子本身

    你认为,一个函数图象里是否有可能包含这个函数本身的“图象”?难以置信的是,还真有人构造了这样一个东西。2001年,Jeff Tupper发表的一篇论文里提到了这样一个有趣的不等式:
  
    在0 <= x <= 105,n <= y <= n + 16的范围内,这个不等式对应的图象是这个样子:
  

其中,n = 96093937991895888497167296212785275471500433966012930665150551927170280239526642
46896428421743507181212671537827706233559932372808741443078913259639413377234878
57735749823926629715517173716995165232890538221612403238855866184013235585136048
82869333790249145422928866708109618449609170518345406782773155170540538162738096
76025656250169814820834187831638491155902256100036523513703438744618483787372381
98224849863465033159410054974700593138339226497249461751545728366702369745461014
655997933798537483143786841806593422227898388722980000748404719

    你会觉得这个很神奇吗?你也许会想,天哪,这个是怎么构造出来的啊!但仔细思考之后,你会发现这个一点都不神奇。事实上明白了道理之后你可以构造出无数个这样的式子来。现在给你一些时间让你思考一下,你能否看出其中的奥秘?

    就像魔术揭秘一样,说穿了真相后上面的这些东西就一点意思都没有了。在这个式子里,涉及到x和y的变量时都加上了取整符号,因此整个图象都是一格一格的。这样,不等式右边的式子就简化为y div 17 * 2^(-17x – y mod 17) mod 2,其中x和y都为整数。接着观察,一个数乘以2的负k次方相当于对应的二进制数右移k位,那么x * 2^(-k) mod 2实质上就是二进制数x右起第k位上的数字。对于某个自然数t,当17t <= y < 17(t+1)时,指数-17x – y mod 17恰好对应所有的负整数,于是位于y=17t和y=17t+16之间的图象的每个像素和t的二进制中的每一位数字一一对应。随着t值的增加,图形的像素会一点一点地变化。当纵坐标足够大时,必然会出现一段高度为17的图象,图象的样子和不等式本身的样子相同。当然,你也可以在里面“找到”任何你想要的图象,只需要把图象还原为二进制数并转换为十进制即可。你甚至可以告诉你的MM,说你发现了一个函数,函数在某个位置的图象正好是某某某我爱你的字样。

Matrix67原创
转贴请注明出处
最近发现了一些很不厚道的人,希望大家注意哦!

神奇的分形艺术(三):Sierpinski三角形

    在所有的分形图形中,Sierpinski三角形可能是大家最熟悉的了,因为它在OI题目中经常出现,OJ上的题目省选题目中都有它的身影。这篇文章将简单介绍Sierpinski三角形的几个惊人性质。如果你以前就对Sierpinski三角形有一些了解,这篇文章带给你的震撼将更大,因为你会发现Sierpinski三角形竟然还有这些用途。

Sierpinski三角形的构造
      
    和之前介绍的两种图形一样,Sierpinski三角形也是一种分形图形,它是递归地构造的。最常见的构造方法如上图所示:把一个三角形分成四等份,挖掉中间那一份,然后继续对另外三个三角形进行这样的操作,并且无限地递归下去。每一次迭代后整个图形的面积都会减小到原来的3/4,因此最终得到的图形面积显然为0。这也就是说,Sierpinski三角形其实是一条曲线,它的Hausdorff维度介于1和2之间。

    Sierpinski三角形的另一种构造方法如下图所示。把正方形分成四等份,去掉右下角的那一份,并且对另外三个正方形递归地操作下去。挖个几次后把脑袋一歪,你就可以看到一个等腰直角的Sierpinski三角形。

      

    Sierpinski三角形有一个神奇的性质:如果某一个位置上有点(没被挖去),那么它与原三角形顶点的连线上的中点处也有点。这给出另一个诡异的Sierpinski三角形构造方法:给出三角形的三个顶点,然后从其中一个顶点出发,每次随机向任意一个顶点移动1/2的距离(走到与那个顶点的连线的中点上),并在该位置作一个标记;无限次操作后所有的标记就组成了Sierpinski三角形。下面的程序演示了这一过程,程序在fpc 2.0下通过编译。对不起用C语言的兄弟了,我不会C语言的图形操作。
{$ASSERTIONS+}

uses graph,crt;

const
   x1=320;  y1=20;
   x2=90;   y2=420;
   x3=550;  y3=420;
   density=2500;
   timestep=10;

var
   gd,gm,i,r:integer;
   x,y:real;

begin
   gd:=D8bit;
   gm:=m640x480;
   InitGraph(gd,gm,'');
   Assert(graphResult=grOk);

   x:=x1;
   y:=y1;
   for i:=1 to density do
   begin
      r:=random(3);
      if r=0 then
      begin
         x:=(x+x1)/2;
         y:=(y+y1)/2;
      end
      else if r=1 then
      begin
         x:=(x+x2)/2;
         y:=(y+y2)/2;
      end
      else begin
         x:=(x+x3)/2;
         y:=(y+y3)/2;
      end;
      PutPixel(round(x),round(y),white);
      Delay(timestep);
   end;
   CloseGraph;
end.

Sierpinski三角形与杨辉三角
    第一次发现Sierpinski三角形与杨辉三角的关系时,你会发现这玩意儿不是一般的牛。写出8行或者16行的杨辉三角,然后把杨辉三角中的奇数和偶数用不同的颜色区别开来,你会发现杨辉三角模2与Sierpinski三角形是等价的。也就是说,二项式系数(组合数)的奇偶性竟然可以表现为一个分形图形!在感到诧异的同时,冷静下来仔细想想,你会发现这并不难理解。
      
    我们下面说明,如何通过杨辉三角奇偶表的前四行推出后四行来。可以看到杨辉三角的前四行是一个二阶的Sierpinski三角形,它的第四行全是奇数。由于奇数加奇数等于偶数,那么第五行中除了首尾两项为1外其余项都是偶数。而偶数加偶数还是偶数,因此中间那一排连续的偶数不断地两两相加必然得到一个全是偶数项的“倒三角”。同时,第五行首尾的两个1将分别产生两个和杨辉三角前四行一样的二阶Sierpinski三角形。这正好组成了一个三阶的Sierpinski三角形。显然它的最末行仍然均为奇数,那么对于更大规模的杨辉三角,结论将继续成立。

Sierpinski三角形与Hanoi塔
    有没有想过,把Hanoi塔的所有状态画出来,可以转移的状态间连一条线,最后得到的是一个什么样的图形?二阶Hanoi塔反正也只有9个节点,你可以自己试着画一下。不断调整节点的位置后,得到的图形大概就像这个样子:
      
    如果把三阶的Hanoi塔表示成无向图的话,得到的结果就是三阶的Sierpinski三角形。下面的这张图说明了这一点。把二阶Hanoi塔对应的无向图复制两份放在下面,然后在不同的柱子上为每个子图的每个状态添加一个更大的盘子。新的图中原来可以互相转移的状态现在仍然可以转移,同时还出现了三个新的转移关系将三个子图连接在了一起。重新调整一下各个节点的位置,我们可以得到一个三阶的Sierpinski三角形。
  
    显然,对于更大规模的Hanoi塔问题,结论仍然成立。

Sierpinski三角形与位运算
    编程画出Sierpinski三角形比想象中的更简单。下面的两个代码(实质相同,仅语言不同)可以打印出一个Sierpinski三角形来。
const
   n=1 shl 5-1;
var
   i,j:integer;
begin
   for i:=0 to n do
   begin
      for j:=0 to n do
         if i and j = j then write('#')
         else write(' ');
      writeln;
   end;
   readln;
end.

#include <stdio.h>
int main()
{
    const int n=(1<<5)-1;
    int i,j;
    for (i=0; i<=n; i++)
    {
        for (j=0; j<=n; j++)
           printf( (i&j)==j ? "#" : " ");
        printf("n");
    }    
    getchar();
 &n
bsp;  return 0;
}

    上面两个程序是一样的。程序将输出:
#                              
##                              
# #                            
####                            
#   #                          
##  ##                          
# # # #                        
########                        
#       #                      
##      ##                      
# #     # #                    
####    ####                    
#   #   #   #                  
##  ##  ##  ##                  
# # # # # # # #                
################                
#               #              
##              ##              
# #             # #            
####            ####            
#   #           #   #          
##  ##          ##  ##          
# # # #         # # # #        
########        ########        
#       #       #       #      
##      ##      ##      ##      
# #     # #     # #     # #    
####    ####    ####    ####    
#   #   #   #   #   #   #   #  
##  ##  ##  ##  ##  ##  ##  ##  
# # # # # # # # # # # # # # # #
################################

    这个程序告诉我们:在第i行第j列上打一个点当且仅当i and j=j,这样最后得到的图形就是一个Sierpinski三角形。这是为什么呢?其实原因很简单。把i和j写成二进制(添加前导0使它们位数相同),由于j不能大于i,因此只有下面三种情况:
    情况一:
    i = 1?????
    j = 1?????
    问号部分i大于等于j
    i的问号部分记作i',j的问号部分记作j'。此时i and j=j当且仅当i' and j'=j'

    情况二:
    i = 1?????
    j = 0?????
    问号部分i大于等于j
    i的问号部分记作i',j的问号部分记作j'。此时i and j=j当且仅当i' and j'=j'

    情况三:
    i = 1?????
    j = 0?????
    问号部分i小于j
    此时i and j永远不可能等于j。i' < j'意味着i'和j'中首次出现数字不同的那一位上前者为0,后者为1,那么i和j做and运算时这一位的结果是0,与j不等。

    注意到,去掉一个二进制数最高位上的“1”,相当于从这个数中减去不超过它的最大的2的幂。观察每一种情况中i,j和i',j'的实际位置,不难发现这三种情况递归地定义出了整个Sierpinski三角形。
    嘿!发现没有,我通过Sierpinski三角形证明了这个结论:组合数C(N,K)为奇数当且仅当N and K=K。这篇文章很早之前就计划在写了,前几天有人问到这个东西,今天顺便也写进来。
    另外,把i and j=j 换成i or j=n也可以打印出Sierpinski三角形来。i and j=j表示j的二进制中有1的位置上i也有个1,那么此时i or (not j)结果一定全为1(相当于程序中的常量n),因此打印出来的结果与原来的输出正好左右镜像。

Matrix67原创
转贴请注明出处

网友Voldemort在12楼和13楼很辛苦地帖了一个杨辉三角模2问题的扩展,大家可以看看