Miller-Rabin算法(转)(2007-3-22 23:42:00)

一.费马小定里
if n is prime and (a,n) equals one ,then a^(n-1) = 1 (mod n)

费马小定理只是个必要条件,符合费马小定理而非素数的数叫做Carmichael.

前3个Carmichael数是561,1105,1729。

Carmichael数是非常少的。

在1~100000000范围内的整数中，只有255个Carmichael数。

为此又有二次探测定理,以确保该数为素数:

二.二次探测定理

二次探测定理 如果p是一个素数，0<x<p,则方程x^2≡1(mod p)的解为x=1,p-1

根据以上两个定理,如到Miller-Rabin算法的一般步骤:

0、先计算出m、j，使得n-1=m*2^j，其中m是正奇数，j是非负整数

1、随机取一个b，2<=b

2、计算v=b^m mod n

3、如果v==1，通过测试，返回

4、令i=1

5、如果v=n-1，通过测试，返回

6、如果i==j，非素数，结束

7、v=v^2 mod n，i=i+1

8、循环到5

说明:

Miller-Rabin是随机算法

得到的结果的正确率为 75%,所以应该多次调用该函数,使正确概率提高为

1-(1/4)^p

//功能实现

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <functional>

typedef __int64 lltype;
const int MAX_COUNT = 6;//测试次数
lltype allFactor[65], nFactor;//质因分解结果(刚返回的时候是无序的)
void fFindFactor(const lltype &num);//主函数(num > 1)

lltype pollard_rho(const lltype &c, const lltype &num);
bool miller_rabin(const lltype &n, int count);//miller-rabin
bool fWitNess(const lltype &a, const lltype &n, char t, const lltype &u);

lltype modular_exp(const lltype &a, lltype b, const lltype &n);//d = a^b(mod n)
lltype fMultiModular(const lltype &a, lltype b, const lltype &n);//(a*b)%num
lltype fGCD(lltype a, lltype b);

void fFindFactor(const lltype &num)
{
    if ( miller_rabin(num, MAX_COUNT) == true )
    {
  allFactor[++nFactor] = num;
        return;
    }
    lltype factor;
    do
    {
        factor = pollard_rho(rand()%(num-1) + 1, num);
    }while ( factor >= num );
 //printf("%I64d ",factor);
    fFindFactor(factor);
    fFindFactor(num/factor);
}

lltype pollard_rho(const lltype &c, const lltype &num)
{
    int i(1), k(2);
    lltype x = rand() % num;
    lltype y = x, comDiv;
    do
    {
        ++i;
        if ( (x = fMultiModular(x, x, num) - c) < 0 )
   x += num;
        if ( x == y )
            break;
        comDiv = fGCD((y-x+num)%num, num);
        if ( comDiv > 1 && comDiv < num )
            return comDiv;
        if ( i == k )
        {
            y = x;
            k <<= 1;
        }
    }while ( true );
    return num;
}

//miller-rabin法测试素数, count 测试次数
bool miller_rabin(const lltype &n, int count)
{
    if ( n == 1 )
        return false;
    if ( n == 2 )
     return true;
    lltype a, u;    //为了谨慎起见用了 lltype, 而没有用 int
    char t(0);      //n-1 == (2^t)*u
    for (u = n-1; !(u & 0x1); u>>=1)
     ++t;
    for (int i = 1; i <= count; ++i)
    {
        a = rand() % (n-1) + 1;
        if ( fWitNess(a, n, t, u) )
     return false;
    }
    return true;
}

bool fWitNess(const lltype &a, const lltype &n, char t, const lltype &u)
{//if n must not be prime return true, else return false
 lltype x, y;
 x = modular_exp(a, u, n);
 while ( t-- )
 {
  y = fMultiModular(x, x, n);
  if ( y == 1 && x != 1  && x != n-1 )
   return true; //must not
  x = y;
 }
 return y != 1;
}

// d == a^b(mod n)
lltype modular_exp(const lltype &a, lltype b, const lltype &n)
{
    lltype d(1), dTemp(a % n);//当前二进制位表示的是进制数值
    while ( b > 0 )
    {
        if ( b & 0x1 )
            d = fMultiModular(d, dTemp, n);//不直接 d*dTemp 是怕溢出
        dTemp = fMultiModular(dTemp, dTemp, n);
        b >>= 1;
    }
    return d;
}

//back = (a*b)%num      ( n <= 2^62)
lltype fMultiModular(const lltype &a, lltype b, const lltype &n)
{
    lltype back(0), temp(a % n);
    //b %= n;
    while ( b > 0 )
    {
        if ( b & 0x1 )
        {
            if ( (back = back + temp) > n )
             back -= n;
  }
        if ( (temp <<= 1) > n )
         temp -= n;
        b >>= 1;
    }
    return back;
}

lltype fGCD(lltype a, lltype b)
{
 lltype c;
 while ( b )
 {
  c = b;
  b = a%b;
  a = c;
 }
 return a;
}

阅读全文(6395) | 评论:0 | 复制链接

(转)Prime Numbers, Factorization and Eule(2007-3-11 23:51:00)

Prime numbers and their properties were extensively studied by the ancient Greek mathematicians. Thousands of years later, we commonly use the different properties of integers that they discovered to solve problems. In this article we’ll review some definitions, well-known theorems, and number properties, and look at some problems associated with them.

A prime number is a positive integer, which is divisible on 1 and itself. The other integers, greater than 1, are composite. Coprime integers are a set of integers that have no common divisor other than 1 or -1.

The fundamental theorem of arithmetic:
Any positive integer can be divided in primes in essentially only one way. The phrase 'essentially one way' means that we do not consider the order of the factors important.

One is neither a prime nor composite number. One is not composite because it doesn’t have two distinct divisors. If one is prime, then number 6, for example, has two different representations as a product of prime numbers: 6 = 2 * 3 and 6 = 1 * 2 * 3. This would contradict the fundamental theorem of arithmetic.

Euclid’s theorem:
There is no largest prime number.

