显示标签为“算法”的博文。显示所有博文
显示标签为“算法”的博文。显示所有博文

2012年9月1日星期六

Knuth-shuffle洗牌算法

要写出一个真正随机(我意思是均匀分布)的洗牌程序并不容易,即使强如Jeff Atwood也是 会写错的。Jeff写错的版本是这样的:

//error!
for(int i=0;i<cards.Length;i++){
    int n = rand.Next(cards.Length);
    Swap(ref cards[i],ref cards[n]);
}

问题在于我们并不能轻易地看出这段程序生成的排列不是均匀分布的(容我想想,我现在看 得出么…)。

实际上,是有个O(n)的洗牌算法 Fisher-Yates,又称为Knuth-shuffle(Knuth在TAOCP介绍过 而得名):

//To shuffle an array a of n elements (indices 0..n-1):
for i from n − 1 downto 1 do
       j <- random integer with 0 ≤ j ≤ i
       exchange a[j] and a[i]

这个算法就地洗牌,但它需要预先知道洗牌数组的长度。它的正确性(也就是生成的排列是 否均匀分布)不难看出。

Knuth-shuffle还有另一种实现:

//To initialize an array a of n elements to a randomly shuffled copy of
//source, both 0-based: 
a[0] <- source[0]
for i from 1 to n − 1 do
    j <- random integer with 0 ≤ j ≤ i
    a[i] <- a[j]
    a[j] <- source[i]

这个实现的正确性可以这么考虑:对于n-1次random调用生成的n!次随机序列,总会对应不 同的排列,所以n!次排列是均匀分布的。

这个实现有一个坏处,它需要额外的数组来存储生成的排列。但它也有一个大大的好处,它 可以不用预先知道数组的长度:

//To initialize an empty array a to a randomly shuffled copy of source whose
//length is not know:
while source.moreDataAvailable
    j <- random integer with 0 ≤ j ≤ a.length
    if j = a.length
        a.append(source.next)
    else
        a.append(a[j])
        a[j] <- source.next

Knuth-shuffle做些修改还可以有别的应用。比如下面这个问题:

  1. 如何从未知长度的序列中均匀地随机选择一个数?要求最多只遍历一遍。
  2. 如何从未知长度的序列中均匀地随机选择k个数?要求最多只遍历一遍。

.

问题(1)参考程序:

import random

def random_element(seq):
    """Return an element chosen at random from `seq`."""
    it = None
    for n, elem in enumerate(seq):
        if random.randint(0, n) == 0:
            it = elem
    return it

问题(2)参考程序:

import random

def random_sample(n, items):
    results = []

    for i, v in enumerate(items):
        r = random.randint(0, i)
        if r < n:
            if i < n:
                results.insert(r, v) # add first n items in random order
            else:
                #at a decreasing rate,replace random items
                results[r] = v 

    if len(results) < n:
        raise ValueError("Sample larger than population.")

    return results

2012年8月19日星期日

神奇的异或

熟练理解和掌握位操作,在有些情况下能使程序的效率提高不少,当然代码的可读性就降低 了,维护也比较困难。相比于与操作(&),或操作(|),异或操作(^)显得更难以理解,甚至看 不出它有什么妙用。下面列举一些用到异或操作的问题,或许能从中学到一些新的东西。
  1. 不使用中间变量,交换两个数的值
      //version 1
      void swap(int &a,int &b){
          a ^= b;
          b ^= a;
          a ^= b;
      }
    
      //version 2
      void swap(int &a,int &b){
          a -= b;
          b += a;
          a = b-a;
      }
    
    
    /*注意这两个版本都有个bug,当a,b指向同一个对象时,会变成0。而且有文章指出,版本
    1位操作未必比使用中间变量的性能好,再考虑到可读性,那这样的代码就更不推荐了。
    姑且就当作一种娱乐吧。
    */
    
  2. 有一组数字,每个数字都不相同,放在大小为N-1的数组里,数字的范围是1-N,找出丢失 的那个数。只能申请少于O(N)的常数空间,时间复杂度O(n)。如果丢失两个数呢?
    /*容易想到的O(n)算法是将所有的数相加,然后用N(N+1)/2去减,得到的差就是丢失的那个
    数。但问题在于所有的数之和可能会溢出。还有人说用 sum(array[i]-i),实际上这也没
    有解决可能的溢出问题。
    
    异或可以解决这个问题。将所有的数array[i]异或得到k,k再与1到N所有的数异或得到
    k',则k'就是丢失的那个数a。
    
    如果丢失两个数,同样用异或的方法。假设丢失的是a,b,则 k' = a^b,k'获得方法如上。
    定义lowbit(i)函数,返回i的二进制表示中不为0的最低位,则a,b在lowbit(k')位上肯定
    是不同的。将原数组array[1..N-1]中lowbit(k')位为1的数异或得到k'',再用k''与1-N中
    lowbit(k')位为1的数异或得到k''',则k'''为a,b其中的一个,再用k''' ^ k'得到另一个。
    */
    
      //查找唯一丢失的数
      //array 长度为N-1,值范围为1-N,数组中每个数都不相同
      int findmiss(int *array,int N){
          int a; // result
          int i;
          a = 0;
          for(i=0;i<N-1;i++)
              a ^= array[i];
          for(i=1;i<=N;i++)
              a ^= i;
          return a;
      }
      
      //lowbit 的实现
      int lowbit(int i){
          return i & -i;
      }
      
      //查找丢失的两个数
      //array 长度为N-2,值范围为1-N,数组中每个数都不相同
      //a,b为丢失的数
      void findmiss2(int *array,int N,int *a,int *b){
          int i,k,lb;
          k = 0;
          *a = 0,*b = 0;
          for(i=0;i<N-2;i++)
              k ^= array[i];
          /* after this recurse, k=a^b */
          for(i=1;i<=N;i++)
              k ^= i;
          lb = lowbit(k); 
          for(i=0;i<N-2;i++){
              if(array[i] & lb != 0)
                  *a ^= array[i];
              else
                  *b ^= array[i];
          }
          for(i=1;i<=N;i++){
              if(i & lb != 0)
                  *a ^= i;
              else
                  *b ^= i;
          }
          return ;
      }
    
  3. n个数中有且仅有一个数出现了奇数次(其他数都出现了偶数次),如果用线性时间常数 空间找出这个数?如果是两个数出现了奇数次呢?如果是三个数呢?或者是k个数?
    /*如果稍微想一想,就会发现这题与题目2有很多相似之处。所以下面的解法也很类似。
      
    第一问:有且仅有一个数出现奇数次,对数组n个数从头到尾异或一遍,得到的值即为要
    求的那个数。
    
    第二问:如果有两个数出现了奇数次。同题目2,先异或所有元素得到k,取lowbit(k),
    按照lowbit(k)将原数组分为两组,对这两组分别求异或,就得到要求的两个数。
    
    第三问: 如果是三个数。这一问就很难了,即使知道应该从异或这个角度去考虑,也很难
    想出好方法。当然,总有牛人会想出好方法,关键是找出第一个来,然后借助第二问结论
    求另外两个。为了不影响意思的理解,转载原文如下:
    
    // Let s = a ^ b ^ c.  We know that s not in (a, b, c), 
    // since if s == a, say, then b ^ c == 0 and b == c.  
    // Let f(x) be the lowest bit where x differs from s.  
    // The algorithm computes flips = f(a) ^ f(b) ^ f(c), 
    // since the numbers appearing an even number of times cancel.  
    // The variable flips has parity 1, so it is non-zero, 
    // and lowbit(flips) is a bit where one or three of a, b, c 
    // differ from s.  It can't be three(这个需要仔细想想), however, so the final 
    // exclusive-or includes exactly one of a, b, c.
    
    第四问:如果是k个数。从第3问我们可以发现,用类似的方法总可以找出第一个来,然
    后将问题转化为k-1规模的,则可以递归解决。复杂度为O(k^2*n).
    */
    
      // get lowest different bit
      int lowbit(int x){
          return x & ~(x - 1);
      }
      
      // Given an array of n integers, such that each number 
      // in the array appears exactly twice, except for two 
      // numbers (say a and b) which appear exactly once.
      // 
      // In O(n) time and O(1) space find a and b. 
      // e.g.
      // { 2 3 3 2 4 6 4 7 8 8 }  ---> a/b = { 6 7}
      void Find2(int seq[], int n, int& a, int& b)
      {
          ////XOR
          int xors = 0;
          for(int i = 0; i < n; i++)
              xors ^= seq[i];
      
          ////get different bit
          int diff = lowbit(xors);
      
          ////
          a = 0;
          b = 0;
          for(int i = 0; i < n; i++){
              if(diff & seq[i])
                  a ^= seq[i];
              else
                  b ^= seq[i];
          }
      }
      
      
      // Given an array of n integers, such that each number 
      // in the array appears exactly twice, except for three 
      // numbers (say a, b and c) which appear exactly once.
      // 
      // In O(n) time and O(1) space find a,b and c. 
      // e.g.
      // { 2 3 3 2 4 6 4 7 8 8 1 }  ---> a/b = { 6 7 1}
      
      
      // Let s = a ^ b ^ c.  We know that s not in (a, b, c), 
      // since if s == a, say, then b ^ c == 0 and b == c.  
      // Let f(x) be the lowest bit where x differs from s.  
      // The algorithm computes flips = f(a) ^ f(b) ^ f(c), 
      // since the numbers appearing an even number of times cancel.  
      // The variable flips has parity 1, so it is non-zero, 
      // and lowbit(flips) is a bit where one or three of a, b, c 
      // differ from s.  It can't be three, however, so the final 
      // exclusive-or includes exactly one of a, b, c.
      
      void Find3(int seq[], int n, int& a, int& b, int& c)
      {
          ////XOR
          int xors = 0;
          for(int i = 0; i < n; i++)
              xors ^= seq[i];
      
          ////
          int flips = 0;
          for(int i = 0; i < n; i++)
              flips ^= lowbit(xors ^ seq[i]);
          flips = lowbit(flips);
      
          ////get one of three
          a = 0;
          for(int i = 0; i < n; i++){
              if(lowbit(seq[i] ^ xors) == flips)
                  a ^= seq[i];
          }
      
          ////swap a with the last element of seq
          for(int i = 0; i < n; i++){
              if(a == seq[i]){
                  int temp = seq[i];
                  seq[i] = seq[n - 1];
                  seq[n - 1] = temp;
              }
          }
      
          ////call Find2() to get b and c
          Find2(seq, n - 1, b, c);
      }
    
  4. n个数中只有一个数出现了一次,其他的数出现了3次,如何用O(n)的时间复杂度,常数空 间的算法找出这个数?如果这个数出现两次呢?
    /*考虑3进制。利用3进制的异或:  
    1 xor 0 = 1
    1 xor 1 = 2
    1 xor 2 = 0
    2 xor 0 = 2
    2 xor 2 = 1
    0 xor 0 = 0
    
    详见参考资料[4]
    */
    
  5. nim游戏.有若干堆石子,每次可以选择一堆石子,从这堆石子中拿走任意数量的石 子,也就是至少拿走一个,最多把这堆石子全部拿走,两人轮流取,谁取走最后一个石子 谁就赢.问是否先手必胜。
    假设有n堆石子,数量分别为a1,a2,a3...an,那么如果石子的异或和不是0, 那么先手必
    胜也就是a1^a2^...^an != 0,那么先手必胜。看起来很神奇吧,感觉很复杂的游戏居然用
    异或就给解决了。是不是觉得自己很笨,怎么没有想到?没啥关系,早在16世纪,nim游
    戏就已经提出来了,但是到了1901年才被哈佛的一个教授解决,说明这世上跟我一样笨的
    人还是很多的。详细的讨论可以看 "http://en.wikipedia.org/wiki/Nim" .那先手要怎么走呢(想想怎么使异或变为0)?
    

