[难题征解] 反求线性相关数组(续)

gxqcn 2006-01-06 10:34:52
加精
>> I 前言
================
为了活跃论坛气氛,抛出一个难题。该话题曾在 2004 年被讨论过,但并没有得到彻底解决(为防论坛强制结帖而提前结帖),有兴趣的朋友请先浏览:http://topic.csdn.net/TopicFiles/2004/06/29/13/3131084.xml

>> II 引子
================
为了引出该问题,特自定义一个概念以方便阐述(可能叙述得不够严谨,若与现有的术语相冲突,欢迎指正):

【线性相关数组】已知一个二维数组 C={c[i][j]}mxn,如果存在两个一维数组 A={a[i]}m、B={b[j]}n,使得 c[i][j]=a[i]+b[j] (0≤i≤m-1,0≤j≤n-1)均成立,则称 C 与 A、B 线性相关,记作“C = A⊙B”。

>> III 问题
================
已知二维数组 C 的某个排列 C'(即其 mxn 个元素被随机打乱了),且其所有元素为非负整数,且两两不等,最小元素为 0。请反求出对应的 A、B,使得 C = A⊙B(其中 A、B 均为非递减的非负整数数组;如果 m==n,则附加要求 A≤B(按字典排列))。

可简化条件:已确知一定存在 A、B,使 C = A⊙B 成立;

备注:不要求给出代码,但须将算法思路阐明。

>> IV 讨论
================
1、如果有人可迅速判定“是否存在 A、B,使得 C = A⊙B 成立”的快速算法,欢迎提供;
2、如果满足上述“简化条件”,(A, B) 是否具有唯一性,也欢迎讨论。

>> V 例子
================
例如,已知 3x3 的矩阵 { { 3, 8, 1 }, { 2, 4, 6 }, { 7, 0, 5 }},
应求得 A={ 0, 1, 2 },B={ 0, 3, 6 }.
...全文
907 62 打赏 收藏 转发到动态 举报
写回复
用AI写文章
62 条回复
切换为时间正序
请发表友善的回复…
发表回复
yuxh 2006-05-18
  • 打赏
  • 举报
回复
ohh!我的想法与zzwu(未名)的是一样的
yuxh 2006-05-18
  • 打赏
  • 举报
回复
int get_next(int *idx, int len, int cur, int n)
{
int i;
for(i=cur-1; i < len && n > 0; n--)
for(i++; i<len && idx[i] == 1; i++);
if(i >= len) return -1;
else return i;
}
int get_val_idx(int *ary, int begin, int len, int *idx, int val)
{
int i;
for(i=begin; i<len && (idx[i] == 1 || ary[i] != val); i++);
if(i >= len) return -1;
else return i;
}
int split(int *ary, int *idx, int *dst, int m, int n, int i, int j)
{
int len = m*n, *idx_bak, next, k, t, val, val_idx, min_depth, n_idx, min;
if((next = get_next(idx, len, 1, 1)) == -1) return 0;

idx_bak = (int *)malloc(sizeof(int) * len);
memcpy(idx_bak, idx, sizeof(int) * len);

if(n==i)
min = dst[n-1];
else if(m==j)
min=dst[(m-1)*n];
else {
min_depth = n-i < m-j ? n-i : m-j;
if((n_idx = get_next(idx, len, next, min_depth)) == -1) goto FAIL;
min = ary[n_idx];
}
if((n_idx = get_next(idx, len, next, m+n-i-j)) == -1) goto FAIL;
min += ary[n_idx] - ary[0];
if(min > ary[m*n-1]) goto FAIL;

if(j < m) {
dst[j*n] = ary[next];
idx[next] = 1;
for(k=1; k<i; k++) {
val = dst[(j-1)*n+k] + dst[j*n] - dst[(j-1)*n];
if((val_idx = get_val_idx(ary, next, len, idx, val)) == -1) goto NEXT;
dst[j*n+k] = ary[val_idx];
idx[val_idx] = 1;
}
if(split(ary, idx, dst, m, n, i, j+1) == 0) goto SUCCESS;
}

NEXT:
if(i < n) {
dst[i] = ary[next];
idx[next] = 1;
for(k=1; k<j; k++) {
val = dst[k*n+i-1] + dst[i] - dst[i-1];
if((val_idx = get_val_idx(ary, next, len, idx, val)) == -1) goto FAIL;
dst[k*n+i] = ary[val_idx];
idx[val_idx] = 1;
}
if(split(ary, idx, dst, m, n, i+1, j) == 0) goto SUCCESS;
}
FAIL:
memcpy(idx, idx_bak, sizeof(int) * len);
free(idx_bak);
return -1;
SUCCESS:
free(idx_bak);
return 0;
}
对于36*36的矩阵,可以很快地得到:
Array A:
0 81 324 486 648 1296 1458 2349 2673 3078 4050 4941 5346 6318 7533 8100 8586 10125 10449 11988 12474 13041 14256 15228 15633 16524 17496 17901 18225 19116 19278 19926 20088 20250 20493 20574
Array B:
0 1 4 6 8 16 18 29 33 38 50 61 66 78 93 100 106 125 129 148 154 161 176 188 193 204 216 221 225 236 238 246 248 250 253 254
zzwu 2006-01-18
  • 打赏
  • 举报