To prove this, let's consider only n prime numbers: p₁, p₂, …, p_n. But no prime p_i divides the number

N = p₁ * p₂ * … * p_n + 1,

so N cannot be composite. This contradicts the fact that the set of primes is finite.

Exercise 1. Sequence a_n is defined recursively:

Prove that a_i and a_j, i ¹ j are relatively prime.

Hint: Prove that a_n+1 = a₁a₂…a_n + 1 and use Euclid’s theorem.

Exercise 2. Ferma numbers F_n (n ≥ 0) are positive integers of the form

Prove that F_i and F_j, i ≠ j are relatively prime.

Hint: Prove that F_{n +1} = F₀F₁F₂…F_n + 2 and use Euclid’s theorem.

Dirichlet’s theorem about arithmetic progressions:
For any two positive coprime integers a and b there are infinitely many primes of the form a + n*b, where n > 0.

Trial division:
Trial division is the simplest of all factorization techniques. It represents a brute-force method, in which we are trying to divide n by every number i not greater than the square root of n. (Why don't we need to test values larger than the square root of n?) The procedure factor prints the factorization of number n. The factors will be printed in a line, separated with one space. The number n can contain no more than one factor, greater than n.

     void factor(int n) 
     { 
       int i; 
       for(i=2;i<=(int)sqrt(n);i++) 
       { 
         while(n % i == 0) 
         { 
           printf("%d ",i); 
           n /= i; 
         } 
       } 
       if (n > 1) printf("%d",n); 
       printf("\n"); 
     } 

Consider a problem that asks you to find the factorization of integer g(-2³¹ < g <2³¹) in the form

g = f₁ x f₂ x … x f_n or g = -1 x f₁ x f₂ x … x f_n

where f_i is a prime greater than 1 and f_i ≤ f_j for i < j.

For example, for g = -192 the answer is -192 = -1 x 2 x 2 x 2 x 2 x 2 x 2 x 3.

To solve the problem, it is enough to use trial division as shown in function factor.

Sieve of Eratosthenes:
The most efficient way to find all small primes was proposed by the Greek mathematician Eratosthenes. His idea was to make a list of positive integers not greater than n and sequentially strike out the multiples of primes less than or equal to the square root of n. After this procedure only primes are left in the list.

The procedure of finding prime numbers gen_primes will use an array primes[MAX] as a list of integers. The elements of this array will be filled so that

At the beginning we mark all numbers as prime. Then for each prime number i (i ≥ 2), not greater than √MAX, we mark all numbers i*i, i*(i + 1), … as composite.

     void gen_primes() 
     { 
       int i,j; 
       for(i=0;i<MAX;i++) primes[i] = 1; 
       for(i=2;i<=(int)sqrt(MAX);i++) 
         if (primes[i]) 
           for(j=i;j*i<MAX;j++) primes[i*j] = 0; 
     } 

For example, if MAX = 16, then after calling gen_primes, the array ‘primes’ will contain next values:

Goldbach's Conjecture:
For any integer n (n ≥ 4) there exist two prime numbers p₁ and p₂ such that p₁ + p₂ = n. In a problem we might need to find the number of essentially different pairs (p₁, p₂), satisfying the condition in the conjecture for a given even number n (4 ≤ n ≤ 2 15). (The word ‘essentially’ means that for each pair (p₁, p₂) we have p₁ ≤ p₂.)

For example, for n = 10 we have two such pairs: 10 = 5 + 5 and 10 = 3 + 7.

To solve this,as n ≤ 2¹⁵ = 32768, we’ll fill an array primes[32768] using function gen_primes. We are interested in primes, not greater than 32768.

The function FindSol(n) finds the number of different pairs (p₁, p₂), for which n = p₁ + p₂. As p₁ ≤ p₂, we have p₁ ≤ n/2. So to solve the problem we need to find the number of pairs (i, n – i), such that i and n – i are prime numbers and 2 ≤ i ≤ n/2.

     int FindSol(int n) 
     { 
       int i,res=0; 
       for(i=2;i<=n/2;i++) 
         if (primes[i] && primes[n-i]) res++; 
       return res; 
     } 

Euler’s totient function
The number of positive integers, not greater than n, and relatively prime with n, equals to Euler’s totient function φ (n). In symbols we can state that

φ (n) ={a Î N: 1 ≤ a ≤ n, gcd(a, n) = 1}

This function has the following properties:

1. If p is prime, then φ (p) = p – 1 and φ (p^a) = p a * (1 – 1/p) for any a.
2. If m and n are coprime, then φ (m * n) = φ (m) * φ (n).
3. If n = , then Euler function can be found using formula:

φ (n) = n * (1 – 1/p 1) * (1 – 1/p 2) * ... * (1 – 1/p k)

The function fi(n) finds the value of φ(n):

     int fi(int n) 
     { 
       int result = n; 
       for(int i=2;i*i <= n;i++) 
       { 
         if (n % i == 0) result -= result / i; 
         while (n % i == 0) n /= i; 
       } 
       if (n > 1) result -= result / n; 
       return result; 
     } 

For example, to find φ(616) we need to factorize the argument: 616 = 2³ * 7 * 11. Then, using the formula, we’ll get:

φ(616) = 616 * (1 – 1/2) * (1 – 1/7) * (1 – 1/11) = 616 * 1/2 * 6/7 * 10/11 = 240.

Say you've got a problem that, for a given integer n (0 < n ≤ 10⁹), asks you to find the number of positive integers less than n and relatively prime to n. For example, for n = 12 we have 4 such numbers: 1, 5, 7 and 11.

The solution: The number of positive integers less than n and relatively prime to n equals to φ(n). In this problem, then, we need do nothing more than to evaluate Euler’s totient function.

Or consider a scenario where you are asked to calculate a function Answer(x, y), with x and y both integers in the range [1, n], 1 ≤ n ≤ 50000. If you know Answer(x, y), then you can easily derive Answer(k*x, k*y) for any integer k. In this situation you want to know how many values of Answer(x, y) you need to precalculate. The function Answer is not symmetric.

For example, if n = 4, you need to precalculate 11 values: Answer(1, 1), Answer(1, 2), Answer(2, 1), Answer(1, 3), Answer(2, 3), Answer(3, 2), Answer(3, 1), Answer(1, 4), Answer(3, 4), Answer(4, 3) and Answer(4, 1).

The solution here is to let res(i) be the minimum number of Answer(x, y) to precalculate, where x, y Î{1, …, i}. It is obvious that res(1) = 1, because if n = 1, it is enough to know Answer(1, 1). Let we know res(i). So for n = i + 1 we need to find Answer(1, i + 1), Answer(2, i + 1), … , Answer(i + 1, i + 1), Answer(i + 1, 1), Answer(i + 1, 2), … , Answer(i + 1, i).

The values Answer(j, i + 1) and Answer(i + 1, j), j Î{1, …, i + 1}, can be found from known values if GCD(j, i + 1) > 1, i.e. if the numbers j and i + 1 are not common primes. So we must know all the values Answer(j, i + 1) and Answer(i + 1, j) for which j and i + 1 are coprime. The number of such values equals to 2 * φ (i + 1), where φ is an Euler’s totient function. So we have a recursion to solve a problem:

res(1) = 1,
res(i + 1) = res(i) + 2 * j (i + 1), i > 1

Euler’s totient theorem:
If n is a positive integer and a is coprime to n, then a ^{φ (n)} º 1 (mod n).

Fermat’s little theorem:
If p is a prime number, then for any integer a that is coprime to n, we have

a ^p ≡ a (mod p)

This theorem can also be stated as: If p is a prime number and a is coprime to p, then

a ^{p -1} ≡ 1 (mod p)

Fermat’s little theorem is a special case of Euler’s totient theorem when n is prime.

The number of divisors:
If n = , then the number of its positive divisors equals to

(a₁ + 1) * (a₂ + 1) * … * (a_k + 1)

For a proof, let A i be the set of divisors . Any divisor of number n can be represented as a product x₁ * x₂ * … * x _k , where x_i Î A_i. As |A_i| = a_i + 1, we have

(a₁ + 1) * (a₂ + 1) * … * (a_k + 1)

possibilities to get different products x₁ * x₂ * … * x_k.

For example, to find the number of divisors for 36, we need to factorize it first: 36 = 2² * 3². Using the formula above, we’ll get the divisors amount for 36. It equals to (2 + 1) * (2 + 1) = 3 * 3 = 9. There are 9 divisors for 36: 1, 2, 3, 4, 6, 9, 12, 18 and 36.

Here's another problem to think about: For a given positive integer n (0 < n < 2³¹) we need to find the number of such m that 1 ≤ m ≤ n, GCD(m, n) ≠ 1 and GCD(m, n) ≠ m. For example, for n = 6 we have only one such number m = 4.

The solution is to subtract from n the amount of numbers, coprime with it (its amount equals to φ(n)) and the amount of its divisors. But the number 1 simultaneously is coprime with n and is a divisor of n. So to obtain the difference we must add 1. If n = is a factorization of n, the number n has (a₁ + 1) * (a₂ + 1) * … * (a_k + 1) divisors. So the answer to the problem for a given n equals to

n – φ(n) – (a₁ + 1) * (a₂ + 1) * … * (a_k + 1) + 1

转自：http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=primeNumbers

阅读全文(1997) | 评论:0 | 复制链接

类设计者的核查表(2007-3-9 0:10:00)

摘自<<c++沉思录>>

1.你的类需要一个构造函数吗？

2.你的数据成员是私有的吗？

3.你的类需要一个无参的构造函数吗？

4.是不是每个构造函数初始化所有的数据成员？

5.类需要析构函数吗？

6.类需要一个虚析构函数吗？

7.你的类需要复制构造函数吗？

8.你的类需要一个赋值操作符吗？

9.你的赋值操作符能正确地将对象赋给对象本身吗？

10.你的类需要定义关系操作符吗？

11.删除数组时你记住了用delete[]吗？

12.记得在复制构造函数和赋值操作符的参数类型中加上const

了吗？

13.如果函数有引用参数，它们应该是const引用吗？

14.记得适当地声明成员函数为const的了吗？

阅读全文(1895) | 评论:0 | 复制链接

杂类集(2007-2-3 21:17:00)

 Great Common Divisor
最大公约数

n1.Euclid算法

 

重复gcd(a,b) = gcd(a,b mod a)至b mod a等于0 ，

gcd(b,0)=b, b最后的取值为a,b初值的最大公约数

 

2.stein算法

1. 

1.如果A=0，B是最大公约数，算法结束

2.   如果B=0，A是最大公约数，算法结束

3.   设置A1 = A、B1=B和C1 = 1

4.    如果An和Bn都是偶数，则An+1 =An /2，Bn+1 =Bn /2，Cn+1 =Cn *2(乘2整数左移一位，除2整数右移一位即可)

5.

5.如果An是偶数，Bn不是偶数，则An+1 =An /2，Bn+1 =Bn ，Cn+1 =Cn (显然，2不可能为任何奇数的约数)

6. 

6.如果Bn是偶数，An不是偶数，则Bn+1 =Bn /2，An+1 =An ，Cn+1 =Cn (同上J)

7.    如果An和Bn都不是偶数，则An+1 =|An -Bn|，Bn+1 =min(An,Bn)，Cn+1 =Cn

8.    n++，转4

注：、Cn为两数最大公约数

 

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

给定两个任意正整数n和任一质数p,求将n!表示成算术基本定理的形式的时候，p的指数．

count=0;

while(n)

{

     count+=n/p;

     n/=p;

}

cout<<count<<endl;

阅读全文(2232) | 评论:0 | 复制链接

hash function(2006-9-29 22:02:00)

经典字符串中出现的hash函数

（转载自http://www.ccw.com.cn/htm/app/aprog/01_8_22_3.asp）

1．  PHP中出现的字符串Hash函数

static unsigned long str_hash_function (char *arKey , unsigned int nKeyLength) {

int h=0;

int g;

char *arEnd=arKey+nKeyLength;

while (arKey<arEnd) {

    h=(h<<4)+*arKey++;

    if ((g=(h&0xF0000000))) {

        h=h^(g>>24);

        h=h^g;

}

}

return h;

}

2．  OpenSSL中出现的字符串Hash函数

unsigned long str_hash_function (char *str) {

int i , l;

unsigned long ret=0;

unsigned short *s;

if (str==NULL) return (0);

l=(strlen(str)+1)/2;

s=(unsigned short *) str;

for (i=0;i<l;i++) {

ret^=(s[i]<<(i&oxof));

return ret;

}

unsinged long str_hash_function2 (const char *str) {

    unsigned long ret=0;

    long n;

    unsigned long v;

    int r;

    if ((c==NULL)||(*c==’\0’)) return ret;

    n=0x100;

            while (*c) {

               v=n| (*c);

               n+=0x100;

               r=(int )((v>>2)^v)&0x0f;

               ret=(ret<>(32-r));

               ret&=0xFFFFFFFFL;

               retu^=v*v;

               c++;

            }

return ((ret>>16)^ret);

}

3．  MySql中出现的字符串Hash函数(完全看不懂就不写了)

4．  很简单又好用的一个字符串的Hash函数

int  hash_function () {

            int h=0;

            int i;

            for (i=0;b[i];i++) h=5*h+b[i];

            return h;

}

阅读全文(3224) | 评论:0 | 复制链接

计算几何算法概览(转)(2006-9-6 21:28:00)

怒火之袍
计算几何算法概览 一、引言 　　计算机的出现使得很多原本十分繁琐的工作得以大幅度简化，但是也有一些在人们直观看来很容易的问题却需要拿出一套并不简单的通用解决方案，比如几何问题。作为计算机科学的一个分支，计算几何主要研究解决几何问题的算法。在现代工程和数学领域，计算几何在图形学、机器人技术、超大规模集成电路设计和统计等诸多领域有着十分重要的应用。在本文中，我们将对计算几何常用的基本算法做一个全面的介绍，希望对您了解并应用计算几何的知识解决问题起到帮助。 二、目录 　　本文整理的计算几何基本概念和常用算法包括如下内容： 　　矢量的概念 　　矢量加减法 　　矢量叉积 　　折线段的拐向判断 　　判断点是否在线段上 　　判断两线段是否相交 　　判断线段和直线是否相交 　　判断矩形是否包含点 　　判断线段、折线、多边形是否在矩形中 　　判断矩形是否在矩形中 　　判断圆是否在矩形中 　　判断点是否在多边形中 　　判断线段是否在多边形内 　　判断折线是否在多边形内 　　判断多边形是否在多边形内 　　判断矩形是否在多边形内 　　判断圆是否在多边形内 　　判断点是否在圆内 　　判断线段、折线、矩形、多边形是否在圆内 　　判断圆是否在圆内 　　计算点到线段的最近点 　　计算点到折线、矩形、多边形的最近点 　　计算点到圆的最近距离及交点坐标 　　计算两条共线的线段的交点 　　计算线段或直线与线段的交点 　　求线段或直线与折线、矩形、多边形的交点 　　求线段或直线与圆的交点 　　凸包的概念 　　凸包的求法 三、算法介绍 　　矢量的概念： 　　如果一条线段的端点是有次序之分的，我们把这种线段成为有向线段(directed segment)。如果有向线段p1p2的起点p1在坐标原点，我们可以把它称为矢量(vector)p2。 　　矢量加减法： 　　设二维矢量P = ( x1, y1 )，Q = ( x2 , y2 )，则矢量加法定义为： P + Q = ( x1 + x2 , y1 + y2 )，同样的，矢量减法定义为： P - Q = ( x1 - x2 , y1 - y2 )。显然有性质 P + Q = Q + P，P - Q = - ( Q - P )。 　　矢量叉积： 　　计算矢量叉积是与直线和线段相关算法的核心部分。设矢量P = ( x1, y1 )，Q = ( x2, y2 )，则矢量叉积定义为由(0,0)、p1、p2和p1+p2所组成的平行四边形的带符号的面积，即：P × Q = x1*y2 - x2*y1，其结果是一个标量。显然有性质 P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q )。一般在不加说明的情况下，本文下述算法中所有的点都看作矢量，两点的加减法就是矢量相加减，而点的乘法则看作矢量叉积。 　　叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系： 　　若 P × Q > 0 , 则P在Q的顺时针方向。 　　若 P × Q < 0 , 则P在Q的逆时针方向。 　　若 P × Q = 0 , 则P与Q共线，但可能同向也可能反向。 　　折线段的拐向判断： 　　折线段的拐向判断方法可以直接由矢量叉积的性质推出。对于有公共端点的线段p0p1和p1p2，通过计算(p2 - p0) × (p1 - p0)的符号便可以确定折线段的拐向： 　　若(p2 - p0) × (p1 - p0) > 0,则p0p1在p1点拐向右侧后得到p1p2。 　　若(p2 - p0) × (p1 - p0) < 0,则p0p1在p1点拐向左侧后得到p1p2。 　　若(p2 - p0) × (p1 - p0) = 0,则p0、p1、p2三点共线。 　　具体情况可参照下图：     　　判断点是否在线段上： 　　设点为Q，线段为P1P2 ，判断点Q在该线段上的依据是：( Q - P1 ) × ( P2 - P1 ) = 0 且 Q 在以 P1，P2为对角顶点的矩形内。前者保证Q点在直线P1P2上，后者是保证Q点不在线段P1P2的延长线或反向延长线上，对于这一步骤的判断可以用以下过程实现： 　　ON-SEGMENT(pi,pj,pk) 　　if min(xi,xj) <= xk <= max(xi,xj) and min(yi,yj) <= yk <= max(yi,yj) 　　then return true; 　　else return false; 　　特别要注意的是，由于需要考虑水平线段和垂直线段两种特殊情况，min(xi,xj)<=xk<=max(xi,xj)和min(yi,yj)<=yk<=max(yi,yj)两个条件必须同时满足才能返回真值。 　　判断两线段是否相交： 　　我们分两步确定两条线段是否相交： 　　(1)快速排斥试验 　　　　设以线段 P1P2 为对角线的矩形为R， 设以线段 Q1Q2 为对角线的矩形为T，如果R和T不相交，显然两线段不会相交。 　　(2)跨立试验 　　　　如果两线段相交，则两线段必然相互跨立对方。若P1P2跨立Q1Q2 ，则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧，即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0。上式可改写成( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0。当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时，说明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共线，但是因为已经通过快速排斥试验，所以 P1 一定在线段 Q1Q2上；同理，( Q2 - Q1 ) ×(P2 - Q1 ) = 0 说明 P2 一定在线段 Q1Q2上。所以判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。同理判断Q1Q2跨立P1P2的依据是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。具体情况如下图所示：     　　在相同的原理下，对此算法的具体的实现细节可能会与此有所不同，除了这种过程外，大家也可以参考《算法导论》上的实现。 　　判断线段和直线是否相交： 　　有了上面的基础，这个算法就很容易了。如果线段P1P2和直线Q1Q2相交，则P1P2跨立Q1Q2，即：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。 　　判断矩形是否包含点： 　　只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。      　　判断线段、折线、多边形是否在矩形中： 　　因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。 　　判断矩形是否在矩形中： 　　只要比较左右边界和上下边界就可以了。 　　判断圆是否在矩形中： 　　很容易证明，圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。 　　判断点是否在多边形中： 　　判断点P是否在多边形中是计算几何中一个非常基本但是十分重要的算法。以点P为端点，向左方作射线L，由于多边形是有界的，所以射线L的左端一定在多边形外，考虑沿着L从无穷远处开始自左向右移动，遇到和多边形的第一个交点的时候，进入到了多边形的内部，遇到第二个交点的时候，离开了多边形，……所以很容易看出当L和多边形的交点数目C是奇数的时候，P在多边形内，是偶数的话P在多边形外。 　　但是有些特殊情况要加以考虑。如图下图(a)(b)(c)(d)所示。在图(a)中，L和多边形的顶点相交，这时候交点只能计算一个；在图(b)中，L和多边形顶点的交点不应被计算；在图(c)和(d) 中，L和多边形的一条边重合，这条边应该被忽略不计。如果L和多边形的一条边重合，这条边应该被忽略不计。      　　为了统一起见，我们在计算射线L和多边形的交点的时候，1。对于多边形的水平边不作考虑；2。对于多边形的顶点和L相交的情况，如果该顶点是其所属的边上纵坐标较大的顶点，则计数，否则忽略；3。对于P在多边形边上的情形，直接可判断P属于多边行。由此得出算法的伪代码如下：     count ← 0;     以P为端点，作从右向左的射线L;      for 多边形的每条边s      do if P在边s上            then return true;         if s不是水平的           then if s的一个端点在L上                  if 该端点是s两端点中纵坐标较大的端点                    then count ← count+1                else if s和L相交                  then count ← count+1;     if count mod 2 = 1        then return true;     else return false; 　　其中做射线L的方法是：设P'的纵坐标和P相同，横坐标为正无穷大（很大的一个正数），则P和P'就确定了射线L。 　　判断点是否在多边形中的这个算法的时间复杂度为O(n)。 　　另外还有一种算法是用带符号的三角形面积之和与多边形面积进行比较，这种算法由于使用浮点数运算所以会带来一定误差，不推荐大家使用。 　　判断线段是否在多边形内： 　　线段在多边形内的一个必要条件是线段的两个端点都在多边形内，但由于多边形可能为凹，所以这不能成为判断的充分条件。如果线段和多边形的某条边内交（两线段内交是指两线段相交且交点不在两线段的端点），因为多边形的边的左右两侧分属多边形内外不同部分，所以线段一定会有一部分在多边形外(见图a)。于是我们得到线段在多边形内的第二个必要条件：线段和多边形的所有边都不内交。 　　线段和多边形交于线段的两端点并不会影响线段是否在多边形内；但是如果多边形的某个顶点和线段相交，还必须判断两相邻交点之间的线段是否包含于多边形内部（反例见图b)。     　　因此我们可以先求出所有和线段相交的多边形的顶点，然后按照X-Y坐标排序(X坐标小的排在前面，对于X坐标相同的点，Y坐标小的排在前面，这种排序准则也是为了保证水平和垂直情况的判断正确)，这样相邻的两个点就是在线段上相邻的两交点，如果任意相邻两点的中点也在多边形内，则该线段一定在多边形内。 　　证明如下： 　 　　命题1： 　　　　如果线段和多边形的两相邻交点P1 ，P2的中点P' 也在多边形内，则P1, P2之间的所有点都在多边形内。      　　证明： 　　　　假设P1,P2之间含有不在多边形内的点，不妨设该点为Q，在P1, P'之间，因为多边形是闭合曲线，所以其内外部之间有界，而P1属于多边行内部，Q属于多边性外部，P'属于多边性内部，P1-Q-P'完全连续，所以P1Q和QP'一定跨越多边形的边界，因此在P1,P'之间至少还有两个该线段和多边形的交点，这和P1P2是相邻两交点矛盾，故命题成立。证毕。 　　由命题1直接可得出推论： 　　推论2： 　　　　设多边形和线段PQ的交点依次为P1,P2,……Pn，其中Pi和Pi+1是相邻两交点，线段PQ在多边形内的充要条件是：P，Q在多边形内且对于i =1, 2,……, n-1，Pi ,Pi+1的中点也在多边形内。 　　在实际编程中，没有必要计算所有的交点，首先应判断线段和多边形的边是否内交，倘若线段和多边形的某条边内交则线段一定在多边形外；如果线段和多边形的每一条边都不内交，则线段和多边形的交点一定是线段的端点或者多边形的顶点，只要判断点是否在线段上就可以了。 　　至此我们得出算法如下：     if 线端PQ的端点不都在多边形内        then return false;     点集pointSet初始化为空;     for 多边形的每条边s       do if 线段的某个端点在s上            then 将该端点加入pointSet;          else if s的某个端点在线段PQ上            then 将该端点加入pointSet;          else if s和线段PQ相交 // 这时候已经可以肯定是内交了            then return false;     将pointSet中的点按照X-Y坐标排序;     for pointSet中每两个相邻点 pointSet[i] , pointSet[ i+1]       do if pointSet[i] , pointSet[ i+1] 的中点不在多边形中            then return false;     return true; 　　这个过程中的排序因为交点数目肯定远小于多边形的顶点数目n，所以最多是常数级的复杂度，几乎可以忽略不计。因此算法的时间复杂度也是O(n)。 　　判断折线是否在多边形内： 　　只要判断折线的每条线段是否都在多边形内即可。设折线有m条线段，多边形有n个顶点，则该算法的时间复杂度为O(m*n)。 　　判断多边形是否在多边形内： 　　只要判断多边形的每条边是否都在多边形内即可。判断一个有m个顶点的多边形是否在一个有n个顶点的多边形内复杂度为O(m*n)。 　　判断矩形是否在多边形内： 　　将矩形转化为多边形，然后再判断是否在多边形内。 　　判断圆是否在多边形内： 　　只要计算圆心到多边形的每条边的最短距离，如果该距离大于等于圆半径则该圆在多边形内。计算圆心到多边形每条边最短距离的算法在后文阐述。 　　判断点是否在圆内： 　　计算圆心到该点的距离，如果小于等于半径则该点在圆内。 　　判断线段、折线、矩形、多边形是否在圆内： 　　因为圆是凸集，所以只要判断是否每个顶点都在圆内即可。 　　判断圆是否在圆内： 　　设两圆为O1,O2，半径分别为r1, r2，要判断O2是否在O1内。先比较r1，r2的大小，如果r1<r2则O2不可能在O1内；否则如果两圆心的距离大于r1 - r2 ，则O2不在O1内；否则O2在O1内。 　　计算点到线段的最近点： 　　如果该线段平行于X轴（Y轴），则过点point作该线段所在直线的垂线，垂足很容易求得，然后计算出垂足，如果垂足在线段上则返回垂足，否则返回离垂足近的端点；如果该线段不平行于X轴也不平行于Y轴，则斜率存在且不为0。设线段的两端点为pt1和pt2，斜率为：k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );该直线方程为：y = k* ( x - pt1.x) + pt1.y。其垂线的斜率为 - 1 / k，垂线方程为：y = (-1/k) * (x - point.x) + point.y 。 　　联立两直线方程解得：x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1) ，y = k * ( x - pt1.x) + pt1.y;然后再判断垂足是否在线段上，如果在线段上则返回垂足；如果不在则计算两端点到垂足的距离，选择距离垂足较近的端点返回。 　　计算点到折线、矩形、多边形的最近点： 　　只要分别计算点到每条线段的最近点，记录最近距离，取其中最近距离最小的点即可。 　　计算点到圆的最近距离及交点坐标： 　　如果该点在圆心，因为圆心到圆周任一点的距离相等，返回UNDEFINED。 　　连接点P和圆心O，如果PO平行于X轴，则根据P在O的左边还是右边计算出最近点的横坐标为centerPoint.x - radius 或 centerPoint.x + radius。如果PO平行于Y轴，则根据P在O的上边还是下边计算出最近点的纵坐标为 centerPoint.y -+radius或 centerPoint.y - radius。如果PO不平行于X轴和Y轴，则PO的斜率存在且不为0，这时直线PO斜率为k = （ P.y - O.y ）/ ( P.x - O.x )。直线PO的方程为：y = k * ( x - P.x) + P.y。设圆方程为:(x - O.x ) ^2 + ( y - O.y ) ^2 = r ^2，联立两方程组可以解出直线PO和圆的交点，取其中离P点较近的交点即可。 　　计算两条共线的线段的交点： 　　对于两条共线的线段，它们之间的位置关系有下图所示的几种情况。图(a)中两条线段没有交点；图 (b) 和 (d) 中两条线段有无穷焦点；图 (c) 中两条线段有一个交点。设line1是两条线段中较长的一条，line2是较短的一条，如果line1包含了line2的两个端点，则是图(d)的情况，两线段有无穷交点；如果line1只包含line2的一个端点，那么如果line1的某个端点等于被line1包含的line2的那个端点，则是图(c)的情况，这时两线段只有一个交点，否则就是图(b)的情况，两线段也是有无穷的交点；如果line1不包含line2的任何端点，则是图(a)的情况，这时两线段没有交点。 　　计算线段或直线与线段的交点: 　　设一条线段为L0 = P1P2，另一条线段或直线为L1 = Q1Q2 ，要计算的就是L0和L1的交点。 　1． 首先判断L0和L1是否相交（方法已在前文讨论过），如果不相交则没有交点，否则说明L0和L1一定有交点，下面就将L0和L1都看作直线来考虑。 　2． 如果P1和P2横坐标相同，即L0平行于Y轴 　　a) 若L1也平行于Y轴， 　　　　i. 若P1的纵坐标和Q1的纵坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； 　　　　ii. 否则说明L0和L1平行，他们没有交点； 　　b) 若L1不平行于Y轴，则交点横坐标为P1的横坐标，代入到L1的直线方程中可以计算出交点纵坐标； 　3． 如果P1和P2横坐标不同，但是Q1和Q2横坐标相同，即L1平行于Y轴，则交点横坐标为Q1的横坐标，代入到L0的直线方程中可以计算出交点纵坐标； 　4． 如果P1和P2纵坐标相同，即L0平行于X轴 　　a) 若L1也平行于X轴， 　　　　i. 若P1的横坐标和Q1的横坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； 　　　　ii. 否则说明L0和L1平行，他们没有交点； 　　b) 若L1不平行于X轴，则交点纵坐标为P1的纵坐标，代入到L1的直线方程中可以计算出交点横坐标； 　5． 如果P1和P2纵坐标不同，但是Q1和Q2纵坐标相同，即L1平行于X轴，则交点纵坐标为Q1的纵坐标，代入到L0的直线方程中可以计算出交点横坐标； 　6． 剩下的情况就是L1和L0的斜率均存在且不为0的情况 　　a) 计算出L0的斜率K0，L1的斜率K1 ； 　　b) 如果K1 = K2  　　　　i. 如果Q1在L0上，则说明L0和L1共线，假如L1是直线的话有无穷交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； 　　　　ii. 如果Q1不在L0上，则说明L0和L1平行，他们没有交点。 　　c) 联立两直线的方程组可以解出交点来 　　这个算法并不复杂，但是要分情况讨论清楚，尤其是当两条线段共线的情况需要单独考虑，所以在前文将求两条共线线段的算法单独写出来。另外，一开始就先利用矢量叉乘判断线段与线段（或直线）是否相交，如果结果是相交，那么在后面就可以将线段全部看作直线来考虑。需要注意的是，我们可以将直线或线段方程改写为ax+by+c=0的形式，这样一来上述过程的部分步骤可以合并，缩短了代码长度，但是由于先要求出参数，这种算法将花费更多的时间。 　　求线段或直线与折线、矩形、多边形的交点： 　　分别求与每条边的交点即可。 　　求线段或直线与圆的交点: 　　设圆心为O，圆半径为r，直线（或线段）L上的两点为P1,P2。 　　1. 如果L是线段且P1，P2都包含在圆O内，则没有交点；否则进行下一步。 　　2. 如果L平行于Y轴， 　　　a) 计算圆心到L的距离dis； 　　　b) 如果dis > r 则L和圆没有交点； 　　　c) 利用勾股定理，可以求出两交点坐标，但要注意考虑L和圆的相切情况。 　　3. 如果L平行于X轴，做法与L平行于Y轴的情况类似； 　　4. 如果L既不平行X轴也不平行Y轴，可以求出L的斜率K，然后列出L的点斜式方程，和圆方程联立即可求解出L和圆的两个交点； 　　5. 如果L是线段，对于2，3，4中求出的交点还要分别判断是否属于该线段的范围内。 　　凸包的概念： 　　点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。下图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。     　 　　凸包的求法： 　　现在已经证明了凸包算法的时间复杂度下界是O(n*logn),但是当凸包的顶点数h也被考虑进去的话，Krikpatrick和Seidel的剪枝搜索算法可以达到O(n*logh)，在渐进意义下达到最优。最常用的凸包算法是Graham扫描法和Jarvis步进法。本文只简单介绍一下Graham扫描法，其正确性的证明和Jarvis步进法的过程大家可以参考《算法导论》。 　　对于一个有三个或以上点的点集Q，Graham扫描法的过程如下： 　　令p0为Q中Y-X坐标排序下最小的点  　　设<p1,p2,...pm>为对其余点按以p0为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距p0最远的点外全部移除 　　压p0进栈S 　　压p1进栈S 　　压p2进栈S     for i ← 3 to m       do while 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧         对S弹栈       压pi进栈S     return S; 　 　　此过程执行后，栈S由底至顶的元素就是Q的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。 四、结语 　　尽管人类对几何学的研究从古代起便没有中断过，但是具体到借助计算机来解决几何问题的研究，还只是停留在一个初级阶段，无论从应用领域还是发展前景来看，计算几何学都值得我们认真学习、加以运用，希望这篇文章能带你走进这个丰富多彩的世界。

阅读全文(3044) | 评论:0 | 复制链接

Bellman-Ford算法(2006-9-5 22:16:00)

该算法用来求边上权值为任意值的单源最短路径问题

伪代码如下：

void bellman_ford(int v)

{

   for 1 to n

   initialize dist[i][v];   //此时只经过一条边

  for 2 to n-1  (i)           //经过的边数不大于i条时

     for  1 to n    (j)          //对于每一个目标顶点

        for 1 to n     (k)     //经过该顶点时，与当前最小值比较，并更新当前值

          if edge[k][j] > 0 && dist[k] > edge[k][j]+dist[j]

            更新当前值

}//

阅读全文(7439) | 评论:3 | 复制链接

LIS(2006-8-20 21:04:00)

#include <iostream>
using namespace std;

#define MAXNUM 1000

char c[MAXNUM];
char b[MAXNUM+1];

void solve()
{
     int i,first,last,mid,len;
    
     len = 1;
     i = 1;
     b[1] = c[0];
    
     while (c[i])
     {
         first = 1;
         last = len+1;
        
         while (first < last)
         {
            mid = (first+last)>>1;
            if (b[mid] < c[i])first = mid+1;
            else last = mid;
         }
         b[first] = c[i++];
         if (len < first)len = first;
     }
     cout << len << endl;
}

int main()
{
    cin >> c;
    solve();
    system("pause");
    return 0;
}

ps:这是在看了别人的博客后自己试着去写的一个关于求最长递增子序列的程序。时间复杂度为O(nlogn).

阅读全文(2685) | 评论:0 | 复制链接

Josephus问题的数学方法(转)(2006-8-20 11:02:00)

    无论是用链表实现还是用数组实现都有一个共同点：要模拟整个游戏过程，不仅程序写起来比较烦，而且时间复杂度高达O(nm)，当n，m非常大(例如上百万，上千万)的时候，几乎是没有办法在短时间内出结果的。我们注意到原问题仅仅是要求出最后的胜利者的序号，而不是要读者模拟整个过程。因此如果要追求效率，就要打破常规，实施一点数学策略。
为了讨论方便，先把问题稍微改变一下，并不影响原意：

问题描述：n个人（编号0~(n-1))，从0开始报数，报到(m-1)的退出，剩下的人继续从0开始报数。求胜利者的编号。