2012年6月20日星期三

String Reduction

来源:Tristan's Collection of Interview Questions

Problem: A string that only contains three types of letters :"a", "b", "c". Our goal is to reduce the string. If two different letters are adjacent to each other, we can reduce the two letters into another letter. The rule is like this: if we encounter "ab" (or "ba"), we reduce to "c"; if we encounter "bc" (or "cb"), we reduce to "a"; if we encounter "ca" (or "ac"), we reduce to "b". If two adjacent letters are the same, we can't reduce them. We need to find an optimal way to reduce the string as much as possible.

For example, for string "acbb", if we try to reduce "ac"first, we get "bbb". But, if we reduce "cb" first, we get "aab" and we can further make "aab" => "ac" => "b". The latter approach gives us an optimal reduction.What is the smallest string which can result by applying this operation repeatedly?























Solution: Here what really matters is the numbers of letter "a", "b" and "c". Let us name them as num_a, num_b and num_c. If two of them are zero, we can't reduce that string. Otherwise, we can reduce the string into a string that has a length of 1 or 2. If num_a, num_b and num_c are all even or odd, we can reduce to a string with length 2; If not, we can reduce to a string with length 1.

Then how to do the reduction? The detail is as follow:
  1. if the string has a length of 3 and contains one "a", one "b" and one "c", whatever order of the three letter are in, we can only reduce the string into a string with length 2; if the string has a length of 2 and contains two different letters, we can only reduce the string into a string with length 1. Let us regard these two cases as base cases.

  2. for a general case, we have num_a "a" , num_b "b" and num_c "c". After every reduction, the sum of num_a, num_b and num_c will decrease by 1, since we substitue two letters with one letter.

  3. Then at each round, which two adjacent letters we choose to reduce? We try to find an adjacent pair that contains a letter which has the highest count. For example, if now, we have 3 "a", 4 "b" and 6 "c" in the string, we choose an adjacent pair that contains "c"since num_c = max(num_a, num_b, num_c) . Can we find such pair? definitely we can. If there are multiple such pairs, choose a random one. Then after we reduce this pair, num_c--, max(num_a, num_b, num_c) may decrease by 1, remain unchanged, or increase by 1. However, since max(num_a, num_b, num_c) is up-bounded by num_a + num_b + num_c. num_a + num_b + num_c is decreasing after every round, then max(num_a, num_b, num_c) will also decrease if we look at a long wrong. Therefore, by using our strategy, max(num_a, num_b, num_c) will eventually decrease to 1, which means we are encounter the base cases in step 1.

  4. Then when the string will be reduced to a length of 1 and when to a length of 2? We observe that is num_a, num_b and num_c are all even, then after one transformation, they will become all odd; similarly, if there are originally all odd, after one transformation, they will become all even. Then according to the analysis in step 3), we know at the end, the max(num_a, num_b, num_c) will eventually decrease to 1. But, they should still be all odd at that time (since "1" is odd). Therefore, at the very end, we will have num_a = 1, num_b = 1 and num_c =1, which will eventually lead to a string of length 2. It is easy to prove that if num_a, num_b and num_c are not all even or odd, the string will be reduced to length 1.

  5. if num_a, num_b and num_c are not all even or odd, there must be only two cases: a) two of the counters are odd and one is even b) two of the counters are even and one is odd. For example, if num_b is odd, num_a and num_c are both even. The string actually will eventually be reduced to "b".

  6. if num_a, num_b and num_c are all even or odd, there might be multiple different final forms.
if((a == 0 && b == 0 && c == 0) ||
   (a == 0 && b == 0 && c != 0) ||
   (a == 0 && b != 0 && c == 0) ||
   (a != 0 && b == 0 && c == 0))
   {
      result = a+b+c;
   }
else if(a != 0 && b != 0 && c != 0)
   {
      if((a%2 == 0 && b%2 == 0 && c%2 == 0) ||
         (a%2 == 1 && b%2 == 1 && c%2 == 1))
           result = 2;
      else
           result = 1;
   }
else if((a == 0 && b != 0 && c != 0) || 
        (a != 0 && b == 0 && c != 0) || 
        (a != 0 && b != 0 && c == 0))
   {
      if(a%2 == 0 && b%2 == 0 && c%2 == 0)
           result = 2;
      else
           result = 1;
   }

2012年6月14日星期四

树状数组

一个基本的问题:求数组中的逆序对数目?

这也是算法导论第一章的最后一道习题。按照算法导论中的提示,该问题可以通过修改合并排序在O(nlgn)内解决。另外一种方法,是使用一种特殊的数据结构称之为“树状数组”,也可以在O(nlgn)内解决。代码如下:

#include<iostream>
#include<stdio.h>
#include<string.h>
using namespace std ;
#define MAXN 100002
#define MAX 1000002
int n,a[MAXN],c[MAX] ;
int main()
{
 int runs ;
 scanf("%d",&runs) ;
 while(runs--)
 {
  scanf("%d",&n) ;
  for(int i = 0;i < n;i++) scanf("%d",&a[i]) ;
  long long ret = 1LL * n * (n - 1) / 2 ;
  memset(c,0,sizeof c) ;
  for(int i = 0;i < n;i++)
  {
   for(int j = a[i];j > 0;j -= j & -j) ret -= c[j] ;
   for(int j = a[i];j < MAX;j += j & -j) c[j]++ ;  
  }
  cout << ret << endl ;
 }
 return 0 ;
}