回复
gxqcn:

不知道本题有什么实用价值? 已发现的.
zzwu 2006-01-17
  • 打赏
  • 举报
回复
to mathe():
你说的有理,行列均匀不均匀并不影响效率,搜索一样快.
我没有看你们前面的讨论,是因为内容太多,看不胜看,
要看懂全部细节也决不容易,你们已经做了极为全面和
深入的考虑,要判别我的方案有否新东西一定很容易,
所以就把它端了出来,没有别的意思.
gxqcn 2006-01-16
  • 打赏
  • 举报
回复
该问题至此已圆满解决,非常感谢 mathe!
本贴将于春节前后结帖,大家若有好的想法可以继续。
mathe 2006-01-16
  • 打赏
  • 举报
回复
to zzwu

2006-1-11 16:27:17
贴出的代码就是你上面的分析方法.不过还是你描述的比较详细.
至于优先搜索"均匀的"方案都没有必要,搜索的算法很快的
zzwu 2006-01-16
  • 打赏
  • 举报
回复
以上是一种想法, 它企图一下子确定整个B数组,为此,
要从C'中,寻找差为A[1]-A[0]的所有数对,
如果C'的元素N*M太大,这种数对的数目可能很多,如何组成循环,要考虑。

下面还有一种想法是: 
将C'中元素,由小到大,一个一个地试放到矩阵中,在此过程中,
为了保证所需性质,会自动用去一些元素。仍举上例说明之.

原C'={0,1,2,3,4,5,8,9,10,12,15,17,18,19,20,29},
过程是:

1: 把C'中最小元0到第一行第一列:
A\B┃ B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0
A[1]┃
A[2]┃
A[3]┃

新C'={1,2,3,4,5,8,9,10,12,15,17,18,19,20,29},


接下来,考察新C'中的第一个元素1,它可以放第一列或第一行,
考虑到对称性,不妨把1放在第一列:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0
A[1]┃ 1
A[2]┃
A[3]┃
新C'={2,3,4,5,8,9,10,12,15,17,18,19,20,29},


3: 接下来是放新C'中的第一个元素2,它可以放第一列或第一行,即放成
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0
A[1]┃ 1
A[2]┃ 2
A[3]┃

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2
A[1]┃ 1
A[2]┃
A[3]┃


我们优先试后一情况,也就是使行与列长度均匀,结果应为:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2
A[1]┃ 1 3
A[2]┃
A[3]┃
新C'={4,5,8,9,10,12,15,17,18,19,20,29},

4.接下来是放其中的第一个元素4,同样, 4可以第一列如果这样,则结果应是:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2
A[1]┃ 1 3
A[2]┃ 4 6
A[3]┃
但这不行,因为没有6;


所以,4只可以放在第一行,这样就得:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4
A[1]┃ 1 3 5
A[2]┃
A[3]┃
这样确实行,因有确实有5。
新C'={8,9,10,12,15,17,18,19,20,29},


5. 接着再确定新C'中的第一个元素8如何放,它可以是第一列或第一行,
  如为第一列,则结果应是:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4
A[1]┃ 1 3 5
A[2]┃ 8 10 12
A[3]┃


如放第一行,则结果应是:

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4 8
A[1]┃ 1 3 5 9
A[2]┃
A[3]┃