我们知道第一个人(编号一定是m%n-1) 出列之后，剩下的n-1个人组成了一个新的约瑟夫环（以编号为k=m%n的人开始）:
  k  k+1  k+2  ... n-2, n-1, 0, 1, 2, ... k-2
并且从k开始报0。

现在我们把他们的编号做一下转换：
k     --> 0
k+1   --> 1
k+2   --> 2
...
...
k-2   --> n-2
k-1   --> n-1

变换后就完完全全成为了(n-1)个人报数的子问题，假如我们知道这个子问题的解：例如x是最终的胜利者，那么根据上面这个表把这个x变回去不刚好就是n个人情况的解吗？！！变回去的公式很简单，相信大家都可以推出来：x'=(x+k)%n

如何知道(n-1)个人报数的问题的解？对，只要知道(n-2)个人的解就行了。(n-2)个人的解呢？当然是先求(n-3)的情况 ---- 这显然就是一个倒推问题！好了，思路出来了，下面写递推公式：

令f[i]表示i个人玩游戏报m退出最后胜利者的编号，最后的结果自然是f[n]

递推公式
f[1]=0;
f[i]=(f[i-1]+m)%i;  (i>1)

有了这个公式，我们要做的就是从1-n顺序算出f[i]的数值，最后结果是f[n]。因为实际生活中编号总是从1开始，我们输出f[n]+1