下面介绍下树状数组(Binary Indexed Tree)。如果搜一下,就会发现有很多文章详细介绍树状数组,比如这篇,所以这里不再赘述。 它能够高效地获取数组中连续n个数的和。概括说,树状数组通常用于解决以下问题:数组{a}中的元素可能不断地被修改,怎样才能快速地获取连续几个数的和? 传统数组(共n个元素)的元素修改和连续元素求和的复杂度分别为O(1)和O(n)。树状数组通过将线性结构转换成伪树状结构(线性结构只能逐个扫描元素,而树状结构可以实现跳跃式扫描),使得修改和求和复杂度均为O(lgn),大大提高了整体效率。 let the code talk.值得提一下的是,树状数组下标是从1开始的。




//求2^k
int lowbit(int t)
{
    return t & ( t ^ ( t - 1 ) );
}
 
//求前n项和
int sum(int end)
{
   int sum = 0;
   while(end > 0)
  {
     sum += in[end];
     end -= lowbit(end);
  }
 
  return sum;
}
 
//增加某个元素的大小
void plus(int pos, int num)
{
   while(pos <= n)
  {
     in[pos] += num;
     pos += lowbit(pos);
  }
}

2012年2月10日星期五

euler 48: 大整数乘法

euler project 48:

The series, 11 + 22 + 33 + ... + 1010 = 10405071317.
Find the last ten digits of the series, 11 + 22 + 33 + ... + 10001000.

我是用python解的,python 的数没有大小的限制,这样就很方便,brute force.

在问题的讨论帖中,有人(grimbal)贴了他用C解决的算法。我觉得很优美,记录如下:

#define mod 10000
void ex48(){
    int n, k;
    unsigned int p2, p1, p0;
    unsigned int s2, s1, s0;
    unsigned int t2, t1, t0;

    // s = 0
    s2 = s1 = s0 = 0;
    for( n=1 ; n<=1000 ; n++ ){
        // p = 1
        p2 = p1 = 0; p0 = 1;
        for( k=1 ; k<=n ; k<<=1 );
        while( k>>=1 ){
            // p *= p
            t0 = p0*p0;
            t1 = 2*p0*p1 + t0/mod;
            t2 = 2*p0*p2 + p1*p1 + t1/mod;
            p0 = t0%mod;
            p1 = t1%mod;
            p2 = t2%mod;
            if( n&k ){
                // p *= n
                t0 = p0*n;
                t1 = p1*n + t0/mod;
                t2 = p2*n + t1/mod;
                p0 = t0%mod;
                p1 = t1%mod;
                p2 = t2%mod;
            }
       }
       // s += p
       t0 = p0+s0;
       t1 = p1+s1 + t0/mod;
       t2 = p2+s2 + t1/mod;
       s0 = t0%mod;
       s1 = t1%mod;
       s2 = t2%mod;
    }

    printf("ex48: %02d%04d%04d\n", s2%100, s1, s0);
}

2012年1月18日星期三

euler 3: Largest prime factor

题出自 Project Euler Problem 3:
The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?

我自己的解法思路是先判断该数(记为A)是否为prime,如果是,则直接返回它本身。否则从该数的一半开始(记为index),判断index是否能被A整数以及index是否是prime,当两者同时满足时,返回index,否则index自减1。 判断某数是否为prime的算法,是从该数的平方根数开始一个一个检查。代码如下:

#include <stdio.h>
#include <math.h>

//#include <limits.h>

int is_prime(long l);
long largest_prime_factor(long l);

int main(int argc, char *argv[])
{
  long result = largest_prime_factor(600851475143L);
  printf("%ld\n",result);
  return 0;
}

int is_prime(long l){
  long l_sqrt = (long)sqrt(l);
  long i = 2L;
  for(; i <= l_sqrt; i++){
    if(l % i == 0)
      return 0;                 /* is not prime */
  }

  /* is prime */
  return 1;
}


long largest_prime_factor(long l){
  if(is_prime(l))
    return l;
 
  long l_mid = l/2;
  long l_ind;
  for(l_ind = l_mid; l_ind >=2; l_ind--){
    if((l%l_ind == 0) && is_prime(l_ind))
      return l_ind;
  }

  return 1; // never come here
}



上面的算法很惨,算了几十分钟都没有算出结果。所以最后我放弃了,google到了下面这个算法:

#include <stdio.h>

long largest_prime_factor(long l);

int main(int argc, char *argv[])
{
  long result = largest_prime_factor(600851475143L);
  printf("%ld\n",result);
  return 0;
}

long largest_prime_factor(long l){
  long factor = 2L;
  while(l > 1L){
    while(l % factor == 0 ){
      l = l / factor;
    }
    factor++;
  }

  return (factor-1);
}

细想一下,这个算法与筛选法类似。但是,我怎么就没想到呢?FAIL!!!

2012年1月1日星期日

判断无符号整数二进制中1的个数(外一则)


三种算法。

算法1:
比较容易想到的算法,逐位比较,最多需要32次。
int count_ones_01(unsigned int n)
{
    int i;
    for(i = 0;n > 0;n >>= 1)
         i += (n&1);
    return i;
}

算法二:
Monster级算法。比较次数比算法1少,最多需要32次(如果int长度为32位的话)。
int count_ones_02(unsigned int n){
  int count=0;
  while(n>0){
    n &= (n-1);
    count++;
  }
  return count;
}

算法三:
二分法。遗憾的是不具有跨平台性。运算次数为5次(\(\log_{2}{32}\))
int count_ones(unsigned int n)
{
    n = (n & 0x55555555) + ((n >> 1) & 0x55555555);
    n = (n & 0x33333333) + ((n >> 2) & 0x33333333);
    n = (n & 0x0f0f0f0f) + ((n >> 4) & 0x0f0f0f0f);
    n = (n & 0x00ff00ff) + ((n >> 8) & 0x00ff00ff);
    n = (n & 0x0000ffff) + ((n >> 16) & 0x0000ffff);

    return n;
}

算法的思路是这样的:2位为一组,相加,看看有几个1。再4位为一组,相加,看看有几个1,
以此类推。
为了简单说明,先看看8位的情形。
x = (x & 0x55) + ((x >> 1) & 0x55);    (1)
x = (x & 0x33) + ((x >> 2) & 0x33);    (2)
x = (x & 0x0f) + ((x >> 4) & 0x0f);    (3)
return x;

对于语句(1)

当x=11111111
x= (x & 0x55) + ((x >> 1) & 0x55);
之后x=10101010。你2位2位地看,10=2, 就是2 2 2
2, 就是说各组都是2个1。

当x=00101001
x= (x & 0x55) + ((x >> 1) & 0x55);
之后x=00010101。你2位2位地看,就是0 1 1 1, 前
1组只有0个1,后面的组都是1个1。

接着看语句(2)

当x=abcdefgh
x=(x & 0x33)+((x >> 2)&0x33);
相当于 00cd00gh + 00ab00ef。

因为语句(1)之后,每两位表示该两位有几个1,相加就指示每4位有多少个1。这样就是4位
4位为一组。注意这样的分组,组与组之间永远都不会产生进位的。正因为不会产生进位,才
可以分开来看。

8位,16位,32位依上类推。

外一则:
利用算法三类似的思想,进行反转一个字节的二进制次序。
unsigned char reverse8( unsigned char c )
{
     c = ( c & 0x55 ) << 1 | ( c & 0xAA ) >> 1;
     c = ( c & 0x33 ) << 2 | ( c & 0xCC ) >> 2;
     c = ( c & 0x0F ) << 4 | ( c & 0xF0 ) >> 4;
     return c;
}

2011年12月28日星期三

最大的子序列和问题


最大的子序列和问题描述:
给定整数A1,A2,…,AN,求SUM(Ak),k从i到j,的最大值。为了方便起见,如果所有整数均为
负数,则最大子序列和为0。

例:
输入-2,11,-4,13,-5,-2时,答案为20(从A2到A4)

