关于单位分数的一些难题

mathe

看看稍微修改以后就可以求出s(7)

1. 
   
2. // shc.cpp : Defines the entry point for the console application.
   
3. //
   
4. 
   
5. #include "stdafx.h"
   
6. #include <math.h>
   
7. #define L 40
   
8. #define N 7
   
9. #define K 5
   
10. unsigned n[K];
    
11. unsigned long long M[K];
    
12. unsigned long long R[K];
    
13. unsigned long long maxM;
    
14. int A[4];
    
15. long long count[K];
    
16. long long gc[L];
    
17. 
    
18. void solve()
    
19. {
    
20.     unsigned long long V=M[K-1]*M[K-1]+R[K-1];
    
21.     unsigned long long r=M[K-1]%R[K-1];
    
22.     r=R[K-1]-r;
    
23.     unsigned long long u=r;
    
24.     for(;u<=M[K-1];u+=R[K-1]){
    
25.         if(V%u==0&&(u+M[K-1]>=R[K-1]*n[K-1])){
    
26.             int i;
    
27.             for(i=0;i<K;i++){
    
28.                 printf("%d ",n[i]);
    
29.             }
    
30.             printf("%lld %lld\n",(u+M[K-1])/R[K-1],(V/u+M[K-1])/R[K-1]);
    
31.         }
    
32.     }
    
33. }
    
34. 
    
35. void output(int level, int overflow)
    
36. {
    
37.     if(level==K){
    
38.         int a3=A[3]%4;
    
39.         int a1=A[1]%4;
    
40. /*        if(A[2]==1&&a1==1)
    
41.             return;
    
42.         if(A[0]==1&&a3==0)
    
43.             return;
    
44.         if(A[0]==0&&A[2]==0){
    
45.             if((a1==1||a1==2)&&a3<2)
    
46.                 return;
    
47.         }
    
48.   */      int c=(int)ceil(log10((double)(M[level-1])));
    
49.         if(c>=L)c=L-1;
    
50.         if(c<0)c=0;
    
51.         gc[c]++;
    
52.         solve();
    
53.     }
    
54.     if(M[level-1]>maxM)maxM=M[level-1];
    
55.     count[level-1]++;
    
56. }
    
57. 
    
58. unsigned long long gcd(unsigned long long x,unsigned long long y)
    
59. {
    
60.     while(y!=0){
    
61.         unsigned long long r=x%y;
    
62.         x=y;
    
63.         y=r;
    
64.     }
    
65.     return x;
    
66. }
    
67. 
    
68. void search(int level)
    
69. {
    
70.     if(level>=K){
    
71.         output(level,0);
    
72.         return;
    
73.     }
    
74.     double ll=(double)M[level-1]/R[level-1];
    
75.     double uu=ll*(N-level);
    
76.     if(uu>=4294967295.5||uu>=18446744073709551615.0/M[level-1]){///overflow
    
77.         output(level,1);
    
78.         return;
    
79.     }
    
80.     unsigned li=(unsigned)ceil(ll);
    
81.     unsigned ui=(unsigned)uu;
    
82.     if(li<=n[level-1])li=n[level-1]+1;
    
83.     unsigned i;
    
84.     for(i=li;i<=ui;i++){
    
85.         if(gcd(M[level-1],i)!=1)
    
86.             continue;
    
87.         n[level]=i;
    
88.         M[level]=M[level-1]*i;
    
89.         R[level]=R[level-1]*i-M[level-1];
    
90.         A[i%4]++;
    
91.         search(level+1);
    
92.         A[i%4]--;
    
93.     }
    
94. }
    
95. 
    
96. void search0()
    
97. {
    
98.     int i;
    
99.     for(i=0;i<4;i++)A[i]=0;
    
100.     for(i=2;i<=7;i++){
     
101.         n[0]=i;
     
102.         R[0]=i-1;
     
103.         M[0]=i;
     
104.         A[i%4]++;
     
105.         search(1);
     
106.         A[i%4]--;
     
107.     }
     
108. }
     
109. 
     
110. int _tmain(int argc, _TCHAR* argv[])
     
111. {
     
112.     search0();
     
113.     int i;
     
114.         return 0;
     
115. }
     
116. 
     

复制代码

mathe

根据45#,我们有340万个平均22位最多28位的整数需要分解因子.
我不知道分解这么大的整数平均需要多长时间.假设分解每一个数需要1秒钟,那么也只需要一个CPU大概1000小时,从计算量上看,是完全可行的.
而如果直接采用试除法,由于已经知道因子的特殊形式,平均每个大概需要试除数百万次,也就是大概需要$10^14$次128比特整数的除法运算,从计算量上看,可能比因子分解算法慢数倍,但是从计算量上也是可行的.而且对于这个方法,我们只要提供一个128比特整数比较有效的加减乘除代码就可以了.
谁先提供一个效率比较高的128比特整数的加减乘除运算代码看看,我们只需要修改52#中的solve程序采用128比特计算,然后修改N为8,K为6就可以了.可以先部分运行一下代码测试一下大概的速度

mathe

快速乘法有了:
x86上128位二进制乘法最快速算法征解
除法有没有好的?

mathe

看了一下solve的代码,只需要128位除64位的除法运算,那么应该直接用汇编指令就可以了.
谁将上面的solve函数改写一下?汇编我长时间不用,已经非常不熟练了

gxqcn

除法要快速，一般是针对除数是固定值时，预先算好一些算法值以加速运算。
对于128位除64位的除法，用一般的试商法即可；
如果在64位系统下，只要将高地位分别移进 RDX:RAX，而后 DIV 除数 即可。

mathe

