这道题大约没答案

.·.·.

我们已知一系列的正整数$n_i,m_i$(i=1,2,...,k)
对集合$J\subset\{1,2,...,k\}$
记$n_J$为，满足对$J$中任意元素$i$，$n_J$除以$m_i$的余数为$n_i$的最小的正整数
再记$f(J)=\frac{n_J}{\prod_{i\in J}m_i}$
问有没有一个能够快速求出J使得f(J)最小的算法呢？


感觉不像有的样子
不过这个题目有一点好玩的地方在于，如果只是给定$m_i$，而$n_i$都是随机选取的，那么我们能否求出f(J)最小值的期望呢？
或者，如果我们能够构造一系列$J_s$使得$f(J_s)$达到最小，那么，我们能构造出来的f(J)究竟有多小呢？

一个栗子，
n1=3,n2=4,n3=5,m1=7,m2=9,m3=11
这时候J至多有7种取值，而显然J={1}的时候，我们得到的f(J)=3/7达到最小

而对
n1=1,n2=1,n3=1,m1=7,m2=9,m3=11
这时候J至多有7种取值，而显然J={1,2,3}的时候，我们得到的f(J)=1/693达到最小

.·.·.

好像应该补充一个估计的
对于这一点：
如果只是给定$m_i$，而$n_i$都是随机选取的，那么我们能否求出$f(J)$最小值的期望呢？
刚刚给出了一个很宽的上下界
首先，注意到每一个$n_i$都是随机选取的，我们可以近似认为$\frac{n_i}{m_i}$服从均匀分布，于是期望的上界一定会小于等于均匀分布的最小的次序统计量的期望，也就是$E Min_J f(J)<\frac 1 {k+1}$
下界，注意到一共$2^k$种组合，假定各个组合之间彼此独立的话，我们应该可以找到$2^k$个均匀分布，从而$E Min_J f(J)<\frac 1 {2^k+1}$

估计有点粗……在想有没有什么更精确的估计

补充内容 (2018-6-24 19:22):
下界的那个不等号是大于……当时写错了。
以及，玩了玩这个，忽然感觉剩余系好有意思……虽然不是随机数但真的跟随机数一样跳得厉害

hujunhua

转述一下，不知道我理解的对不对。就是

给定一个同余方程集
           \(x\equiv n_i\pmod{m_i}\cdots\cdots(i)\)
              \(i\in I=\{1,2,\cdots,k\}\)
取其中的若干个元素构成的方程组
           \(x\equiv n_i\pmod{m_i}\cdots\cdots(i)\)
              \(i\in J\subset{I}\)
的最小正解`n_J`和模组积`m_J`(`=\D\prod_{i\in{J}}{m_i}`), 定义
              \(f(J)=\frac{n_J}{m_J}\)
要求编程计算怎样选择`J`可得到最小的`f(J)`.

.·.·.

hujunhua 发表于 2018-6-27 20:14
转述一下，不知道我理解的对不对。就是

给定一个同余方程集

是的，以及
问题甚至可以简化成
每一个m_i都是某一个质数的整数次幂
反正中国剩余定理就是干这个活儿的……
以及，暂时的发现是，最终结果（在随机意义下）可能远比我之前想象得要小。
借助R，python以及pari/gp，简单实验了一下
python写的pari/gp的代码生成程序（为了防止排版问题我把空格替换成了!&^%!，只要替换回去就可以送python3执行了）：

1. str="searched(a,b,c,d,e,f)=vecmin(["
   
2. modlist=[3,5,7,11,13,17]
   
3. for!&^%!i!&^%!in!&^%!range(1,64):
   
4. !&^%!!&^%!if!&^%!bin(i).count('1')==1:
   
5. !&^%!!&^%!!&^%!!&^%!!&^%!j=len(bin(i))-bin(i).index('1')-1
   
6. !&^%!!&^%!!&^%!!&^%!!&^%!str+="lift(%c)/%d,"%(j+97,modlist[j])
   
7. !&^%!!&^%!!&^%!!&^%!!&^%!print([i,j,k])
   
8. !&^%!!&^%!else:
   
9. !&^%!!&^%!!&^%!g=[j!&^%!for!&^%!j!&^%!in!&^%!bin(i)[2:]];g.reverse()
   
10. !&^%!!&^%!!&^%!s="lift(chinese(["
    
11. !&^%!!&^%!!&^%!mul=1
    
12. !&^%!!&^%!!&^%!for!&^%!j,k!&^%!in!&^%!enumerate(g):
    
13. !&^%!!&^%!!&^%!!&^%!if!&^%!k=='1':
    
14. !&^%!!&^%!!&^%!!&^%!!&^%!s+="%c,"%(j+97)
    
15. !&^%!!&^%!!&^%!!&^%!!&^%!mul*=modlist[j]
    
16. !&^%!!&^%!!&^%!str+=s[:-1]+"]))/%d,"%mul
    
17. 
    
18. print(str[:-1]+"])")
    

