将1到N^2表示成两个长度为N的数列的正交和

最新推荐文章于 2022-10-25 22:38:26 发布

原创最新推荐文章于 2022-10-25 22:38:26 发布 · 1.3k 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#list #ui #c #output #input

数学同时被 2 个专栏收录

21 篇文章

订阅专栏

算法

16 篇文章

订阅专栏

我们可以找到两个长度为N的数列，使得从两个数列中各自挑选一个数相交得出的N^2个数正好是1到N^2之间的所有整数。

对于正样的数列，如果我们将其中一个数列中所有数都加上一个常数，另外一个数列中所有数都减去一个常数，那么必然也得到一个符合条件的数列，所以我们可以认为这样得到的数列对和原先的数列对是等价的，不在我们考虑范围之内。

问题是，相互不等价的数列对有多少个？

这个问题来源于CSDN中一道题目的讨论：

http://community.csdn.net/Expert/TopicView.asp?id=4501259

而最终的计算程序由于使用大整数，将用到GMP库。

最终结果挺有意思，同N的因子分解结果有关系。比如

N=p1^a1 *p2^a2*..*pt^at

那么结果仅仅同a1,a2,...,at有关系，而同p1,p2,...,pt没有关系。

对于整数N,对于一个序列
n0=N, n1, n2,...,nk
其中n(i+1)|n(i)而且n(i+1)>n(i), nk=1
我们称这个序列为N一个长度为k的因子序列.
比如N长度为1的因子序列只有一个为N,1.
记N长度为s的因子序列的数目为L(N,s)
那么L(N,0)=0,L(N,1)=1
上面计数问题的结果即
f(N)=Sum{L(N,k)*[L(N,k)+L(N,k-1)], for k=1,2,3,...}
比如对于N=6=2*3
长度为1的因子序列1个 {6,1}
2 2个 {6,2,1}, {6,3,1}
长度大于3的因子序列没有
所以f(6)=1*(1+0)+2*(2+1)=7
而对于N=p^k, 长度为t的因子序列相当于从p,p^2,...,p^(k-1)中选择t-1个数,共有C(k-1,t-1)个
所以L(p^k, t)=C(k-1,t-1)
f(p^k)=Sum{C(k-1,t-1)*[C(k-1,t-1)+C(k-1,t-2)], t=1,2,...,k}
=Sum{C(k-1,t-1)*C(k,t-1), t=1,2,...,k}
=Sum{C(k-1,s)*C(k,s), s=0,1,...,k-1}
=Sum{C(k-1,s)*C(k,k-s), s=0,1,...,k-1}
=C(2k-1,k) (从2k-1个样品中取k个的方案正好是从前k-1个中取s个,后k个中取k-s个,对所有可能的s进行求和)
=(2k-1)!/(k!*(k-1)!),即gxqcn的结论.

我现在写的程序就是对给定的N,先求出所有的L(N,s),然后计算出f(N).
还没有找到更加好的办法.

其实这个方法不仅仅对于N*N的方阵可用,对于将0,1,...,M*N-1填入一个M*N的矩阵的计数方案也可用.


其中L(N,s)也可以想象成下面的模型
假设
N=p1^r1*p2^r2*...*pt^rt
我们想象有一个袋子中总共装了r1+r2+...+rt个球,这些球总共有t种颜色,r1个为第一种颜色,r2个为第二种颜色,...,rt个为第t种颜色.
现在要分s次将袋子种的球取光,请问有多少种方案?
这个方案数就是L(N,s).

找到一种时间复杂度为O(L^2)次乘法的算法.
其中L=r1+r2+...+rt

定义C(m,n)=m!/(n!*(m-n)!)是组合数.

i)计算组合数C(m,n),其中m=1,2,...,L+max(r)-1. 0<=n<=m (max(r)=max(r1,r2,...,rt))
计算方法
C(1,0)=C(1,1)=1
for(m=2;m<L+max(r);m++){
C(m,0)=1;
for(n=1;n<=m;n++)C(m,n)=C(m-1,n-1)+C(m-1,n);
}
总共花费O(L^2)次加法.

然后我们定义S(N,1)=1,S(0,k)=1 S(N,k)=S(N,k-1)+S(N-1,k),
其中S(N,k)就是k个非负整数和为N的表示方法数目(k个整数有序,也就是0+1和1+0不同)
ii)计算S(n,m),其中n=0,1,2,...,L; 1<=m<=L-n+1
使用上面递推式计算,时间复杂度也是O(L^2)次加法

最后一步直接计算L(N,k)
其中L(N,1)=1
iii) L(N,k)=C(r1+k-1,k-1)*C(r2+k-1,k-1)*...*C(rt+k-1,k-1)
-Sum{ S(j+1,k-j)*L(N,k-j-1), for j=0 to k-1}
从L(N,1)开始递推计算到L(N,L),其中每步max(k,t)<=L次乘法,共L步,所以时间复杂度为O(L^2)次乘法.

有了L(N,k)后,我们就可以使用O(L)次乘法计算出f(N)了.

发现S(n,m)其实就是C(n+m-1,m-1),所以只要事先计算好组合数就够了

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <gmp.h>