两种方案都行,先试前一种,这时
新C'={9,15,17,18,19,20,29},

6. 接着确定其中的第一个元素9, 可以放成
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4
A[1]┃ 1 3 5
A[2]┃ 8 10 12
A[3]┃ 9 11 13
但此不行,因没有11,也没13

可以改放
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4 9
A[1]┃ 1 3 5 10
A[2]┃ 8 10 12 17
A[3]┃
显然也不行,因为没有10,17


这样,也就是否定了第5步的两种方案。
但第5步是第4步的一种方案继承下来的,
所以要再退到第3步,第3步还有一种方案是:

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0
A[1]┃ 1
A[2]┃ 2
A[3]┃

这时新C'={3,4,5,8,9,10,12,15,17,18,19,20,29},

确定新C'中的第一个元素3如何放,它可以是第一列或第一行,
优先放在行上,使行列均匀,则结果应为:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 3
A[1]┃ 1 4
A[2]┃ 2 5
A[3]┃

这时新C'={8,9,10,12,15,17,18,19,20,29},

按类似方法,可以把整个A[],B[]都找出,如果没有解,同样可以知道。









mathe 2006-01-16
  • 打赏
  • 举报
回复

#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;
}

mathe 2006-01-16
  • 打赏
  • 举报
回复
试着使用一下gmp,使用我上面提到的O(L^2)的算法,可以解决在N非常非常巨大时,使用
0..N*N-1填充符合条件的正方形的数目:
例子:
Please input the large number:
54432000000000
L[1]=1
L[2]=1798
L[3]=410403
L[4]=31857996
L[5]=1214965005
L[6]=27422340066
L[7]=408775871903
L[8]=4320097285208
L[9]=33984291814461
L[10]=206058188845110
L[11]=988087528491591
L[12]=3820069601860788
L[13]=12082187648562993
L[14]=31606430046799986
L[15]=68938028561606715
L[16]=126077161204620528
L[17]=194005583715214194
L[18]=251550075342977964
L[19]=274681887121543958
L[20]=251924874608623160
L[21]=193081337853227502
L[22]=122676066960417036
L[23]=63856316507059626
L[24]=26770562107016112
L[25]=8816110383325300
L[26]=2195812819915000
L[27]=388805935131900
L[28]=43610194227600
L[29]=2329089562800
The result is:
623834616899829926853986500935358200
zzwu 2006-01-16
  • 打赏
  • 举报
回复
我有一个想法.
但前面内容太多,我只看了题目,其余完全未看,不知道我的想法是否新的.
通过前面的例子来说明.

已知:
C'={0,1,2,3,4,5,8,9,10,12,15,17,18,19,20,29},
求A(,,),B(,,),使C=A⊙B

1. 从C'中找出第1,2元素,作为A(0),A(1),即:
A\B┃ B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0
A[1]┃ 1
A[2]┃
A[3]┃

2.下面来确定第1,2行的其余元素,为了简单,把这些元素记为
C01 C02 C03和C11 C12 C13,即:
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 C01 C02 C03
A[1]┃ 1 C11 C12 C13
A[2]┃
A[3]┃

显然,任意一列,第2行与第1行对应列元素之差必须=A[1]-A[0]=1-0=1,
所以,我们从C'除0,1之外的其余元素中寻找差=1的所有数对,它们是
(2,3)(3,4)(4,5)(8,9)(9,10)(17,18)(18,19)(19,20) ---(A)
第1,2行的元素只能从以上数对中取.
为此,我们先取最前面3对:
(2,3)(4,5)(8,9) ---(B)[注]
将它们试放到矩阵中,得

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 2 4 8
A[1]┃ 1 3 5 9
A[2]┃
A[3]┃

[注]取数对时,不允许共同元素,所以,取(2,3),就不能取(3,4)

3. 下面要在C'中寻找第2,3行的元素,它们相邻列元素之差应和第0行(或
第1行)相邻列元素之差相同,即2(=2-0),2(=4-2),4(8-4).

因C'其余元素
{10,12,15,17,18,19,20,29},

A[2]必须其中最小元素10, 故第2行其余的元素必须是10+2=12,12+2=14,14+4=18,
但其中元素14在上面C'的剩余元素集中并不存在,所以,我们知道,从(A)中
取(B)行不通,故要改取其他的数对来试,