复制代码

这是机器生成的难看的gp代码：

1. searched(a,b,c,d,e,f)=vecmin([lift(a)/3,lift(b)/5,lift(chinese([a,b]))/15,lift(c)/7,lift(chinese([a,c]))/21,lift(chinese([b,c]))/35,lift(chinese([a,b,c]))/105,lift(d)/11,lift(chinese([a,d]))/33,lift(chinese([b,d]))/55,lift(chinese([a,b,d]))/165,lift(chinese([c,d]))/77,lift(chinese([a,c,d]))/231,lift(chinese([b,c,d]))/385,lift(chinese([a,b,c,d]))/1155,lift(e)/13,lift(chinese([a,e]))/39,lift(chinese([b,e]))/65,lift(chinese([a,b,e]))/195,lift(chinese([c,e]))/91,lift(chinese([a,c,e]))/273,lift(chinese([b,c,e]))/455,lift(chinese([a,b,c,e]))/1365,lift(chinese([d,e]))/143,lift(chinese([a,d,e]))/429,lift(chinese([b,d,e]))/715,lift(chinese([a,b,d,e]))/2145,lift(chinese([c,d,e]))/1001,lift(chinese([a,c,d,e]))/3003,lift(chinese([b,c,d,e]))/5005,lift(chinese([a,b,c,d,e]))/15015,lift(f)/17,lift(chinese([a,f]))/51,lift(chinese([b,f]))/85,lift(chinese([a,b,f]))/255,lift(chinese([c,f]))/119,lift(chinese([a,c,f]))/357,lift(chinese([b,c,f]))/595,lift(chinese([a,b,c,f]))/1785,lift(chinese([d,f]))/187,lift(chinese([a,d,f]))/561,lift(chinese([b,d,f]))/935,lift(chinese([a,b,d,f]))/2805,lift(chinese([c,d,f]))/1309,lift(chinese([a,c,d,f]))/3927,lift(chinese([b,c,d,f]))/6545,lift(chinese([a,b,c,d,f]))/19635,lift(chinese([e,f]))/221,lift(chinese([a,e,f]))/663,lift(chinese([b,e,f]))/1105,lift(chinese([a,b,e,f]))/3315,lift(chinese([c,e,f]))/1547,lift(chinese([a,c,e,f]))/4641,lift(chinese([b,c,e,f]))/7735,lift(chinese([a,b,c,e,f]))/23205,lift(chinese([d,e,f]))/2431,lift(chinese([a,d,e,f]))/7293,lift(chinese([b,d,e,f]))/12155,lift(chinese([a,b,d,e,f]))/36465,lift(chinese([c,d,e,f]))/17017,lift(chinese([a,c,d,e,f]))/51051,lift(chinese([b,c,d,e,f]))/85085,lift(chinese([a,b,c,d,e,f]))/255255])
   
2. c=0;for(i3=1,2,for(i5=1,4,for(i7=1,6,for(i11=1,10,for(i13=1,12,for(i17=1,16,c++))))));c
   
3. a=Col(0,92160);
   
4. cc=0;for(i3=1,2,for(i5=1,4,for(i7=1,6,for(i11=1,10,for(i13=1,12,for(i17=1,16,cc++;a[cc]=(searched(Mod(i3,3),Mod(i5,5),Mod(i7,7),Mod(i11,11),Mod(i13,13),Mod(i17,17)))))))))
   
5. a
   
6. \w /mnt/b/output.txt
   

复制代码

这里是R的绘图部分结果

1. vec=source("B:\\output.txt")#这里更改了output使得其符合R的语法，然而仍然读了15秒……
   
2. runif_count=48;len=length(vec);g=apply(matrix(runif(len*runif_count),len),1,min);plot(sort(vec),sort(g),type='l',main="g=apply(matrix(runif(len*runif_count=48),len),1,min)");abline(0,1,col=2)
   

复制代码

对6个质数（3,5,7,11,13,17）的测试结果表明，我们的目标值大约跟$U_(48,1)$的分布差不多，这里$U_(j,i)$代表的是，$j$个服从均匀分布$U(0,1)$的随机变量之中第$i$小的随机变量
当然这里只是头部如此，尾巴上面会有偏差，毕竟我们选择的都是小质数，如果质数大一些，或许这个出现在尾部的偏差理应会减小……
所以感觉是，最优值就在那里，但我们就是找不到……

.·.·.

这是用了两个小时才读到R里面的数据
用的质数是5,7,11,13,17,19，可以看出如果质数足够大，分布或许真的会接近估计出来的那个下界……
所以问题只是，最小值在哪里取到

.·.·.

小发现
可以借助CRT复原一个大数N，使得$N=n_i mod m_i$
之后连分数计算$\frac N{\prod_i m_i}$
如果我们得到的连分数之中，某一个渐进分数的分母恰好可以被$m_i$因子分解掉，那么去掉这些$m_i$，选取剩下的那些$m_i$，那么我们有一定概率得到一个相当小的结果。