最大子序列和问题的解。
算法一:
穷举,复杂度为O(N^3)
int maxsubsequencesum(const int A[],int N){
    int thissum,maxsum,i,j,k;
   
    maxsum=0;
    for(i=0;i<N;i++)
        for(j=i;j<N;j++){
                thissum=0;
                for(k=i;k<=j;k++)
                        thissum += A[k];
               
                if(thissum > maxsum)
                        maxsum = thissum;
        }
    return maxsum;
}
算法二:
对算法一进行优化,减少内循环,复杂度为O(N^2)
int maxsubsequencesum_2(const int A[],int N){
    int thissum,maxsum,i,j;
   
    maxsum=0;
    for(i=0;i<N;i++){
        thissum=0;
        for(j=i;j<N;j++){
             thissum += A[j];
                    
             if(thissum > maxsum)
                   maxsum = thissum;
        }
    }
    return maxsum;
}
算法三:
分治:最大子序列可能在三处出现,或者整个出现在数据的左半部,或者整个出现在右半部,
或者跨越输入数据的中部而占据左右两半部分。前两种情况可以递归求解,第三种情况的最
大和可以通过求出前半部分的最大和(包含前半部分的最后一个元素)以及后半部分的最大
和(包含后半部分的第一个元素),将这两个和相加而得到。
复杂度为O(NlgN)
int maxsubsum(const int A[], int left,int right){
    int maxleftsum,maxrightsum;
    int maxleftbordersum,maxrightbordersum;/*包含边界的最大值*/
    int leftbordersum,rightbordersum;
    int center,i;

    if(left == right){
         if(A[left] > 0)
             return A[left];
         else
             return 0; /*由题给条件,负数的最大子序列长度为0*/
    }

    /* divide */
    center = (left + right)/2;
    maxleftsum = maxsubsum(A,left,center);
    maxrightsum = maxsubsum(A,center+1,right);

   /*conquer*/
   maxleftbordersum=0;
   leftbordersum=0;
   for(i=center; i>=left;i--){
       leftbordersum += A[i];
       if(leftbordersum > maxleftbordersum)
          maxleftbordersum = leftbordersum;
   }

  maxrightbordersum=0;
  rightbordersum=0;
  for(i=center+1; i <= right; i++){
      rightbordersum += A[i];
      if(rightbordersum > maxrightbordersum)
          maxrightbordersum = rightbordersum;
  }

  /* max3 取三者最大值 */
  return max3(maxleftsum, maxrightsum, maxleftbordersum + maxrightbordersum);
}

int maxsubsequencesum_3(const int A[], int N){
    return maxsubsum(A,0,N-1);
}
算法四:
在线算法,复杂度O(N),
int maxsubsequencesum(const int A[],int N){
    int thissum,maxsum,j;
    thissum = maxsum = 0;

    for(j=0; j < N; j++){
        thissum += A[j];
       
        if(thissum > maxsum)
            maxsum = thissum;
        else if(thissum < 0)
            thissum = 0;
    }
    return maxsum;
}

2011年12月7日星期三

POJ 1009

复制至此有两个原因:
1. 这题我看完毫无头绪,佩服会解该题的人。
2. 下面的代码写的很干净

原文地址:这里


只需考虑原图中数值变化的点,其他点的编码与其左侧的点相同。
最开始用brute force,优化了输出编码中有连续0的情况,sample input可以较快通过,但提交结果超时。又换成上面的方法,可以AC,但代码改得很乱,待重写。


/*
 *  只需考虑原图中数值发生变化的点
 */
#include <stdio.h>
#include <math.h>

typedef struct _Ref {
    int index;
    int code;
} Ref;

int IN[1000][2];
Ref OUT[8000];

int encode(int n, int width, int total);
int getcode(int n);

void sort(Ref *list, int n);

int
main()
{
    int width,i,j,m,n,t,k,cnt,total,current,row,col;
    Ref cur;
   
    while(scanf("%d", &width) != EOF) {
        if(width == 0)
            break;
       
        // input  
        i = total = 0;
        scanf("%d%d", &IN[i][0], &IN[i][1]);
        while(IN[i][0] != 0 || IN[i][1] != 0) {
            total += IN[i][1];
            i++;
            scanf("%d%d", &IN[i][0], &IN[i][1]);  
        }
       
        printf("%d\n", width);
       
        cnt = i;
        n = 1;
        k = 0;
        for(m=0; m<=cnt; m++) {
            row = (n - 1) / width;
            col = (n - 1) % width;
           
            for(i=row-1; i<=row+1; i++) {
                for(j=col-1; j<=col+1; j++) {
                    t = i * width + j;
                    if(i<0 || j<0 || j>width-1 || t>total-1)
                        continue;
                   
                    OUT[k].index = t + 1;
                    OUT[k++].code = encode(t+1, width, total);                    
                }
            }
            n += IN[m][1];
        }
       
        sort(OUT, k);
       
        cur = OUT[0];
        for(i=0; i<k; i++) {
            if(OUT[i].code == cur.code)
                continue;
               
            printf("%d %d\n", cur.code, OUT[i].index - cur.index);
            cur = OUT[i];
        }
       
        printf("%d %d\n", cur.code, total - cur.index + 1);
        printf("0 0\n");
    }
    printf("0\n");
   
    return 0;
}

// 返回第n个元素的编码,n从1开始
int
getcode(int n)
{
    int t,i;
   
    i = t = 0;
    while(t < n) {
        t += IN[i++][1];
    }
   
    return IN[i-1][0];
}

// 对第n个元素进行编码,n从1开始
int
encode(int n, int width, int total)
{
    int i,j,t,code,result,row,col;
   
    code = getcode(n);
   
    row = (n - 1) / width;
    col = (n - 1) % width;
   
    result = 0;
    for(i=row-1; i<=row+1; i++) {
        for(j=col-1; j<=col+1; j++) {
            t = i * width + j;
            if(i<0 || j<0 || j>width-1 || t == n-1 || t>total-1)
                continue;
               
            if(abs(getcode(t+1) - code) > result)
                result = abs(getcode(t+1) - code);
        }
    }
   
    return result;
}

void sort(Ref *list, int n)
{
    Ref tmp;
    int i,j;
   
    for(i=0; i<n-1; i++) {
        for(j=i+1; j<n; j++) {
            if(list[i].index > list[j].index) {
                tmp = list[j];
                list[j] = list[i];
                list[i] = tmp;
            }
        }
    }
}

2011年11月8日星期二

如何心算是星期几


这个心算算法的技巧在于利用了这么一个结论:对于任意给定的一年,某些特定的日期总是属于相同的星期几。我们称这些日期为“锚点“。方便记忆的锚点有:5月9号,9月5号,7月11号,11月7号,4月4号,6月6号,8月8号,10月10号,12月12号,以及2月的最后一天(平年时为28号,闰年时为29号)。前面四个日期可以通过"在7-11朝9晚5地工作"这句顺口溜来记忆。

锚点在给定的年份属于星期几是不固定的,但变化是有规律可循的。比如今年2011年,锚点都是星期一,而在2010年,锚点都是星期日。每过一年,锚点星期就往前移一天。这很好理解,每年有365天,365 % 7 = 1,所以锚点星期会移一天。因为闰年有366天,所以当年是闰年的话,跟闰年的上一年比较,锚点是向前移了两天的。比如2012年是闰年,所以锚点都是星期三。

这样只要记住一些锚点的星期,就可以很快心算出任意日期的星期了。


参考资料

2011年4月11日星期一

一次谷歌面试趣事

本文原始地址:一次谷歌面试趣事


很多年前我进入硅谷人才市场,当时是想找一份高级工程师的职位。如果你有一段时间没有面试过,根据经验,有个非常有用的提醒你应该接受,就是:你往往会在前几次面试中的什么地方犯一些错误。简单而言就是,不要首先去你梦想的公司里面试。面试中有多如牛毛的应该注意的问题,你可能全部忘记了,所以,先去几个不太重要的公司里面试,它们会在这些方面对你起教育(再教育)作用。

我第一家面试的公司叫做gofish.com,据我所知,gofish这家公司如今的情况跟我当时面试时完全的不同。我几乎能打保票的说,当时我在那遇到的那些人都已不再那工作了,这家公司的实际情况跟我们这个故事并不是很相关。但在其中的面试却是十分相关的。对我进行技术性面试的人是一个叫做Guy的家伙。

Guy穿了一条皮裤子。众所周知,穿皮裤子的面试官通常是让人"格外"恐怖的。而Guy也没有任何让人失望的意思。他同样也是一个技术难题终结者。而且是一个穿皮裤子的技术难题终结者 —— 真的,我做不到他那样。

我永远不会忘记他问我的一个问题。事实上,这个问题是非常的普通 —— 在当时也是硅谷里标准的面试题。

问题是这样的:

假设这有一个各种字母组成的字符串,假设这还有另外一个字符串,而且这个字符串里的字母数相对少一些。从算法是讲,什么方法能最快的查出所有小字符串里的字母在大字符串里都有?

比如,如果是下面两个字符串:

String 1: ABCDEFGHLMNOPQRS

String 2: DCGSRQPOM