我们将(B)中第1对(2,3)去掉,改为(3,4),
第2对必须是(8,9),(注:因4已用去,故不能取(4,5)),
同理,第3对必须是(17,18),填入矩阵:

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 3 8 17
A[1]┃ 1 4 9 18
A[2]┃
A[3]┃

下面要在C'中寻找第2,3行的元素,它们相邻列元素之差应和第0行(或
第1行)相邻列元素之差相同,即3(=3-0),5(=8-3),9(17-8).

因C'其余元素{2,5,10,12,15,19,20,29},

A[2]必须=2, 故第2行的其余元素必须是2+3=5,5+5=10,10+9=19,
但这些元素2,5,10,19在上面都存在,所以,我们可以将它们填入,得
A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 3 8 17
A[1]┃ 1 4 9 18
A[2]┃ 2 5 10 19
A[3]┃

C'剩下元素{12,15,20,29},

显然,它们的差依次也是3,5,9,所以正好可以放入:

A\B┃B[0] B[1] B[2] B[3]
━━╋━━━━━━━━━━
A[0]┃ 0 3 8 17
A[1]┃ 1 4 9 18
A[2]┃ 2 5 10 19
A[3]┃ 12 15 20 29

程序结束.


mathe 2006-01-13
  • 打赏
  • 举报
回复
改成用long long,可以计算更大的范围了:
Please input N (N>1):
11520
11520 has total 54 factors,maximal factor chain length 11
Result:4357807740
Please input N (N>1):
34560
34560 has total 72 factors,maximal factor chain length 12
Result:135778680180
Please input N (N>1):
172800
172800 has total 108 factors,maximal factor chain length 13
Result:10084988190420
mathe 2006-01-13
  • 打赏
  • 举报
回复
现在我找到一种算法,速度还是挺快的.
不过对于3840已经有85649850种结果了,对于更大的N就不好办了,要溢出了
Please input N (N>1):
128
128 has total 8 factors,maximal factor chain length 7
Result:1716
Please input N (N>1):
3840
3840 has total 36 factors,maximal factor chain length 10
Result:85649850
Please input N (N>1):
1920
1920 has total 32 factors,maximal factor chain length 9
Result:15291276
Please input N (N>1):
960
960 has total 28 factors,maximal factor chain length 8
Result:2628780
mathe 2006-01-13
  • 打赏
  • 举报
回复
发现S(n,m)其实就是C(n+m-1,m-1),所以只要事先计算好组合数就够了
mathe 2006-01-13
  • 打赏
  • 举报
回复
找到一种时间复杂度为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)了.
mathe 2006-01-13
  • 打赏
  • 举报
回复
其中L(N,s)也可以想象成下面的模型
假设
N=p1^r1*p2^r2*...*pt^rt
我们想象有一个袋子中总共装了r1+r2+...+rt个球,这些球总共有t种颜色,r1个为第一种颜色,r2个为第二种颜色,...,rt个为第t种颜色.
现在要分s次将袋子种的球取光,请问有多少种方案?
这个方案数就是L(N,s).
mathe 2006-01-13
  • 打赏
  • 举报
回复
一点也不错.
我找到的结果是这样的:

对于整数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的矩阵的计数方案也可用.


gxqcn 2006-01-13
  • 打赏
  • 举报
回复
两个方阵中如果一个方阵可通过有限次的行对换、列对换、转置可转化为另一个方阵,则称这两个方阵【同构】;否则称【不同构】。

为严谨起见,上贴相应处应修订为:
“记 N 阶连续整数方阵可重排成互不同构“加法表”方阵的数目为 f(N)”
gxqcn 2006-01-13
  • 打赏
  • 举报
回复
提出这个计数问题很有意义,但通项公式确很难给出,它似乎涉及到数论(整除性)和组合数学两方面。

如果将数字
0,1,2,...,N*N-1
填充如一个N*N的方阵中,使其复原成为某个“加法表”(C=A⊙B),
则有如下性质:
1、A[0] = B[0] = 0
2、A[N-1] + B[N-1] = N*N-1
3、∑A + ∑B = N*(N*N-1)/2