存在一些很大的M/R,所以还是要用对V因子分解的方法
然后对V每个小于sqrt(V)的因子d,如果(d+M)是R倍数,就可以构着一个解.
这个还是谁用Mathematica来做比较好.如果能够运行两个月左右估计就能解决这个问题了

数学星空

呵,这个还是谁用Mathematica来做比较好.如果能够运行两个月左右估计就能解决这个问题了
这个难题可能对wayne来说不是问题,可惜我的mathematica编程不行.....

数学星空

请教mathe:
51#楼的程序文件格式为.cpp吧?
在visual  studio 2008中运行怎么不成功呢?
见下图:

gxqcn

看图示似乎是编译未通过，尚未生成 exe 文件。

你可以用向导产生一个工程（目的：生成stdafx.h及staafx.cpp），
然后再将 mathe 的代码粘贴进来进行编译。

mathe

试着用LiLDA求解,发现代码有bug,做因子分解时经常崩溃掉,而且稍微添加一点代码,崩溃的地方都不同,不知道什么原因

1. 
   
2. #include <time.h>
   
3. #include <LiDIA/rational_factorization.h>
   
4. #include <iostream.h>
   
5. using namespace LiDIA;
   
6. using namespace std;
   
7. #include <math.h>
   
8. #define L 40
   
9. #define N 8
   
10. #define K 6
    
11. unsigned long n[K];
    
12. bigint M[K];
    
13. bigint R[K];
    
14. bigint maxM;
    
15. int A[4];
    
16. long long c;
    
17. long long gc[L];
    
18. int start_time;
    
19. 
    
20. #define MAX_INDEX 20
    
21. int iindex[MAX_INDEX];
    
22. void solve()
    
23. {
    
24.         bigint V,MM,RR,uu;
    
25.         MM=M[K-1];
    
26.         RR=R[K-1];
    
27.         V=MM*MM+RR;
    
28.         rational_factorization f;
    
29.         f.assign(V);
    
30.         f.factor();
    
31.         if(!f.is_prime_factorization()){
    
32.                 printf("cannot factor %lld solution\n",c);
    
33.         }else{
    
34.                 if(f.no_of_comp()>=MAX_INDEX){
    
35.                         printf("too much compondents for %lld solution\n",c);
    
36.                 }else{
    
37.                         int i;
    
38.                         for(i=0;i<f.no_of_comp();i++){
    
39.                                 iindex[i]=0;
    
40.                         }
    
41.                         do{
    
42.                                 bigint v=1;
    
43.                                 for(i=0;i<f.no_of_comp();i++){
    
44.                                         int j;
    
45.                                         for(j=0;j<iindex[i];j++){
    
46.                                                 v*=f.base(i);
    
47.                                         }
    
48.                                 }
    
49.                                 if((v+MM)%RR==0){///solution found
    
50.                                         for(i=0;i<K;i++){
    
51.                                                 cout<<"n"<<i+1<<"="<<n[i]<<";";
    
52.                                         }
    
53.                                         cout<<"n"<<K+1<<'='<<(v+MM)/RR<<";n"<<K+2<<"="<<(V/v+MM)/RR<<";\n";
    
54.                                 }
    
55.                                 for(i=f.no_of_comp()-1;i>=0;i--){
    
56.                                         if(iindex[i]==f.exponent(i)){
    
57.                                                 iindex[i]=0;
    
58.                                                 continue;
    
59.                                         }else{
    
60.                                                 iindex[i]++;
    
61.                                                 break;
    
62.                                         }
    
63.                                 }
    
64.                                 if(i<0)break;
    
65.                         }while(1);
    
66.                 }
    
67.         }
    
68.         fflush(stdout);
    
69. }
    
70. 
    
71. void output(int level, int overflow)
    
72. {
    
73.     if(level==K){
    
74.         solve();
    
75.     }
    
76.     if(M[level-1]>maxM)maxM=M[level-1];
    
77.     c++;
    
78.     fprintf(stderr,"processed %lld by %ds\n",c,time(NULL)-start_time);
    
79. }
    
80. 
    
81. void search(int level)
    
82. {
    
83.     if(level>=K){
    
84.         output(level,0);
    
85.         return;
    
86.     }
    
87.     double ll=(double)dbl(M[level-1])/dbl(R[level-1]);
    
88.     double uu=ll*(N-level);
    
89.     if(uu>=4294967295.5){///overflow
    
90.         output(level,1);
    
91.         return;
    
92.     }
    
93.     unsigned long li=(unsigned long)ceil(ll);
    
94.     unsigned long ui=(unsigned long)uu;
    
95.     if(li<=n[level-1])li=n[level-1]+1;
    
96.     unsigned long i;
    
97.     for(i=li;i<=ui;i++){
    
98.         if(gcd(M[level-1],i)!=1)
    
99.             continue;
    
100.         n[level]=i;
     
101.         M[level]=M[level-1]*i;
     
102.         R[level]=R[level-1]*i-M[level-1];
     
103.         A[i%4]++;
     
104.         search(level+1);
     
105.         A[i%4]--;
     
106.     }
     
107. }
     
108. 
     
109. void search0()
     
110. {
     
111.     int i;
     
112.     for(i=0;i<4;i++)A[i]=0;
     
113.     for(i=2;i<=7;i++){
     
114.         n[0]=i;
     
115.         R[0]=i-1;
     
116.         M[0]=i;
     
117.         A[i%4]++;
     
118.         search(1);
     
119.         A[i%4]--;
     
120.     }
     
121. }
     
122. 
     
123. int main()
     
124. {
     
125.         start_time=time(NULL);
     
126.     search0();
     
127.     int i;
     
128.         return 0;
     
129. }
     
130. 
     

复制代码