答案是true,所有在string2里的字母string1也都有。如果是下面两个字符串:

String 1: ABCDEFGHLMNOPQRS

String 2: DCGSRQPOZ

答案是false,因为第二个字符串里的Z字母不在第一个字符串里。

当他问题这个问题时,不夸张的说,我几乎要脱口而出。事实上,对这个问题我很有信心。(提示:我提供的答案对他来说显然是最糟糕的一种,从面试中他大量的各种细微表现中可以看出来)。

对于这种操作一种幼稚的做法是轮询第二个字符串里的每个字母,看它是否同在第一个字符串里。从算法上讲,这需要O(n*m)次操作,其中n是string1的长度,m是string2的长度。就拿上面的例子来说,最坏的情况下将会有16*8 = 128次操作。

一个稍微好一点的方案是先对这两个字符串的字母进行排序,然后同时对两个字串依次轮询。两个字串的排序需要(常规情况)O(m log m) + O(n log n)次操作,之后的线性扫描需要O(m+n)次操作。同样拿上面的字串做例子,将会需要16*4 + 8*3 = 88加上对两个字串线性扫描的16 + 8 = 24的操作。(随着字串长度的增长,你会发现这个算法的效果会越来越好)

最终,我告诉了他一个最佳的算法,只需要O(n+m)次操作。方法就是,对第一个字串进行轮询,把其中的每个字母都放入一个Hashtable里(成本是O(n)或16次操作)。然后轮询第二个字串,在Hashtable里查询每个字母,看能否找到。如果找不到,说明没有匹配成功。这将消耗掉8次操作 —— 这样两项操作加起来一共只有24次。不错吧,比前面两种方案都要好。

Guy没有被打动。他把他的皮裤子弄的沙沙响作为回应。"还有没有更好的?"他问道。

我的天?这个家伙究竟想要什么?我看看白板,然后转向他。"没有了,O(n+m)是你能得到的最好的结果了 —— 我是说,你至少要对每个字母至少访问一次才能完成这项操作 —— 而这个方案是刚好是对每个字母只访问一次"。我越想越确信我是对的。

他走到白板前,"如果这样呢 —— 假设我们有一个一定个数的字母组成字串 —— 我给每个字母分配一个素数,从2开始,往后类推。这样A将会是2,B将会是3,C将会是5,等等。现在我遍历第一个字串,把每个字母代表的素数相乘。你最终会得到一个很大的整数,对吧?然后 —— 轮询第二个字符串,用每个字母除它。如果除的结果有余数,这说明有不匹配的字母。如果整个过程中没有余数,你应该知道它是第一个字串恰好的子集了。这样不行吗?"

每当这个时候 —— 当某个人的奇思异想超出了你的思维模式时,你真的需要一段时间来跟上他的思路。现在他站在那里,他的皮裤子并没有帮助我理解他。

现在我想告诉你 —— Guy的方案(不消说,我并不认为Guy是第一个想出这招的人)在算法上并不能说就比我的好。而且在实际操作中,你很可能仍会使用我的方案,因为它更通用,无需跟麻烦的大型数字打交道。但从"巧妙水平"上讲,Guy提供的是一种更、更、更有趣的方案。

我没有得到这份职位。也许是因为我拒绝了他们提供给我的一些讨厌的工作内容和其它一些东西,但这都无所谓了。我还有更大更好的目标呢。

接着,我应聘了become.com。在跟CTO的电话面试中,他给我布置了一道"编程作业"。这个作业有点荒唐 —— 现在回想起来,大概用了我3天的时间去完成。我得到了面试,得到了那份工作 —— 但对于我来说,最大的收获是这道编程作业强迫我去钻研并有所获。我需要去开发一个网页爬虫,一个拼写检查/纠正器,还有一些其它的功能。不错的东西。然而,最终,我拒绝了这份工作。

终于,我来到了Google面试。我曾说过Google的面试过程跟外面宣传的很一致。冗长 —— 严格,但诚实的说,相当的公平。他们在各种面试过程中尽最大的努力去了解你、你的能力。并不是说他们在对你做科学研究,但我确信他们是努力这样做。

我在Google的第四场面试是一个女工程师,老实话,是一场很无聊的面试。在前面几场面试中我表现的很好,感觉到我的机会非常的大。我相信如果不做出什么荒唐事情来,十拿九稳我能得到这份工作。

她问了我一些关于排序或设计方面的非常简单的问题,我记不清了。但就在45分钟的面试快要结束时,她对我说"我还有一个问题。假设你有一个一定长度的由字母组成的字符串。你还有另外一个,短些。你如何才能知道所有的在较短的字符串里的字母在长字符串里也有?"

哇塞。Guy附身了。

现在,我完全可以马上结束这场面试。我可以对她说"哈哈,几个星期前我就知道答案啦!",这是事实。但就是在几个星期前被问到这个问题时 —— 我给出的也是正确的答案。这是我本来就知道答案的问题。看起来就好像是Guy为我的这次面试温习过功课一样。而且,可恶,人们通常是通过上网来搜集面试问题 —— 而我,我可以毫不客气的说,对于这些问题,我不需要任何"作弊"。我自己知道这些答案!

现在你们可能认为——就在她问出了问题之后,在我准备开始说出在脑海里构思完成的最后的演讲之前——你们可能会想,我应该是,当然该,从情理上讲,镇定的回答出这个问题,并且获得赞赏。可糟糕的是,事实并不是这样。打个比喻,就像是她问出来问题后,我在闹子里立即举起了手,并大叫着"我!嗨!嗨!我知道!让我来回答吧!"我的大脑试图夺走我对嘴巴的控制权(这事经常发生),幸亏我坚强的毅力让我镇定下来。

于是我开始回答。平静的。带着不可思议的沉着和优雅。带着一种故意表现出来的 —— 带着一种,我认为,只有那种完全的渊博到对古今中外、不分巨细的知识都精通的人才能表现出来的自信。

我轻描淡写的说出来那种很幼稚的方案,就好象是这种方案毫无价值。我提到了给它们排序,就好像是在给早期的《星际迷航》中的一个场景中的人物穿上红T恤似的。最后,平淡的,就好像是我决定了所有事情的好坏、算法上的效率,我说出了O(n+m)一次性方案。

我要告诉你——尽管我表明上的平静——这整个过程我却在做激烈的挣扎,内心里我在对自己尖着——"你个笨蛋,赶紧告诉她素数方案!"

当我完成了对一次性算法的解释后,她完全不出意外的认可的点了下头,并开始在笔记本上记录。这个问题她以前也许问过了一百次,我想大部分的人都能回答上来。她也许写的是"回答正确。无聊的面试。正确的回答了无聊的字符串问题。没有惊喜。无聊的家伙,但可以留下。"

我等了一会。我让这种焦灼的状态持续的尽可能的长。我可以发誓的说,如果再耽搁一分钟,我一定会憋出脑血栓、脱口说出关于素数的未解之谜。

我打破了沉默。"你知道吗,还有另外一个,可能是更聪明的算法。"

她二目空空的抬头看了一眼,仅在瞬间闪现过一丝希望。

"假设我们有一定长度的字符串。我们可以给每个字母分配一个素数,从2开始。然后我们把大字串中的每个字母代表的素数相乘得出一个数,用小字串中的每个字母代表的素数去除它。如果除的过程中没有产生余数,则小字串是大字串的一个子集。"

在此时,我猜,她看起来就像是Guy当时把相同的话说给我听时我表现出来的样子。而我演讲时泰然自若的表情没了,眼睛瞪大,说话时稍微带出来一些唾沫星子。

一会儿后,她不得不说了,"可是…等一下,有可能…是的,可以这样!可是如何…如果…噢,噢,可行!简洁!"

我得意洋洋的吸了一口气。我在我的面试记录里写下了"她给了我一个'简洁'的评语!"在她提出这个问题之前我就确信能得到这份工作,现在我更加确信了。还有一点我十分确信的是,我(更准确的说是Guy)给了她今天的好心情。

我在Google干了3年,生活的十分愉快。我在2008年辞职去到一个小公司里做CTO,之后又开办了一个自己的公司。大概是一年前,我偶然的在一个创业论坛会上遇到了Guy,他记不得我了,当我向他细述这段往事时,他对他那条皮裤子大笑不已。

话说回来,如果这个故事里有什么教育意义的话——永远不要冒失的首先去应聘你梦想的公司,应先去应聘那些你不看好的职位。你除了能从这些面试中获得经验外,你指不定能遇到某个能为你的更重要的面试铺路的人呢。事实上,这个经验在你生活中的很多其它事情上也适应。