若记N阶连续整数方阵可重排成“加法表”方阵的数目为 f(N),
则我还总结了如下规律:
1、f(p) = 1;
2、f(p^2) = 3;
3、f(p^3) = 10
4、f(p^4) = 35
5、一般地,f(p^r) = (2r-1)! / ((r-1)! * r!)
6、f(p1 * p2) = 7
7、f(p1^2 * p2) = 42
8、f(p1^3 * p2) = 230
9、f(p1 * p2 * p3) = 115
10、。。。
(其中 r 为任意自然数,p、p1、p2、p3、...均为素数,且 p1<p2<p3<...)
所以感觉与 N 分解质因数后的指数序列相关,而与具体的素因子无关。

to mathe,
“现在我找到一种算法,速度还是挺快的.”,愿闻其祥。
我原本以为仅仅是去除数组索引(如a[i]、ra[i]必=i);但进行这些优化后速度仍不理想。
不知你又采取了什么更好的算法?让我们共同学习学习。。。
gxqcn 2006-01-12
  • 打赏
  • 举报
回复
临近年终,部门有不少活动,所以回复迟了些。

运行了一下最新的代码,确实非常快!
代码写得非常精炼,算法也很巧妙(虽然我还不甚明了),
我原本说双休日写个程序验证自己的一个想法,现在看来没什么必要了。

前面曾提到“乘法表的复原”问题,不知 mathe 上述算法可否作适当修改而解决它,为方便验证程序,下面给出一个被打乱了(实际是被排序过了)的“乘法表”:
C[256] = { 11, 13, 17, 19, 22, 25, 26, 27, 29, 31, 33, 34, 37, 38, 39, 41, 43, 44, 47, 49, 50, 51, 52, 53, 54, 55, 57, 58, 61, 62, 65, 66, 68, 73, 74, 75, 76, 77, 78, 81, 82, 85, 86, 87, 88, 91, 93, 94, 95, 98, 99, 100, 102, 104, 106, 108, 110, 111, 114, 116, 117, 119, 121, 122, 123, 124, 125, 129, 130, 132, 133, 135, 136, 141, 143, 145, 146, 147, 148, 150, 152, 153, 154, 155, 156, 159, 162, 164, 165, 170, 171, 172, 174, 175, 182, 183, 185, 186, 187, 188, 189, 190, 195, 196, 200, 203, 204, 205, 209, 212, 215, 216, 217, 219, 220, 222, 225, 228, 231, 232, 235, 238, 243, 244, 245, 246, 248, 250, 255, 258, 259, 260, 261, 265, 266, 270, 273, 275, 279, 282, 285, 287, 290, 292, 294, 296, 297, 300, 301, 305, 310, 318, 319, 324, 328, 329, 333, 340, 341, 343, 344, 348, 350, 357, 365, 366, 369, 370, 371, 372, 375, 376, 378, 380, 387, 392, 399, 405, 406, 407, 410, 423, 424, 427, 430, 434, 435, 438, 441, 444, 451, 465, 470, 473, 477, 488, 490, 492, 500, 511, 516, 517, 518, 525, 530, 539, 540, 549, 555, 564, 567, 574, 580, 583, 584, 588, 602, 609, 610, 615, 620, 636, 645, 651, 657, 658, 671, 686, 705, 730, 732, 735, 740, 742, 777, 795, 803, 820, 854, 860, 861, 876, 903, 915, 940, 980, 987, 1022, 1029, 1060, 1095, 1113, 1220, 1281, 1460, 1533 };
请反求出整数数组:A[16]、B[16],使得 A[i]*B[j](0≤i,j≤15)遍历数组 C[].

希望 mathe 能再贴出“复原乘法表”的源代码,这样该类问题就得到了完美解决,本贴的最终目的就达到了。
mathe 2006-01-12
  • 打赏
  • 举报
回复
其实还有一个问题可以问:
对于任意一个N*N解方阵C,如果我们知道至少有一组分解方案:
C=A⊙B (当然要求C[0][0]=A[0]=B[0]=0)
那么是否分解方案的数目不多于0..N*N-1的分解方案数目?
也就是比如C是36*36阶的,那么分解方案数目必然不多于393
加载更多回复(42)

33,008

社区成员

发帖
与我相关
我的任务
社区描述
数据结构与算法相关内容讨论专区
社区管理员
  • 数据结构与算法社区
加入社区
  • 近7日
  • 近30日
  • 至今
社区公告
暂无公告

试试用AI创作助手写篇文章吧