由于是逐级递推，不需要保存每个f[i]，程序也是异常简单：

＃i nclude <stdio.h>

main()
{
  int n, m, i, s=0;
  printf ("N M = "); scanf("%d%d", &n, &m);
  for (i=2; i<=n; i++) s=(s+m)%i;
  printf ("The winner is %d\n", s+1);
}

这个算法的时间复杂度为O(n)，相对于模拟算法已经有了很大的提高。算n，m等于一百万，一千万的情况不是问题了。可见，适当地运用数学策略，不仅可以让编程变得简单，而且往往会成倍地提高算法执行效率。

阅读全文(3359) | 评论:2 | 复制链接

费马最後定理（转）(2006-8-19 20:13:00)

被公认执世界报纸牛耳地位地位的纽约时报於1993年6月24日在其一版头题刊登了一则有关数学难题得以解决的消息，那则消息的标题是「在陈年数学困局中，终於有人呼叫『我找到了』」。时报一版的开始文章中还附了一张留着长发、穿着中古世纪欧洲学袍的男人照片。这个古意盎然的男人，就是法国的数学家费马（Pierre de Fermat）（费马小传请参考附录）。费马是十七世纪最卓越的数学家之一，他在数学许多领域中都有极大的贡献，因为他的本行是专业的律师，为了表彰他的数学造诣，世人冠以「业余王子」之美称，在三百六十多年前的某一天，费马正在阅读一本古希腊数学家戴奥芬多斯的数学书时，突然心血来潮在书页的空白处，写下一个看起来很简单的定理这个定理的内容是有关一个方程式 x² + y² =z²的正整数解的问题，当n=2时就是我们所熟知的毕氏定理（中国古代又称勾股弦定理）：x² + y² =z²，此处z表一直角形之斜边而x、y为其之两股，也就是一个直角三角形之斜边的平方等於它的两股的平方和，这个方程式当然有整数解（其实有很多），例如：x=3、y=4、z=5；x=6、y=8、z=10；x=5、y=12、z=13…等等。

    费马声称当n>2时，就找不到满足xⁿ +yⁿ = zⁿ的整数解，例如：方程式x³ +y³=z³就无法找到整数解。

    当时费马并没有说明原因，他只是留下这个叙述并且也说他已经发现这个定理的证明妙法，只是书页的空白处不够无法写下。始作俑者的费马也因此留下了千古的难题，三百多年来无数的数学家尝试要去解决这个难题却都徒劳无功。这个号称世纪难题的费马最後定理也就成了数学界的心头大患，极欲解之而後快。

    十九世纪时法国的法兰西斯数学院曾经在一八一五年和一八六０年两度悬赏金质奖章和三百法郎给任何解决此一难题的人，可惜都没有人能够领到奖赏。德国的数学家佛尔夫斯克尔（P_‧Wolfskehl）在1908年提供十万马克，给能够证明费马最後定理是正确的人，有效期间为100年。其间由於经济大萧条的原因，此笔奖额已贬值至七千五百马克，虽然如此仍然吸引不少的「数学痴」。

    二十世纪电脑发展以後，许多数学家用电脑计算可以证明这个定理当n为很大时是成立的，1983年电脑专家斯洛文斯基借助电脑运行5782秒证明当n为2⁸⁶²⁴³-1时费马定理是正确的（注2⁸⁶²⁴³-1为一天文数字，大约为25960位数）。

    虽然如此，数学家还没有找到一个普遍性的证明。不过这个三百多年的数学悬案终於解决了，这个数学难题是由英国的数学家威利斯（Andrew Wiles）所解决。其实威利斯是利用二十世纪过去三十年来抽象数学发展的结果加以证明。

    五０年代日本数学家谷山丰首先提出一个有关椭圆曲现的猜想，後来由另一位数学家志村五郎加以发扬光大，当时没有人认为这个猜想与费马定理有任何关联。在八０年代德国数学家佛列将谷山丰的猜想与费马定理扯在一起，而威利斯所做的正是根据这个关联论证出一种形式的谷山丰猜想是正确的，进而推出费马最後定理也是正确的。这个结论由威利斯在1993年的6月21日於美国剑桥大学牛顿数学研究所的研讨会正式发表，这个报告马上震惊整个数学界，就是数学门墙外的社会大众也寄以无限的关注。不过威利斯的证明马上被检验出有少许的瑕疵，於是威利斯与他的学生又花了十四个月的时间再加以修正。1994年9月19日他们终於交出完整无瑕的解答，数学界的梦魇终於结束。1997年6月，威利斯在德国哥庭根大学领取了佛尔夫斯克尔奖。当年的十万法克约为两百万美金，不过威利斯领到时，只值五万美金左右，但威利斯已经名列青史，永垂不朽了。

证明全过程：

http://cgd.best.vwh.net/home/flt/flt08.htm

阅读全文(2272) | 评论:0 | 复制链接