说正经的,如果你有机会想找一个解决问题的高手——雇佣Guy比谁都强。那个家伙很厉害。

(在这些陈年旧账里发现的一点技术瑕疵:字母有可能重复而字符串可能会很长,所以必须要有统计。用那个最幼稚的解决方案时,当在大字符串里找到一个字符后就把它删掉,当这样仍然是 O(n*m)次。在Hashtable里我们会有一个key->value的计数。Guy的方案在这种情况下仍然好用。)

修改:11/30/10 —— 本文中提到的Guy看到了这篇文章,并在评论中做了澄清。值得一读。


2011年4月3日星期日

如何计算某一天是星期几?—— 蔡勒(Zeller)公式

来源:原文地址




历史上的某一天是星期几?未来的某一天是星期几?关于这个问题,有很多计算公式(两个通用计算公式和一些分段计算公式),其中最著名的是蔡勒(Zeller)公式。即w=y+[y/4]+[c/4]-2c+[26(m+1)/10]+d-1



公式中的符号含义如下,w:星期;c:世纪-1;y:年(两位数);m:月(m大于等于3,小于等于14,即在蔡勒公式中,某年的1、2月要看作上一年的13、14月来计算,比如2003年1月1日要看作2002年的13月1日来计算);d:日;[ ]代表取整,即只要整数部分。(C是世纪数减一,y是年份后两位,M是月份,d是日数。1月和2月要按上一年的13月和 14月来算,这时C和y均按上一年取值。)算出来的W除以7,余数是几就是星期几。如果余数是0,则为星期日。



以2049年10月1日(100周年国庆)为例,用蔡勒(Zeller)公式进行计算,过程如下:

蔡勒(Zeller)公式:w=y+[y/4]+[c/4]-2c+[26(m+1)/10]+d-1

=49+[49/4]+[20/4]-2×20+[26× (10+1)/10]+1-1

=49+[12.25]+5-40+[28.6]

=49+12+5-40+28

=54 (除以7余5)

即2049年10月1日(100周年国庆)是星期5。



你的生日(出生时、今年、明年)是星期几?不妨试一试。



不过,以上公式只适合于1582年10月15日之后的情形(当时的罗马教皇将恺撒大帝制订的儒略历修改成格里历,即今天使用的公历)。过程的推导:(对推理不感兴趣的可略过不看)

星期制度是一种有古老传统的制度。据说因为《圣经·创世纪》中规定上帝用了六天时间创世纪,第七天休息,所以人们也就以七天为一个周期来安排自己的工作和生活,而星期日是休息日。从实际的角度来讲,以七天为一个周期,长短也比较合适。所以尽管中国的传统工作周期是十天(比如王勃《滕王阁序》中说的“十旬休暇”,即是指官员的工作每十日为一个周期,第十日休假),但后来也采取了西方的星期制度。



在日常生活中,我们常常遇到要知道某一天是星期几的问题。有时候,我们还想知道历史上某一天是星期几。通常,解决这个方法的有效办法是看日历,但是我们总不会随时随身带着日历,更不可能随时随身带着几千年的万年历。假如是想在计算机编程中计算某一天是星期几,预先把一本万年历存进去就更不现实了。这时候是不是有办法通过什么公式,从年月日推出这一天是星期几呢?



答案是肯定的。其实我们也常常在这样做。我们先举一个简单的例子。比如,知道了2004年5月1日是星期六,那么2004年5月31日“世界无烟日”是星期几就不难推算出来。我们可以掰着指头从1日数到31日,同时数星期,最后可以数出5月31日是星期一。其实运用数学计算,可以不用掰指头。我们知道星期是七天一轮回的,所以5月1日是星期六,七天之后的5月8日也是星期六。在日期上,8-1=7,正是7的倍数。同样,5月15日、5月22日和5月29日也是星期六,它们的日期和5月1日的差值分别是14、21和28,也都是7的倍数。那么5月31日呢?31-1=30,虽然不是7的倍数,但是31除以7,余数为2,这就是说,5月31日的星期,是在5月1日的星期之后两天。星期六之后两天正是星期一。



这个简单的计算告诉我们计算星期的一个基本思路:首先,先要知道在想算的日子之前的一个确定的日子是星期几,拿这一天做为推算的标准,也就是相当于一个计算的“原点”。其次,知道想算的日子和这个确定的日子之间相差多少天,用7除这个日期的差值,余数就表示想算的日子的星期在确定的日子的星期之后多少天。如果余数是0,就表示这两天的星期相同。显然,如果把这个作为“原点”的日子选为星期日,那么余数正好就等于星期几,这样计算就更方便了。



但是直接计算两天之间的天数,还是不免繁琐。比如1982年7月29日和2004年5月1日之间相隔7947天,就不是一下子能算出来的。它包括三段时间:一,1982年7月29日以后这一年的剩余天数;二,1983-2003这二十一个整年的全部天数;三,从2004年

元旦到5月1日经过的天数。第二段比较好算,它等于21*365+5=7670天,之所以要加5,是因为这段时间内有5个闰年。第一段和第三段就比较麻烦了,比如第三段,需要把5月之前的四个月的天数累加起来,再加上日期值,即31+29+31+30+1=122天。同理,第一段需要把7月之后的五个月的天数累加起来,再加上7月剩下的天数,一共是155天。所以总共的相隔天数是122+7670+155=7947天。



仔细想想,如果把“原点”日子的日期选为12月31日,那么第一段时间也就是一个整年,这样一来,第一段时间和第二段时间就可以合并计算,整年的总数正好相当于两个日子的年份差值减一。如果进一步把“原点”日子选为公元前1年12月31日(或者天文学家所使用的公元0年12月31日),这个整年的总数就正好是想算的日子的年份减一。这样简化之后,就只须计算两段时间:一,这么多整年的总天数;二,想算的日子是这一年的第几天。巧的是,按照公历的年月设置,这样反推回去,公元前1年12月31日正好是

星期日,也就是说,这样算出来的总天数除以7的余数正好是星期几。那么现在的问题就只有一个:这么多整年里面有多少闰年。这就需要了解公历的置闰规则了。



我们知道,公历的平年是365天,闰年是366天。置闰的方法是能被4整除的年份在2月加一天,但能被100整除的不闰,能被400整除的又闰。因此,像1600、2000、2400年都是闰年,而1700、1800、1900、2100年都是平年。公元前1年,按公历也是闰年。



因此,对于从公元前1年(或公元0年)12月31日到某一日子的年份Y之间的所有整年中的闰年数,就等于



[(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400],



[...]表示只取整数部分。第一项表示需要加上被4整除的年份数,第二项表示需要去掉被100整除的年份数,第三项表示需要再加上被400整除的年份数。之所以Y要减一,这样,我们就得到了第一个计算某一天是星期几的公式:



W = (Y-1)*365 + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400] + D.              (1)



其中D是这个日子在这一年中的累积天数。算出来的W就是公元前1年(或公元0年)12月31日到这一天之间的间隔日数。把W用7除,余数是几,这一天就是星期几。比如我们来算2004年5月1日:



W = (2004-1)*365 + [(2004-1)/4] - [(2004-1)/100] + [(2004-1)/400] +

   (31+29+31+30+1)

  = 731702,



731702 / 7 = 104528……6,余数为六,说明这一天是星期六。这和事实是符合的。



上面的公式(1)虽然很准确,但是计算出来的数字太大了,使用起来很不方便。仔细想想,其实这个间隔天数W的用数仅仅是为了得到它除以7之后的余数。这启发我们是不是可以简化这个W值,只要找一个和它余数相同的较小的数来代替,用数论上的术语

来说,就是找一个和它同余的较小的正整数,照样可以计算出准确的星期数。



显然,W这么大的原因是因为公式中的第一项(Y-1)*365太大了。其实,



(Y-1)*365 = (Y-1) * (364+1)

          = (Y-1) * (7*52+1)

          = 52 * (Y-1) * 7 + (Y-1),



这个结果的第一项是一个7的倍数,除以7余数为0,因此(Y-1)*365除以7的余数其实就等于Y-1除以7的余数。这个关系可以表示为:



(Y-1)*365 ≡ Y-1 (mod 7).



其中,≡是数论中表示同余的符号,mod 7的意思是指在用7作模数(也就是除数)的情况下≡号两边的数是同余的。因此,完全可以用(Y-1)代替(Y-1)*365,这样我们就得到了那个著名的、也是最常见到的计算星期几的公式:



W = (Y-1) + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400] + D.                  (2)



这个公式虽然好用多了,但还不是最好用的公式,因为累积天数D的计算也比较麻烦。是不是可以用月份数和日期直接计算呢?答案也是肯定的。我们不妨来观察一下各个月的日数,列表如下:



