Archive for Technology

Using LAPACK in C

简单的说就是如何在C程序中调用Fortran子程序。我不想再学一门语言,所以直接调用库函数是最好的办法。google了一下,在VS2005里调试成功。简单来说你只需要:

  • 编译好的LAPACK库
  • 相应的头文件
  • 知道Fortran数组内存布局的特点
  • 知道如何在C程序中调用DLL
  • 就没问题。dll和头文件可以在这里下载。建议动手前查看一下头文件中的f2c.h,里面定义了一些常用的数据类型,混个眼熟。

    怎么在C程序中链接外部DLL库,这个和编译器有关,就不在这儿罗嗦了。

    主要的问题是怎样在C的数据类型和Fortran中的类型间进行转换。基本的类型,比如REAL, DOUBLEREAL之类的定义,都在f2c.h里写好了。最要紧的一点:Fortran中二维数组是按列存储的,而C是按行存储。当然,我个人的意见是干脆不要用C提供的多维数组,直接用一维数组,元素按Fortran惯例存放就行了,顺便定义一套宏什么的,用起来省事,看起来顺眼。此外,Fortran中函数传参是call by reference,和C的call by value截然相反,所以所有参数都必须是指针。记住这两条就差不多了。

    我要用LAPACK对角化Hamiltonian,感觉应该用zheev, compute all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix

    附精美范例一则,在VS 2005下编译通过:

    #include <stdio.h>
    #include "f2c.h"
    #include "zgeev.h"
    
    #define SIZE 3
    
    int main(int argc, char *argv[]) {
      doublecomplex A[SIZE][SIZE], b[SIZE];
      doublecomplex DUMMY[1][1],WORK[2*SIZE];
      doublecomplex AT[SIZE*SIZE];
      int i, j, ok, c1, c2, c3;
      char c4;
    
      /* initializing A */
    
      for (i = 0; i < SIZE; i++) {
        for (j = 0; j < SIZE; j++) {
          AT[j + SIZE*i].r = A[j][i].r;
          AT[j + SIZE*i].i = A[j][i].i;
        }
      }
    
      c1 = SIZE;
      c2 = 2*SIZE;
      c3 = 1;
      c4 = ‘N’;
    
      zgeev_(&c4, &c4, &c1, AT, &c1, b,
            DUMMY, &c3, DUMMY, &c3, WORK, &c2, WORK, &ok);
    
      if (ok == 0)
        for (i = 0; i < SIZE; i++)
          printf("%f + %fi\n", b[i].r, b[i].i);
      else
        printf("An error occured");
    
      return 0;
    }
    

    并非原创,来自Using LAPACK Routines in C Programs,略微做了些调整。

    阅读(351 次)

    Comments

    暴力求解24点

    是真正的暴力求解……穷举法。不过这个问题似乎也没有更好的方法了,各种算法除了暴力的程度,似乎也没有其他本质的区别。

    好久没写程序,今天写得这个75%是从网上参考过来的。那个产生所有可能的后缀表达式的算法很漂亮,产生所有排列的函数也很简洁。唯一的缺点是,这个方法会搜索大量重复的表达式,总的搜索空间是全排列数乘上Catalan数,比指数增长还要彪悍。但是对一般4个数的24点,这个低效率还体现不出来。此外,输出的表达式是后缀的,看起来不怎么顺眼。但是后缀转中缀又很麻烦,我就懒得再折腾了。

    import sys
    
    def add(a,b):
      return a + b
    def sub(a,b):
      return a - b
    def mul(a,b):
      return a*b
    def div(a,b):
      return float(a)/b
    
    ops = {add:’+‘, sub:’-‘, mul:’*‘, div:’/‘}
    
    def permute(seq):
      seqn = [[seq.pop()]]
      while seq:
        newseq = []
        new = seq.pop()
        for i in range(len(seqn)):
          item = seqn[i]
          for j in range(len(item) + 1):
            newseq.append(item[:j] + [new] + item[j:])
        seqn = newseq
      return seqn
    
    def _exprs(expr, stk_depth, vals, ops):
      if not vals and stk_depth == 1:
        yield expr
      if stk_depth > 1:
        for op in ops:
          for e in _exprs(expr + [op], stk_depth - 1, vals, ops):
            yield e
      if vals:
        for e in _exprs(expr + [vals[0]], stk_depth + 1, vals[1:], ops):
          yield e
    
    def exprs(vals, ops):
      return _exprs([], 0, vals, ops)
    
    def eval(expr):
      stack = []
      for e in expr:
        if callable(e):
          b = stack.pop()
          a = stack.pop()
          try:
            stack.append(e(a, b))
          except ArithmeticError:
            return None
        else:
          stack.append(e)
      return stack[0]
    
    def tostr(expr):
      ret = []
      for e in expr:
        if callable(e):
          ret.append(ops[e])
        else:
          ret.append(str(e))
        ret.append(’ ‘)
      return ”.join(ret)
    
    if __name__ == "__main__":
      vals = [int(x) for x in sys.argv[1]]
      goal = int(sys.argv[2])
      for p in permute(vals):
        for e in exprs(p, ops.keys()):
          if eval(e) == goal:
            print tostr(e)
    

    阅读(359 次)

    Comments

    龙书第二版下载

    compiler-aho.pdf
    今天终于把龙书第二版的电子书传上去了。不知道链接什么时候会失效。
    失误,一开始写成第三版了

    阅读(722 次)

    Comments (3)

    Evil of Memory Leak

    本篇含技术内容,不喜者勿入。

    这是一个关于内存泄漏的凄惨故事。内存泄漏,乃是C/C++ Programmer最恐惧的Bug之一;来不知其来,去往往也不知其去,潜伏在层层花括号深处,当你发现它的时候,往往已经是程序崩溃之时。我几天前还被一个诡异的memory bug困扰了一晚,此bug举止异于常bug,在调试时单步执行没有任何问题,但一运行就出一些稀奇古怪的结果。以至于我都无法定位bug的位置,最后不得不借助于石器时代的printf调试大法,在可能出错的位置加上一摞一摞的printf才找到bug的可能位置。就算我已经定位到某一行代码,在我盯着屏幕看了半天,看得眼前快要出现幻觉,还是想不通bug从何而来。果然,从哪里来,到哪里去,是个终极关怀级别的问题。这一块代码和内存分配有关,无奈之下只好重写一遍,用了个稍微不同的逻辑。然后bug消失,程序再次和谐。

    我个人喜欢用C写东西,尽管已经有无数与memory leak搏斗的经验,每次仍然免不了犯这类错误,经验虽然不停增长,错误好像未见减少。所以有时候觉得Java, C#这样的带Garbage Collection的语言也不错,但昨天看了一篇文章,讲述了一个C#中的内存泄漏怎样导致200万美金旁落的悲惨故事。涉及到这么大的支票,我们的兴趣总会提高一些。

    文章作者是Princeton大学电子工程系的学生,他(当然还有一堆牛人)参加了一项汽车自动驾驶比赛,奖金是2个million。他们控制汽车的程序是用C#写的,有1万行代码。在比赛的最后关头,汽车老是开一会就失去响应,然后要么一头撞到障碍物(例如灌木丛),或者就茫然地在大地上行走直到没油。分析起来,问题出在代码的障碍物探测模块,它负责寻找汽车前进道路上的障碍。当然,如果车已经开过去了,那么刚才路上的障碍信息就没用了,所以会被垃圾回收器处理掉。最后的Bug是,这些对象,理应被丢弃掉的废物,根本没有被回收。一开始只见内存占用像滚雪球一样越来越大,开了一会汽车就歇了。可想而知,2百万肯定是黄了,最郁闷的是都这样了还不知道为什么垃圾回收器就是不干活,正宗的死不瞑目。

    比赛完了,牛人们坐下来开会总结经验教训。其中某牛从网上下了个ANTS Profier,分析了他们代码的内存使用,终于看到了问题所在。前面说过这是一个C#式的内存泄漏,并不是说垃圾回收器出了问题,没有处理该回收的内存。虽然微软很不招人喜欢,但要相信它还不至于犯这种错误。垃圾回收器判断一个对象是否该被回收,是根据它能否被当前代码所引用。如果代码里已经不能访问到这个对象,那么它自然是垃圾。出问题的程序里,虽然程序员已经清除了到这些废物对象的引用,但这些对象作为事件的订阅者(subscriber,基本上是C#程序员的黑话),还在被其他活跃对象默默地引用。最后的结果自然是内存耗尽。所以那篇回顾文章的作者在标题里就痛心疾首地喊:”If Only We’d Used ANTS Profiler Earlier…”,言下之意,还在为那几乎到手的200万美金心疼。

    这个故事告诉我们一个道理:垃圾回收机制虽然很方便,但它实际上是把一个问题变成了另外一个问题(手动释放内存到去掉所有引用)。相对来说,发生错误的可能性减少了,“去掉所有引用”这句话是我在这里说说的,听起来有些别扭,实际上在代码中是一件比较自然的事,不需要特别操心。但是一旦出错,也更加难以察觉——因为我们几乎不可能跟踪垃圾收集器的动作。所谓一饮一啄,天下没有免费的午餐。

    阅读(641 次)

    Comments (2)

    龙书第二版

    说到龙书(Dragon Book),喜欢玩编译器和程序语言的老大们肯定是人手一册,在编译器技术这一块的教科书中地位几乎有如圣经。说到这儿又要鄙视一下国产的什么编译原理教材,简直是不知所云,抄都抄得这么潦草。

    计算机领域的经典著作,常常会有些“小名”或俗名,如无人不知的TAOCP,C语言圣经K&R,算法入门宝典CLR。TAOCP是书名的缩写,此书堪称最有名的算法书籍,我估计也是被完整阅读次数最少的计算机书籍,之一都省掉了。K&R和CLR都是作者名字的开头第一个字母。龙书算是比较有个性的,这个外号来自其封面。1977年,Alfred V. Aho和Jeffrey D. Ullman写了一本Principles of Compiler Design,封面是一位骑士枪挑一头绿色恐龙,因此得名龙书。到了1986年,此书升级换代,书名变成了Compilers: Principles, Techniques and Tools,沿用至今,作者增加了Ravi Sethi。封面依然是一位骑士和一头恐龙,只不过恐龙变成了红色,骑士则坐在了电脑前。仔细看的话,发现龙身上纹身“Complexity of Compiler Design”,骑士全副武装,盔甲上书“Data Flow Analysis”,夹着一把写着”LALR Parser Generator”的宝剑,盾牌上则写着”Syntax Directed Translation”。

    一晃过了20年,龙书终于出了新版。封面当然是骑士斗恶龙。作者又加了一位,篇幅猛增了200多页,变化十分巨大。除了一些古典内容如词法分析,语法分析,语法制导翻译,中间代码生成没怎么动,后面的运行环境,代码优化部分可谓焕然一新。从Run-time Environment这一部分就可以看出程序语言这些年来的变化:Java,.NET,Python,Ruby…自动内存管理已经不是什么新鲜事,所以龙书里开了近半章的篇幅给垃圾收集。内存管理是个大题目,完全铺开来一本书都不够,从目录上看,龙书已经包括了现在主流的一些技术,估计作者安排这一块内容时也死了不少脑细胞。变化最大的当属代码优化,原来只有一章,现在一下子撑成了四章,洋洋洒洒近400页,可以单独成一本书了,较之第一版,作者似乎不满足于做一本入门教材,要一揽子把从入门到高级全包了。

    其他内容上的变化,暂时只能从目录上看一下。Flow Analysis增加了一些较新的算法,比如Region-Based Analysis。优化部分花了不少力气讲Parallel Machine,现在看来也是大势所趋。

    书肯定是不可能仔细看了;也不知道什么时候会有中文版;这么一本1000页的大砖头翻成中文后也不知道厚成什么样。我个人的一点观察,最近几年动态语言与函数式语言开始逐渐成为主流,用马克思的话,一个函数式编程的幽灵开始在C++/C#/Java中徘徊。现在的编译技术,似乎应该包括Dynamic Typing, JIT等等,特别是种种针对动态语言的优化(比如方法缓存)。作为一本教材,不可能面面俱到,更不能跟风,编者总有自己的原则和取舍,大概龙书的四位作者认为这些东西还不值得进入大学的编译技术课程。但我觉得现在正是一个计算环境变革的前夜;或者说变革已经在进行时了。这样的话,也许下一个版本,不用再让我们等20年。

    说到这儿,又要损一下中国这帮写教科书的人。为了避免打击面过广,不失一般性,限定于计算机行业。老外写书,一版二版三版,真是与时俱进。当年Dirichlet写了一本数论讲义,每次重版他的学生Dedekind都会在后面添个附录反映一下最新的研究进展,到最后附录比正文还厚得多。在我国,就少见这种延续性。人家真是把教材当作一项事业来做的。中国最有名的计算机书估计是谭老师的C语言,去年貌似也出了个新版,除了换个封面,里面新瓶装旧酒,继续误人子弟。对谭老师C语言书的声讨网上已经进行多年,我也没有重复的必要。只是想到这些状况,不能不有点觉得外国月亮有点圆。

    本来想把电子版放上来,但尺寸实在吓人(近50M),我估计赛族上对这个感兴趣的人也不会太多,如果想要电子版的话,私人联系吧。当然直接google也可以。

    下面是三版龙书封面展览:
    db_2.jpgdb_3.jpg

    db_1.jpg

    阅读(1172 次)

    Comments (9)

    Youtube Video Downloader in Ruby

    #Written by Cheng Meng
    #Created on 2007-7-27
    
    requireopen-uri‘
    
    url_video = ""
    youtube = "http://www.youtube.com/"
    ARGV[0] =~ /(?:#{youtube}(?:watch\\?v=|v\\/))?(.*)$/
    video_id = $1
    open("#{youtube}watch\\?v=#{video_id}") do |f|
      f.each_line do |line|
        if line =~ /SWFObject\\("\\/player2.swf\\?([^"]+)"/
          url_video << "#{youtube}get_video?#{$1}"
          break
        end
      end
    end
    #puts "video:#{url_video}"
    video = File.new("video_#{video_id}.flv","wb")
    open(url_video) do |f|
      f.each_byte { |c| video.putc c }
    end
    

    This program is based on the idea explained in Peteris Krumins’s article Downloading Youtube videos with gawk.But I don’t know any gawk(巨拗口,in Chinese),and I just learn a little Ruby.So I decided to hack my own.

    usage:save the above code as .rb file.This program accepts one argument to specify the url of the youtube video page.Three forms of urls are all permitted:the watch page,such as http://www.youtube.com/watch?v=En0A8KGMgq8,or http://www.youtube.com/v/En0A8KGMgq8,which usually contained in embedded code;even just the ID,En0A8KGMgq8,is OK.Then the script parses the argument,finds the actual url of video,and saves video to a “video_ID.flv” file.

    Note:Linux users may need to change the mode option “wb” to “w”(I’m not sure about this,for I have not a linux system in my laptop).Only DOS/Windows needs this “b” to prevent Ruby from converting ‘\a’ to ‘\d\a’.I waste a lot of time in this problem.

    Update:according to Corsair,the task of downloading should be given to more professional tools,such as wget.Just remove the last 3 lines and add “puts url_video”.Use pipe to connect my script to wget,like this:

    wget -O output_filename `video_downloader.rb url`
    

    video_downloader.rb is the modified script.

    阅读(1109 次)

    Comments (10)

    Creating library with GCC

    A Short Definition

    A static library is basically a collection of object files.During the linking phase of compilation,the static library,or the object files in it, will be simply copied into the target executable file.Static libraries usually ends with a “.a” suffix in UNIX variants(including MinGW) and “.lib” in Windows.

    A shared library,also called dynamic library(on Windows,they are the familiar dlls,the origin of “DLL Hell”) is “dynamically” linked into program in runtime.Why I use “dynamically” here is that during compilation the object files contained in shared library are not inserted into the final executable file;instead, when the program is started, a program in the system (called a dynamic loader) checks out which shared libraries were linked with the program, loads them to memory, and attaches them to the copy of the program in memory.

    Creating Static Library

    Suppose you have some object files with the name util1.o,util2.o,….Then then following command will create a static library “libutil.a” and puts all the objs into it:

    ar rcs libutil.a util1.o util2.o ...
    

    Note that the static library must start with “lib” and have the suffix “.a

    ar can be used to creating,modifying and checking static libraries.
    Read the rest of this entry »

    阅读(662 次)

    Comments (1)

    Tip:Code Layout

    用一个很原始的办法解决了代码显示的问题:把代码放在pre标签里头,就能正常显示空格了。这个技术在其他一些blog上也不很管用,比如blogger,似乎可视化的HTML编辑器会对输入内容动很多手脚,把空格什么的都删除了,解决办法是关掉可视化编辑器,直接写html。所幸在Scinese上pre标签还能正确显示,这样只要把代码在Vim里转成html再贴过来就行了。剩下的问题是字体大小,因为页面宽度比较小,一旦有较长的行就只能显示一部分,变成了残缺美。更可恨的是,代码在firefox下和IE下显示的字体大小居然大不相同!我不是很懂html和css,碰到这种情况只能哀叹一下并坚定以后找时间学css的决心,别无他法了。

    贴一个煞笔然而无处不在的Hello World:

    #include <stdio.h>
    
    int main(int argc,char **argv) {
        printf(“hello,world!n);
        return 0;
    }

    按照corsair的意见,color这种史前文物不符合xhtml标准,应该毫不留情地抛弃……但是wordpress又不让inline css,两难之下,还是觉得高亮效果比符合xhtml标准更重要,继续使用color。

    阅读(879 次)

    Comments (1)

    How GCC Divide

    While I’m working on the module of my project which translating arithmetic operations into machine code,I think I should take a look at how gcc does it,especially (one of)the most wierd operation in i386 instruction set,div.After examining gcc’s output(using -S option when compiling),I find something interesting.When divisor is a constant integer,for example,an expression like x/3 (assuming x is a local variable with type int)in C source code,gcc will translate it into a seemingly strange sequence of instructions:

    movl  $1431655766, %eax
    imull x
    movl  %edx, %ecx
    movl  x, %eax
    sarl  $31, %eax
    subl  %eax, %ecx
    

    (the result is stored in %ecx)

    Read the rest of this entry »

    阅读(847 次)

    Comments (2)

    Age of Abundance

    装逼文字一篇,厚颜地blog出来:

    现在,每天相当于有1000G的数据被传输到互联网上。这相当于5万个视频短片,2.5亿个杂志故事,或者5亿个博客。重点不在于供给,而变成了选择的智慧。

    是的,现在的问题是:太多了。信息不断的被制造出来,然后被遗忘在网络的某个不知名角落,呆在某块磁盘的一条不起眼的磁道上,直到某天一条delete命令扫过,重新归于零。五年前——甚至两年前——虽然网络已经是个庞然大物,但仍然可以掌握,除了几个大的门户网站和论坛,别的什么个人主页之类的,自生自灭,和传统媒体世界在结构上,或者说拓扑上没有太大的区别——我们可以说两者是同构的,同胚的(在这里我卖弄了两个数学词汇)。慢慢的情况有了变化,信息的积累到了一定程度,就需要有更高一层进行处理的工具出现:筛选,过滤,组织。这就是搜索引擎。在我的字典里Search Engine就等于Google,Baidu只是用来找找mp3。好了,就算这样,“太多了“问题仍然没有解决:Google原则上让我们能够找到网络上的一条狗,同时会找出无数条差不多的狗。虽然Google有伟大的PageRank算法,但仍然只是机器,还不能和人肉搜索相比。那么只好人多出一点力,先把信息组织一下。当然,这样的努力早就有了,比如Yahoo,比如新浪。

    但是新时代的特征是“You”,是slashdot,digg,wiki,flickr,youtube,……,是blog,是长尾理论,是三个臭皮匠。让所有人都成为信息生产者(本来他们只是被动的消费),这就把原来星形的信息传播方式——从少数拥有者向多数消费者,变成了一个复杂的连通图。这个想法似曾相识,仿佛就是自由贸易和市场经济,看不见的手。我每天在网上的主要活动之一就是在reddit,slashdot(及它的中国后裔solidot),cnbeta上刷新,在Google Reader上看blog(rss真是伟大的发明)。这些聚合类站点的基础是民主,从某种意义上和人民代表大会一个性质。我认为本质上,它们就是事实上的人肉搜索引擎。当然民主有它的缺点,那就是多数人打压少数人,学究点叫多数人的暴政。必须承认,阳春白雪,曲高和寡是存在的;好东西往往不被大众所认知;大众肯定的不一定就好,等等很多”往往“,“不一定”,大量的例外。但是没有办法,完美是不存在的,我们需要的是一个可行解,并不一定最优。

    我想这个变化可以称得上是革命性的。我们找到了一种办法,来解决“太多了”这个难题。这个办法并不陌生,几百年来已经成为人类社会的基础。我准备了几个神秘的说法,比如:在系统包含足够多的成员时,自组织就可以实现。相互作用导致了复杂结构。等等。总之,我们处在一个信息的暴涨时代,一个age of abundance,一个age of uncertainty。

    p.s. 昨天把自己一些乱七八糟的想法写上来之后,看到一篇专业人士的文章,The Role of Humans in Google Search by Matt Cutts,比我从一个用户的角度泛泛而谈要深入的多。Google黑板报上有中文翻译

    阅读(802 次)

    Comments (3)

    « Previous entries