博文
中国剩余定理与大衍求一术(2006-11-05 15:27:00)
摘要：前面提到的‘中国剩余定理’是老外对我们中国古代数学家对同余理论作的贡献的肯定而命名的，其实它包括两部分，一是《孙子算经》中的定理，二是这里所谫的‘大衍求一术’，但是总的来说，它是对一次同余式组的一类通用解法。
《孙子算经》中有物不知其数问题，它表示成现代数学形式是：
N=R1(Mod 3)
N=R2(Mod 5)
N=R3(Mod 7)
书中给出了这组方程的通解：
N=70R1+21R2+15R3-105P
并副上一首诗歌：
“
三人同行七十（７０）稀，
五树梅花二一（２１）枝。
七子团圆正半月（１５），
除百零五（１０５）便得知。
”
至于为什么要这么算，书中未作任何说明，实际上这一问题还是未解决。真正从完整的计算理论上解决这个问题的，是南宋时期的数学家秦九韶，在他的《数书九章》中提出的所谓‘大衍求一术’就是对这一问题的一般解法。
其实原理非常简单，拿前面的通式说明，为什么R1前的系数是70，用手算一下，它被5和7整除，但是被3除余1，也就是70=1 (mod 3)，于是70R1=R1(mod 3)，而对其它模为0，作成的和不影响其它的数，所以N是问题的解，可以从空间的角度理解，70,21,15就是那个通式方程组的基，随着R1，R2，R3取得不同，解都可以由他们的线性组合构成，而105是3，5，7的最小公倍数，即模他们都为0，是解的循环元。
而所谓的‘大衍求一术’就是求满足x=1(mod M1)且x|M2,x|M3...的x，方法是：先取(M2,M3,...Mk)的最小公倍数a，然后解方程ax=1(mod M1)，那么ax就是我们要求的R1前的系数，整体思想就是这样的，而在秦九韶的方法中非常复杂，但是拿现代的眼光看，这个问题已经解决，即我们熟悉的欧几里德算法（->EUCLID算法与RSA）。不再赘述。......
Romberg算法求数值积分的原理与实现(2006-10-18 22:52:00)
摘要：数值积分实际上都是基于插值的，我们算的都是离散的一些个点与点对应的函数值，于是要求连续的积分是不可能的。
记等距分点a<x0<x1<...<b，距离h=(b-a)/n，n为分段数，设它的插值函数存在且为L(x)，我们用SL(x)dx的积分代替Sf(x)dx的积分值，这就是数值积分的最基本原理。
Newton-Cores公式，给出上面的积分为：
In=(b-a)∑[ C(n,k)f(xk) ]
其中C(n,k)是Cores系数，设x=a+th，则有
C(n,k)=(-1)^(n-k)/(nk!(n-k)! S∏(t-j)dt j=0 to n,j!=k
当n=1时，算出C(1,0)=C(1,1)=1/2，即熟悉的梯形公式
I1=(b-a)(f(a)+f(b))/2
当n=2时，算出C(2,0)=1/6,C(2,1)=4/6,C(2,2)=1/6，这时就是Simpson公式，记得原来数学分析书上有提到过，它是按二次曲线插值的，有比较好的效率
I2=(b-a)(f(a)+4f((a+b)/2)+f(b))/6
这样一直算下去，会有更好的公式，且当n为偶数时，具有n+1次代数精度，但是当n>=8时，系数出现负项，将会不稳定，因此实际中只使用低次的Newton_Cores公式。
解决的办法有几种，一种是进行复化的求积公式，跟低次分段插值一样，将积分区间分成若干小段，每小段做低次积分，再求和；另一种办法就是用所谓的变步长积分。
其实所谓变步长法，就是我们熟知的递推法，在这里将看到递推法的妙处。我们拿到一个积分，和被积区间(a,b)，如何确定h呢？递推，先用较大的h试，如果不满意就将它细分，直到满意为止。
如果一开始只有一个分点，再将它们二分得四个，一直二分下去，那么原梯形公式将得到递推梯形公式，它是从n个分段二分到2n个的递推关系：
T2n=Tn/2+h/2∑f(xk+1/2) [1]
右边求和部分是新插入的结点，于是就很方便程序控制，误差控制就考虑|T2n-Tn|的值，事后估计误差。但是直接用梯形法，收敛速度很慢，进一步优化发现，梯形法的余项展开为：
T(h)=I+a1h^2+a2h^4+a3h^6+...
如果h比较小，那T(h)=I+O(h^2)，为了提......
初等数论基础知识(2006-10-15 21:55:00)
摘要：以前听老师说，初等数论要到大四学，今天翻到一本书，不知道讲得全不全，原来就是余数定理那一块知识，回想起来初中的时候老爹就教过我，那个时候他说的是‘中国余数定理’。
这个定理呢是老外的叫法，叫“Chinese remainder theorem”，在中国叫“孙子定理”，源于《孙子算经》（这本书的作者是不是叫孙子，不清楚，应该不会，谁愿意叫这么个名儿），书中有这样一个题：“今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？”，小时候老爹跟我说的时候，我觉得这还不简单，不就列几个方程呗，然后就不知道怎么列，还要设商，三个商不知道，加一个未知数，四个未知数三个方程，怎么解？
我们记这样一个运算：Mod，a Mod b的结果表示a除以b的余数，且是一个非负的小于b的数（a,b皆非负）。简单的理解就是，模运算是将整数打个折都射到[0,1,...b-1]这样一个区间上。
另外还有这样一个关系：≡，表示同余关系，a≡b(Mod m)就是说a Mod m=b Mod m，比如7≡2(Mod 5)。我们在自然数集上考虑这个关系，它满足：
1）自反性
a≡a(Mod m)
2）对称性
if a≡b(Mod m)，then b≡a(Mod m)
3）传递性
if a≡b(Mod m)，b≡c(Mod m)，then a≡c(Mod m)
也就是说它是一种‘等价’关系，它划分了[0],[1],...[m-1]这样m个等价类，也叫同余类。
同余式有一些非常有用的运算定律：
1）if a≡b(Mod m) and c≡d(Mod m) , then a+c≡b+d(Mod m)
2）if a≡b(Mod m) and c≡d(Mod m) , then ac≡bd(Mod m)
3）if a≡b(Mod m) k>0 , then ak≡bk(Mod mk)
4）if a≡b(Mod m) and d|a,d|b,d|m , then a/d≡b/d(Mod m/d)
5）if a≡b(Mod m) and d|m , then a≡b(Mod d)
6）if ca≡cb(Mod m) and d=gcd(c,m) , then a≡b(Mod m/d)......
Euclid算法与RSA(2006-10-11 09:48:00)
摘要：历史上第一个称得上算法的好像就是这个欧几里得算法，其实就是地球人都知道的辗转相除，不要小看她，她是很美的。
简单的描述就是，记gcd(a,b)表示非负整数a,b的最大公因数，那么：gcd(a,b)=gcd(b,a%b)或者gcd(a,0)=gcd(0,a)=a。
写成程序很简单，不管是用递归还是循环：
int gcd(int a,int b)
{
if(a==0)
return b;
if(b==0)
return a;
return gcd(b,a%b);
}
不仅算法形式简单，而且效率很高，我不知道具体是多少复杂度的，我只知道效率很高；）
前天看ＲＳＡ算法，是非对称加密的标准算法，其实算法很简单：
找到两个素数p,q，再找一个数r，使gcd(r,(p-1)(q-1))=1，也就是说互素，然后再找一个数m，使rm=1(mod (p-1)(q-1))，然后再作乘法n=pq，然后把pq丢掉，最好是让任何人都不知道，包括自己（免得说梦话的时候被人听到），然后手里拿到r,m,n，r就是Private Key，只有你知道，而m,n就是Public Key。设信息为a,加密过程：a^r=b (mod n)，b就是密文，解密过程：b^m=a(mod n)，反过来用m加密，用r解密是一样的。
书上说由gcd(r,(p-1)(q-1))=1到求m，使rm=1(mod (p-1)(q-1))是很容易的，就用辗转相除，我想了好久才想到一个方法。
问题：如果gcd(a,b)=1，求x，使ax=1(mod b)
......
STL与泛型(2006-09-28 20:06:00)
摘要：我真的不是一个好的C++爱好者，可能我还不算真正懂C++，我常想什么是C++，以前看到有人这样写程序：
#include <iostream>
using namespace std;
int main(int argc, char* argv[])
{
cout<<"Hello!World!"<<endl;
return 0;
}
然后就觉得好酷，我不用像C那样printf()，于是有段时间就以为C++就是cout,cin，后来又学了写class，于是就认为C++就是在一个程序里写几个自己的class，输入输出用iostream，呵，C++！
再后来学了一点MFC，了解到Microsoft的人是怎么写基础类库的，被那种完整的类体系迷倒，老实说，MFC对SDK的封装简直太彻底了！更重要的是，从MFC的设计者角度，终于明白了一点点，什么是OOD，什么是OOA，如何设计类，如何有效的封装。于是站在MFC肩上，做Windows程序开发，就是从MFC中选择适当的类，继承并设计自己的子类，再外加上VC++提供了Class Wizard，太方便了。于是就以为那就是C++。
感觉从C，C++，每一次大的变革，都是从思想上的变革，软件史上的革命要么是数学上的革命，比如人工智能的突破性进展等，要么就是思想上的革命，思想的飞跃造就了一个IT时代。同样的问题，思考的角度不同，那它本身就不同了。
泛型程序设计，这个词对我还很新，如果说过程化程序设计面向的是过程是函数，如果说C++面向的是对象，那泛型面向的会是什么呢？泛型，我感觉就是更广义的抽象，泛化，数组，链表等泛化成容器，指针泛化成迭代器，函数泛化成算法，突然一下整个C++完全变了样，原来那些东西都不见，感觉像一门新的语言！像这样：
// stl.cpp : Defines the entry point for the console application.
//
#include "StdAfx.h"
#include <iostream>
#include <algorithm>
#include <functional>
#include <str......
香农编码(2006-09-20 20:52:00)
摘要：今天的作业，想到可以直接由浮点数格式中取出尾数得到ShannanCode，麻烦点的地方是需要重新对阶，因为在尾数规格化省掉了一位，最终效率比较高。
#ifndef SHANNAN_CODE_DEF_200609
#define SHANNAN_CODE_DEF_200609
typedef unsigned long U32;
typedef signed char S8;
class CShannan
{
float *m_p;//消息概率
float *m_sp;//累计概率P
float *m_i;//信息量
int *m_length;//码长
U32 *m_code;
int m_n;
public:
static int cpux86();
public:
CShannan(int n);
~CShannan();
void InputPro();
void OutputCode();
};
#endif
#include "stdio.h"
#include "assert.h"
#include "math.h"
#include "shannan.h"
#ifndef NULL
#define NULL 0
#endif
CShannan::CShannan(int n)
{
assert(n>0);
m_n=n;
m_p=new float[n];
m_sp=new float[n];
m_i=new float[n];
m_length=new int[n];
m_code =new U32[n];
assert(m_p!=NULL);
assert(m_sp!=NULL);
assert(m_i!=NULL);
assert(m......
火焰效果（1）(2006-09-19 20:49:00)
摘要：昨天晚上，三本书，翻过来翻过去，写下了下面的程序（部分）：
int CMyFlashDlg::OnCreate(LPCREATESTRUCT lpCreateStruct)
{
if (CDialog::OnCreate(lpCreateStruct) == -1)
return -1;
// TODO: Add your specialized creation code here
CClientDC dc(this);
CRect rect;
GetClientRect(&rect);
m_nWidth=rect.Width();
m_nHeight=rect.Height();
m_mbm1.CreateCompatibleBitmap(&dc,m_nWidth/3,m_nHeight/3);
m_mbm2.CreateCompatibleBitmap(&dc,m_nWidth/3,m_nHeight/3);
m_mdc1.CreateCompatibleDC(&dc);
m_mdc2.CreateCompatibleDC(&dc);
m_mdc1.SelectObject(m_mbm1);
m_mdc2.SelectObject(m_mbm2);
SetTimer(IDT_MAIN_TIMER,30,NULL);
return 0;
}
void CMyFlashDlg::myfire()
{
CClientDC dc(this);
int sw=m_nWidth/3,sh=m_nHeight/3;
for(int y=1;y<sh;++y)
for(int x=1;x<sw;++x)
{
COLOR......
牛顿插值算法与实现(2006-09-15 00:09:00)
摘要：牛顿真是牛，拉格朗日插值法只能算是数学意义上的插值，从插值基函数的巧妙选取，已经构造性的证明了插值法的存在性和惟一性，但是从实现的角度看并不很好，而牛顿很好的解决了这个问题。
牛顿插值是基于下面这些的公式：
f[x0,x1,...xk]=(f[x1,...xk]-f[x0,...xk-1])/(xk-x0)
f[x]=f(x)
f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0)...(x-xn-1)+Rn(x)
前两个是均差的递推关系式，而后一个就是牛顿插值公式，其中N(x)=f(x)-Rn(x)，即目标多项式，Rn(x)是n阶插值余项，我们就是用N(x)去近似f(x)。
可以构造这样一个均方差表：
xk f(xk) 一阶均差 二阶均差 ...
x0 f(x0)
x1 f(x1) f[x0,x1]
x2 f(x2) f[x1,x2] f[x0,x1,x2]
...
如果有n个点插值，表会有(n*n)/2+n个表项，如果直接编程会有O(n*n)的空间复杂度，编程时做个简单的改进，不难发现在这个表中只有部分数据有用，对角线（斜行）它们是目标值，用来表示多项式的，左边的两纵行（实际上只需要x一行）以及最底下的一行，表示当前插值的状态。经过改进后只需要O(n)的空间复杂度。
两个过程：
1，新增加一个点时的更新。只须更新最底下一行数据，其递推关系由均差公式给出，最后算出高一队的均差值，需时O(n)
2，插入点完成后如何计算多项式在另外给定点的值N(x)。
由牛顿插值公式，最终的表达式为：
N(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0)...(x-xn-1)
如果直接将它展开，再算实在麻烦，实际上大可不必这样做，还记得多项式求值的秦九......
八数码问题（一）(2006-09-12 22:55:00)
摘要：原来粗略的考虑，h函数选各个数码离目标位置的最短距离的和，原则上是可行的，但启发性不大，举个简单例子：
{2,1,3,
8,0,4,
7,6,5}
和目标状态只差1和2的位置交换一下就对了，h函数给出的估计值只有2，实际上却差得太远，如果你小时候有玩过这样的‘数字拼盘’游戏，就会知道，像这样的是没有解的，至少那时候的感觉是这样，就算有解，也需要移很多次。
基于这样的h函数，感觉启发性不够，在海量搜索的时候效果不是很好，像有的CASE，搜到OPEN表+CLOSE表共好几万的结点了还在不停的搜。可见八数码问题有难有易，虽然9!只有三十多万，也比较难搜。
今天写的程序：
#include <iostream.h>
struct Snode
{
int key;
int map[9];
int g,f;
int last;
public:
Snode(){key=0;}
void TransformIn(const int *d);
void TransformOut(int *d);
};
void Snode::TransformIn(const int *d)
{
key=0;
for(int i=0;i<9;++i)
key=key*9+(map[i]=d[i]);
}
void Snode::TransformOut(int *d)
{
for(int i=0;i<9;++i)
{
d[i]=map[i];
}
}
int spath[9][9]=
{{0},
{0,1,2,1,2,3,2,3,4},
{1,0,1,2,1,2,3,2,3},
{2,1,0,3,2,1,4,3,2},
{3,2,1,2,1,0,3,2,1},
{4,3,2,3,2,1,2,1,0},
{3,2,3,2,1,2,1,0,1},
{2,3,4,1,2,3......
A*高效搜索算法(2006-09-11 17:31:00)
摘要：A*高效搜索算法 2006/09/11 rickone
了解了基本搜索算法，下面就来看A*，神奇的A*。（--->搜索方法小结）
A*是一种启发式搜索，一种有序搜索，它之所以特殊完全是在它的估价函数上，如果我要求的是从初始结点到目的结点的一个最短路径（或加权代价）的可行解，那对于一个还不是目标结点的结点，我对它的评价就要从两个方面评价：第一，离目标结点有多近，越近越好；第二，离起始结点有多远，越近越好。记号[a,b]是表示结点a到结点b的实际最短路径代价。设起始结点为S，当前结点为n，目标结点为G，于是n的实际代价应该是f*(n)=g*(n)+h*(n)，其中g*(n)=[S,n],h*(n)=[n,G]，对于是g*(n)是比较容易得到的，在搜索的过程中我们可以按搜索的顺序对它进行累积计算，当然按BFS和DFS的不同，我们对它的估价g(n)可以满足g(n)>=g*(n)，大多可以是相等的。但是对于h*(n)我们却了解得非常少，目标结点正是要搜索的目的，我们是不知道在哪，就更不知道从n到目标结点的路径代价，但是或多或少我们还是可以估计的，记估价函数f(n)=g(n)+h(n)。
我们说如果在一般的图搜索算法中应用了上面的估价函数对OPEN表进行排序的，就称A算法。在A算法之上，如果加上一个条件，对于所有的结点x，都有h(x)<=h*(x)，那就称为A*算法。如果取h(n)=0同样是A*算法，这样它就退化成了有序算法。
A*算法是否成功，也就是说是否在效率上胜过蛮力搜索算法，就在于h(n)的选取，它不能大于实际的h*(n)，要保守一点，但越接近h*(n)给我们的启发性就越大，是一个难把握的东西。
A*算法流程：
首先将起始结点S放入OPEN表，CLOSE表置空，算法开始时：
1、如果OPEN表不为空，从表头取一个结点n，如果为空算法失败
2、n是目标解吗？是，找到一个解（继续寻找，或终止算法）；不是到3
3、将n的所有后继结点展开，就是从n可以直接关联的结点（子结点），如果不在CLOSE表中，就将它们放入OPEN表，并把S放入CLOSE表，同时计算每一个后继结点的估价值f(n)，将OPEN表按f(x)排序，最小的放在表头，重复算法到1
最短路径问题，Dijkstra算法与A*
A*是求这样一个和最短路径有关的问......