月  份:1月 2月  3月 4月 5月 6月 7月 8月 9月 10月 11月 12月

--------------------------------------------------------------------------

天  数: 31  28(29)  31   30   31   30   31   31   30   31    30    31



如果把这个天数都减去28(=4*7),不影响W除以7的余数值。这样我们就得到另一张表:



月  份:1月 2月 3月 4月 5月 6月 7月 8月 9月 10月 11月 12月

------------------------------------------------------------------------

剩余天数: 3   0(1)  3    2    3    2    3    3    2    3     2     3

平年累积: 3   3     6    8   11   13   16   19   21   24    26    29

闰年累积: 3   4     7    9   12   14   17   20   22   25    27    30



仔细观察的话,我们会发现除去1月和2月,3月到7月这五个月的剩余天数值是3,2,3,2,3;8月到12月这五个月的天数值也是3,2,3,2,3,正好是一个重复。相应的累积天数中,后一月的累积天数和前一月的累积天数之差减去28就是这个重复。正是因为这种规律的存在,平年和闰年的累积天数可以用数学公式很方便地表达:



    ╭ d;                 (当M=1)

D = {  31 + d;                 (当M=2)           (3)

    ╰ [ 13 * (M+1) / 5 ] - 7 + (M-1) * 28 + d + i.  (当M≥3)



其中[...]仍表示只取整数部分;M和d分别是想算的日子的月份和日数;平年i=0,闰年i=1。对于M≥3的表达式需要说明一下:[13*(M+1)/5]-7算出来的就是上面第二个表中的平年累积值,再加上(M-1)*28就是想算的日子的月份之前的所有月份的总天数。这是一个很巧妙的办法,利用取整运算来实现3,2,3,2,3的循环。比如,对2004年5月1日,有:



D = [ 13 * (5+1) / 5 ] - 7 + (5-1) * 28 + 1 + 1

  = 122,



这正是5月1日在2004年的累积天数。



假如,我们再变通一下,把1月和2月当成是上一年的“13月”和“14月”,不仅仍然符合这个公式,而且因为这样一来,闰日成了上一“年”(一共有14个月)的最后一天,成了d的一部分,于是平闰年的影响也去掉了,公式就简化成:



D = [ 13 * (M+1) / 5 ] - 7 + (M-1) * 28 + d.        (3≤M≤14)        (4)



上面计算星期几的公式,也就可以进一步简化成:



W = (Y-1) + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400] + [ 13 * (M+1) / 5 ] - 7

   + (M-1) * 28 + d.



因为其中的-7和(M-1)*28两项都可以被7整除,所以去掉这两项,W除以7的余数不变,公式变成:



W = (Y-1) + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400] + [ 13 * (M+1) / 5 ] + d.

(5)



当然,要注意1月和2月已经被当成了上一年的13月和14月,因此在计算1月和2月的日子的星期时,除了M要按13或14算,年份Y也要减一。比如,2004年1月1日是星期四,用这个公式来算,有:



W = (2003-1) + [(2003-1)/4] - [(2003-1)/100] + [(2003-1)/400] + [13*(13+1)/5]

   + 1

  = 2002 + 500 - 20 + 5 + 36 + 1

  = 2524;

2524 / 7 = 360……4.这和实际是一致的。



公式(5)已经是从年、月、日来算星期几的公式了,但它还不是最简练的,对于年份的处理还有改进的方法。我们先来用这个公式算出每个世纪第一年3月1日的星期,列表如下:



年份:  1(401,801,…,2001)                   101(501,901,…,2101)

--------------------------------------------------------------------

星期: 4                                      2

====================================================================

年份:201(601,1001,…,2201)                  301(701,1101,…,2301)

--------------------------------------------------------------------

星期:  0                                      5



可以看出,每隔四个世纪,这个星期就重复一次。假如我们把301(701,1101,…,2301)年3月1日的星期数看成是-2(按数论中对余数的定义,-2和5除以7的余数相同,所以可以做这样的变换),那么这个重复序列正好就是一个4,2,0,-2的等差数列。据此,我们可以得到下面的计算每个世纪第一年3月1日的星期的公式:



W = (4 - C mod 4) * 2 - 4.                                              (6)



式中,C是该世纪的世纪数减一,mod表示取模运算,即求余数。比如,对于2001年3月1日,C=20,则:



W = (4 - 20 mod 4) * 2 - 4

  = 8 - 4

  = 4.



把公式(6)代入公式(5),经过变换,可得:



(Y-1) + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400] ≡ (4 - C mod 4) * 2 - 1

(mod 7).                                                               (7)



因此,公式(5)中的(Y-1) + [(Y-1)/4] - [(Y-1)/100] + [(Y-1)/400]这四项,在计算每个世纪第一年的日期的星期时,可以用(4 - C mod 4) * 2 - 1来代替。这个公式写出来就是:



W = (4 - C mod 4) * 2 - 1 + [13 * (M+1) / 5] + d.                       (8)



有了计算每个世纪第一年的日期星期的公式,计算这个世纪其他各年的日期星期的公式就很容易得到了。因为在一个世纪里,末尾为00的年份是最后一年,因此就用不着再考虑“一百年不闰,四百年又闰”的规则,只须考虑“四年一闰”的规则。仿照由公式(1)

简化为公式(2)的方法,我们很容易就可以从式(8)得到一个比公式(5)更简单的计算任意

一天是星期几的公式:



W = (4 - C mod 4) * 2 - 1 + (y-1) + [y/4] + [13 * (M+1) / 5] + d.       (9)



式中,y是年份的后两位数字。



如果再考虑到取模运算不是四则运算,我们还可以把(4 - C mod 4) * 2进一步改写成只含四则运算的表达式。因为世纪数减一C除以4的商数q和余数r之间有如下关系:



4q + r = C,



其中r即是 C mod 4,因此,有:



r = C - 4q

  = C - 4 * [C/4].                                                     (10)



则



(4 - C mod 4) * 2 = (4 - C + 4 * [C/4]) * 2

                  = 8 - 2C + 8 * [C/4]

                  ≡ [C/4] - 2C + 1 (mod 7).                           (11)



把式(11)代入(9),得到:



W = [C/4] - 2C + y + [y/4] + [13 * (M+1) / 5] + d - 1.                 (12)



这个公式由世纪数减一、年份末两位、月份和日数即可算出W,再除以7,得到的余数是几就表示这一天是星期几,唯一需要变通的是要把1月和2月当成上一年的13月和14月,C和y都按上一年的年份取值。因此,人们普遍认为这是计算任意一天是星期几的最好的公式。这个公式最早是由德国数学家克里斯蒂安·蔡勒(Christian Zeller, 1822-1899)在1886年推导出的,因此通称为蔡勒公式(Zeller’s Formula)。为方便口算,式中的[13 * (M+1) / 5]也往往写成[26 * (M+1) / 10]。



现在仍然让我们来算2004年5月1日的星期,显然C=20,y=4,M=5,d=1,代入蔡勒公式,有:



W = [20/4] - 40 + 4 + 1 + [13 * (5+1) / 5] + 1 - 1

  = -15.



注意负数不能按习惯的余数的概念求余数,只能按数论中的余数的定义求余。为了方便计算,我们可以给它加上一个7的整数倍,使它变为一个正数,比如加上70,得到55。再除以7,余6,说明这一天是星期六。这和实际是一致的,也和公式(2)计算所得的结果一致。



最后需要说明的是,上面的公式都是基于公历(格里高利历)的置闰规则来考虑的。对于儒略历,蔡勒也推出了相应的公式是:



W = 5 - C + y + [y/4] + [13 * (M+1) / 5] + d - 1.                      (13)



这样,我们终于一劳永逸地解决了不查日历计算任何一天是星期几的问题。

2011年2月27日星期日

打印质数的各种算法

来源:原文链接


打印质数的算法应该是学习计算机编程的一个经典的问题,在这里想给大家展示一些方法,相信这些方法会对你的编程有一定的启发作用。请你注意几点,
  • 实际应用和教学应用有很大的差别。
  • 最后的那个使用编译时而不是运行时的方法大家可以重点看看。

教科书的示例

首先,先给一个教科书的示例。下面这个示例应该是教科书(至少是我上大学时的教科学)中算法复杂度最好的例子了。其想法很简单,先写一个判断是否是质数的函数isPrime(),然后从1到n分别调用isPrime()函数来检查。检查是否是质数的算法是核心,其简单的使用从2到n的开根的数作为除数。这样的算法复杂度几乎是O(n*log(n)),看上去不错,但其实很不经济。

