关于单位分数的一些难题

mathe

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

2. // shc.cpp : Defines the entry point for the console application.
3. //
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];
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. }
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. }
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. }
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. }
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. }
110. int _tmain(int argc, _TCHAR* argv[])
111. {
112. search0();
113. int i;
114. return 0;
115. }

复制代码

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,做因子分解时经常崩溃掉,而且稍微添加一点代码,崩溃的地方都不同,不知道什么原因

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;
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. }
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. }
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. }
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. }
123. int main()
124. {
125. start_time=time(NULL);
126. search0();
127. int i;
128. return 0;
129. }

复制代码