int *factor_list;
int factor_count;
int L, LAR;
mpz_t *triangle;
mpz_t *f;
mpz_t N;

#define TRIANGLE(x,y) triangle[(x)*LAR+(y)]

void set_triangle()
{
int i,j;
triangle = (mpz_t *)malloc(sizeof(mpz_t)*LAR*LAR);
for(i=0;i<LAR;i++)for(j=0;j<=i;j++)mpz_init(TRIANGLE(i,j));
mpz_set_ui(TRIANGLE(0,0),1);
mpz_set_ui(TRIANGLE(1,0),1);
mpz_set_ui(TRIANGLE(1,1),1);
for(i=2;i<LAR;i++){
mpz_set_ui(TRIANGLE(i,0),1);
mpz_set_ui(TRIANGLE(i,i),1);
for(j=1;j<i;j++){
mpz_add(TRIANGLE(i,j),TRIANGLE(i-1,j-1),TRIANGLE(i-1,j));
}
}
}

void calc_f(){
int i,k;
mpz_t M,R;
mpz_init(M);
mpz_init(R);
f= (mpz_t *)malloc(sizeof(mpz_t)*(L+1));
for(i=0;i<=L;i++)mpz_init(f[i]);
mpz_set_ui(f[1],1);
printf("L[1]=1/n");
for(k=2;k<=L;k++){
mpz_set_ui(M,1);
for(i=0;i<factor_count;i++){
mpz_mul(M,M,TRIANGLE(factor_list[i]+k-1,k-1));
}
for(i=0;i<k-1;i++){
mpz_set(R,TRIANGLE(k,i+1));
mpz_mul(R,R,f[k-i-1]);
mpz_sub(M,M,R);
}
mpz_set(f[k],M);
printf("L[%d]=",k);
mpz_out_str(stdout,10,f[k]);
printf("/n");
}
mpz_clear(R);
mpz_clear(M);
}

void output_result()
{
mpz_t M,R;
int i;
mpz_init(M);
mpz_init(R);
mpz_set_ui(M,1);
for(i=2;i<=L;i++){
mpz_set(R,f[i]);
mpz_add(R,R,f[i-1]);
mpz_mul(R,R,f[i]);
mpz_add(M,M,R);
}
printf("The result is:/n");
mpz_out_str(stdout,10,M);
printf("/n");
mpz_clear(R);
mpz_clear(M);
}

#define PRIME_LIMIT 1048576
#define SQR_LIMIT 1024
int prime_flag[PRIME_LIMIT];
int *prime_list;
int prime_count;
void init_prime(){
int i,count;
memset(prime_flag,-1,sizeof(prime_flag));
prime_flag[0]=prime_flag[1]=0;
for(i=2;i<=SQR_LIMIT;i++){
if(prime_flag[i]){
int j;
for(j=i*i;j<PRIME_LIMIT;j+=i)
prime_flag[j]=0;
}
}
for(i=2,count=0;i<PRIME_LIMIT;i++){
if(prime_flag[i])count++;
}
prime_count=count;
prime_list=(int *)malloc(sizeof(int)*count);
for(i=2,count=0;i<PRIME_LIMIT;i++){
if(prime_flag[i]){
prime_list[count++]=i;
}
}
}

void defactor()
{
int i,pc,maxr;
mpz_t R,M;
mpz_init(R);
mpz_init(M);
mpz_set(M,N);
for(i=0,pc=0;i<prime_count;i++){
if(!mpz_mod_ui(R,M,prime_list[i])){
pc++;
do{
mpz_fdiv_q_ui(M,M,prime_list[i]);
}while(!mpz_mod_ui(R,M,prime_list[i]));
if(!mpz_cmp_si(M,1))break;
}
}
if(mpz_cmp_si(M,1)){
int r;
if(r=mpz_probab_prime_p(M,5)){//Using prime test to verify the left is a prime
pc++;
if(r==1){
fprintf(stderr,"WARNING:The left factor could not be verified to be prime but looks like to a prime. It is treated as a prime/n");
}
}else{
fprintf(stderr,"Cannot find all prime factors of input/n");
exit(-1);
}
}
factor_list=(int *)malloc(sizeof(int)*pc);
factor_count=pc;L=0,maxr=0;
mpz_set(M,N);
for(i=0,pc=0;i<prime_count;i++){
if(!mpz_mod_ui(R,M,prime_list[i])){
factor_list[pc]=0;
do{
mpz_fdiv_q_ui(M,M,prime_list[i]);
factor_list[pc]++;
}while(!mpz_mod_ui(R,M,prime_list[i]));
L+=factor_list[pc];
if(maxr<factor_list[pc])maxr=factor_list[pc];
pc++;
if(!mpz_cmp_si(M,1))break;
}
}
if(mpz_cmp_si(M,1)){
factor_list[pc++]=1;
L++;
if(maxr<1)maxr=1;
}
LAR=L+maxr+1;
mpz_clear(R);
mpz_clear(M);
}

int
main(void)
{
mpz_init(N);
printf("Please input the large number:/n");
mpz_inp_str(N,stdin,10);
init_prime();
defactor();
set_triangle();
calc_f();
output_result();
mpz_clear(N);
return 0;
}