#include <iostream>
using namespace std;
bool isPrime(int nr) {
for (int d = 2; (d * d) <= (nr + 1); ++d){
if (!(nr % d)){
return false;
}
}
return true;
}

int main (int argc, char * const argv[]) {
for (int i = 0; i < 50; ++i){
if (isPrime(i)){
cout << i << endl;
}
}
}


较好的算法


我们知道,我们的算法如果写成线性算法,也就是O(n),已经算是不错了,但是最好的是O(Log(n))的算法,这是一个级数级的算法,著名的二分取中(Binary Search)正是O(Log(n))的算法。通常来说,O(Log(n))的算法都是以排除法做为手段的。所以,找质数的算法完全可以采用排除法的方式。如下所示,这种算法的复杂度是O(n(log(logn)))。
示例:打印30以内的质数
一、初始化如下列表。
2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
二、把第一个数(2)取出来,去掉所有可以被2整除的数。
2  3     5     7     9    11    13    15    17    19    21    23    25    27    29
三、取第二个数(3),去掉所有可以被 3整除的数。
2  3     5     7          11    13          17    19          23    25          29
四、取第三个数(5),因为4已经被去除了,再去掉所有可以被5整除的数。
2  3     5     7          11    13          17    19          23                29
接下来的数是7,但是7的平方是49,其大于了30,所以我们可以停止计算了。剩下的数就是所有的质数了。

实际应用的算法

实际应用中,我们通常不会使用上述的两种算法,因为那是理论学院派的算法。实际中的算法是,我把质数事先就计算好,放在一个文件中,然后在程序启动时(注意是在启动时读这个文件,而不是运行时每调用一次就读一次文件),读取这个文件,然后打印出来就可以了。如果需要查找的化,二分查找或是hash表查找将会获得巨大的性能提升。当然,这样的方法对于空间来说比前面两个都要消耗得大,但是你可以有O(log(n))或是O(1)的时间复杂度。
所以,我想在这里提醒大家——实际和理论的的方法很不一样的,千万不要读书读成书呆子。在游戏编程的世界里,大量的数据都不是运行计算的,而都是写在文件中的。比如,一个火焰效果,一个人物跑动的动作,都是事先写在文件中的。

使用编译时而不是运行时

下面这个例子(本例参考于这里)你需要注意了,这是一个高级用法,使用模式来在编译时计算质数,而不是运行时。这种技术使用了C++编译器对模板的特化时的处理来生成自己相要的结果。这种方法在技术上是相当Cool的,但并不一定实用,这里只是想像大家展示这种用法。这是C++的最骨灰级的用法了。
请看下面的两个模板类,第一个模板以递归的方式检查是否是质数,第二个方法是递归的退出条件(当N=1时),对于模板的重载,请参看相关的C++书籍。
template<int N, int D = N - 1> 
struct isPrime {
enum {
result = (N % D) && isPrime<N, D-1>::result
};
};

template<int N>
struct isPrime<N, 1> {
enum {
result = true
};
};

于是,通过这个模板,我们可以使用下面的代码来检查是否是质数:
if (isPrime<3>::result)
cout << "Guess what: 3 is a prime!";

下一步,我们需要打出一个区间内的质数,所以,我们需要继续设计我们的print模板。
template<int N, bool ISPRIME> 
struct printIfPrime {
static inline void print() {}
};
template <int N>
struct printIfPrime<N, true> {
static inline void print() {
std::cout << N << endl;
}
};

从上面的代码中,我们可以看到,我们的第一个实际是什么也没做,而第二个有输出,注意第二个的模板参数中有一个true,其意味着那个质数的判断。于是我们就可以给出下面的代码来尝试着打印出一段区间内的质数:(请不要编译!!因为那会让编译器进入无限循环中,原因是printPrimes会不停地调用自己永不停止)

template<int N, int MAX>
struct printPrimes {
static inline void print() {
printIfPrime<N, isPrime<N>::result>::print();
printPrimes<N + 1, MAX>::print();
}
};

为了避免这个问题,你需要再加一个模板类,如下所示。这样当N变成MAX的时候,递归就结束了。

template<int N>
struct printPrimes<N, N> {
static inline void print() {
printIfPrime<N, isPrime<N>::result>::print();
}
};


最后,让我们来看看最终的调用:
int main (int argc, char * const argv[]) {
printPrimes<2, 40>::print();
return 0;
}

这个方法很NB,但是有两个问题:
  • 比较耗编译时间。
  • 不能在运行时输入MAX的值。
不过,相信这种玩法会启动你很多的编程思路。
(全文完)

2011年2月10日星期四

给定一个数组,找出不在数组中的最小的那个数字

来源:原文链接


这是在TL讨论中Liu xinyu给出的一个例子,觉得思路挺有启发的,所以整理记录一下。
给定一个数组,其内容是一些随机的、不重复的正整数,如:
{4, 23, 1, 8, 9, 21, 6, 12}
要求找出不在数组中出现的最小的那个数,比如这个数组中未在数组中出现的最小值是:2
这个问题实际应用的原型可以是一个ID分配系统,其使用一个数组来保存已分配的ID,每次回收就从数组中删除一个元素(O(n)),而分配则需要找到最小的那个可用的ID,就是这个算法要做的事情。
这个问题从naive的解法到快速的解法的思路转换是十分巧妙的,当然,如果之前没有接触过类似的题,注意到这个特性应该不是一件很容易的事。
设数组为A,大小为n,下标从1开始,下面是一系列逐步改进的算法:
一、穷举查找
一般的问题都可以通过这种很暴力的方式来做,从1到n逐个判断是否在数组中:
MIN-AVAILABLE-NUM(A, n)
for i = 1 to n
do if i not in A
   then return i
      return n+1
显然,这里的算法复杂度是O(n^2)
二、先排序再二分查找
第一种方法,每次查找都是线性查找,要改进最先想到的自然是二分查找,二分查找的前提是有序, 所以:
  1. 先排序,用O(nlgn)的快速排序、归并排序或者堆排序;因为数组中的元素是一些自然数,我们甚至可以使用O(n) 的基数排序,当然,需要更多的内存。
  2. 对1..n进行判断,复杂度也为O(nlgn)
所以,整体的算法复杂度为O(nlgn)
三、该数组的一个特性
其实仔细观察该数组A[1]..A[n],我们可以得出一个结论:如果该数组中存在未被使用的数,那么Max(A) > n。
证明很简单,假设Max(A) <= n,由于该数组大小为n,那么该数组中的元素只能是从1到n的某个排列,从而得出该数组中不存在未被使用的数,矛盾。
这个特性和抽屉原理有些类似之处。
从而我们可以有另外一个方法:
  1. 先排序
  2. 再利用该特性搜索
MIN-AVAILABLE-NUM(A, n)
for i = 1 to n
do A[i] > i
   then return i
      return n+1

注意到,如果我们使用基数排序,可以将复杂度降低到O(n)。
四、一个线性时间,线性空间的算法
第三个算法虽然能达到理论意义上的O(n),但是基数排序隐含的常数因子较大,而且不是原地排序,这里给出一个不需要排序的算法:
MIN-AVAILABLE-NUM(A, n)
for i = 1 to n
B[i] = 0
for i = 1 to n
do if A[i] < n
   then B[A[i]] = 1
      for  i = 1 to n
if(B[i] == 0) return i;
return n+1;
这里使用一个辅助数组B来表示1到n这些数是否存在在数组A中,只要不存在就将其标为0,最后在B中找到第一个值为0的便是我们要找的那个元素;如果B中元素全为1,这说明A使用了所有1到n这些数,那么返回的便是下一个n+1.
此处无须排序,且复杂度为O(n),但需要一个额外的O(n)的数组。
五、一个线性时间、常数空间的算法
利用快速排序的原理,我们可以在不使用额外数组的情况下达到O(n)的效率,原理为:
取1到n的中间值m = (1 + n)/2,用m将数组分成A1, A2两个部分,A1中的元素全部小于等于m,A2中的元素全部大于m(注意此处用的是下标,而不是A[m]),如果A1的大小为m,则空闲元素在A2中,这在前面证明过,然后就在A2中应用同样的方法。
MIN-AVAILABLE-NUM(A, low, up)
if(low == up) return low
m = (low + up) / 2
split = partition(A, low, up, m)
if a[split] == m
  then return MIN-AVAILABLE-NUM(A, low, split)
      else return MIN-AVAILABLE-NUM(A, split+1, up)
这里递归式为:T(n) = T(n/2) + O(n),根据主定理的第三种情况,复杂度为O(n),其实也就是一个等比数列:n + n/2 + n/4...