看看稍微修改以后就可以求出s(7)
// shc.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <math.h>
#define L 40
#define N 7
#define K 5
unsigned n;
unsigned long long M;
unsigned long long R;
unsigned long long maxM;
int A;
long long count;
long long gc;
void solve()
{
unsigned long long V=M*M+R;
unsigned long long r=M%R;
r=R-r;
unsigned long long u=r;
for(;u<=M;u+=R){
if(V%u==0&&(u+M>=R*n)){
int i;
for(i=0;i<K;i++){
printf("%d ",n);
}
printf("%lld %lld\n",(u+M)/R,(V/u+M)/R);
}
}
}
void output(int level, int overflow)
{
if(level==K){
int a3=A%4;
int a1=A%4;
/* if(A==1&&a1==1)
return;
if(A==1&&a3==0)
return;
if(A==0&&A==0){
if((a1==1||a1==2)&&a3<2)
return;
}
*/ int c=(int)ceil(log10((double)(M)));
if(c>=L)c=L-1;
if(c<0)c=0;
gc++;
solve();
}
if(M>maxM)maxM=M;
count++;
}
unsigned long long gcd(unsigned long long x,unsigned long long y)
{
while(y!=0){
unsigned long long r=x%y;
x=y;
y=r;
}
return x;
}
void search(int level)
{
if(level>=K){
output(level,0);
return;
}
double ll=(double)M/R;
double uu=ll*(N-level);
if(uu>=4294967295.5||uu>=18446744073709551615.0/M){///overflow
output(level,1);
return;
}
unsigned li=(unsigned)ceil(ll);
unsigned ui=(unsigned)uu;
if(li<=n)li=n+1;
unsigned i;
for(i=li;i<=ui;i++){
if(gcd(M,i)!=1)
continue;
n=i;
M=M*i;
R=R*i-M;
A++;
search(level+1);
A--;
}
}
void search0()
{
int i;
for(i=0;i<4;i++)A=0;
for(i=2;i<=7;i++){
n=i;
R=i-1;
M=i;
A++;
search(1);
A--;
}
}
int _tmain(int argc, _TCHAR* argv[])
{
search0();
int i;
return 0;
}
根据45#,我们有340万个平均22位最多28位的整数需要分解因子.
我不知道分解这么大的整数平均需要多长时间.假设分解每一个数需要1秒钟,那么也只需要一个CPU大概1000小时,从计算量上看,是完全可行的.
而如果直接采用试除法,由于已经知道因子的特殊形式,平均每个大概需要试除数百万次,也就是大概需要$10^14$次128比特整数的除法运算,从计算量上看,可能比因子分解算法慢数倍,但是从计算量上也是可行的.而且对于这个方法,我们只要提供一个128比特整数比较有效的加减乘除代码就可以了.
谁先提供一个效率比较高的128比特整数的加减乘除运算代码看看,我们只需要修改52#中的solve程序采用128比特计算,然后修改N为8,K为6就可以了.可以先部分运行一下代码测试一下大概的速度
快速乘法有了:
x86上128位二进制乘法最快速算法征解
除法有没有好的?
看了一下solve的代码,只需要128位除64位的除法运算,那么应该直接用汇编指令就可以了.
谁将上面的solve函数改写一下?汇编我长时间不用,已经非常不熟练了
除法要快速，一般是针对除数是固定值时，预先算好一些算法值以加速运算。
对于128位除64位的除法，用一般的试商法即可；
如果在64位系统下，只要将高地位分别移进 RDX:RAX，而后 DIV 除数 即可。
存在一些很大的M/R,所以还是要用对V因子分解的方法
然后对V每个小于sqrt(V)的因子d,如果(d+M)是R倍数,就可以构着一个解.
这个还是谁用Mathematica来做比较好.如果能够运行两个月左右估计就能解决这个问题了
呵,这个还是谁用Mathematica来做比较好.如果能够运行两个月左右估计就能解决这个问题了
这个难题可能对wayne来说不是问题,可惜我的mathematica编程不行.....
请教mathe:
51#楼的程序文件格式为.cpp吧?
在visualstudio 2008中运行怎么不成功呢?
见下图:
看图示似乎是编译未通过，尚未生成 exe 文件。
你可以用向导产生一个工程（目的：生成stdafx.h及staafx.cpp），
然后再将 mathe 的代码粘贴进来进行编译。
试着用LiLDA求解,发现代码有bug,做因子分解时经常崩溃掉,而且稍微添加一点代码,崩溃的地方都不同,不知道什么原因
#include <time.h>
#include <LiDIA/rational_factorization.h>
#include <iostream.h>
using namespace LiDIA;
using namespace std;
#include <math.h>
#define L 40
#define N 8
#define K 6
unsigned long n;
bigint M;
bigint R;
bigint maxM;
int A;
long long c;
long long gc;
int start_time;
#define MAX_INDEX 20
int iindex;
void solve()
{
bigint V,MM,RR,uu;
MM=M;
RR=R;
V=MM*MM+RR;
rational_factorization f;
f.assign(V);
f.factor();
if(!f.is_prime_factorization()){
printf("cannot factor %lld solution\n",c);
}else{
if(f.no_of_comp()>=MAX_INDEX){
printf("too much compondents for %lld solution\n",c);
}else{
int i;
for(i=0;i<f.no_of_comp();i++){
iindex=0;
}
do{
bigint v=1;
for(i=0;i<f.no_of_comp();i++){
int j;
for(j=0;j<iindex;j++){
v*=f.base(i);
}
}
if((v+MM)%RR==0){///solution found
for(i=0;i<K;i++){
cout<<"n"<<i+1<<"="<<n<<";";
}
cout<<"n"<<K+1<<'='<<(v+MM)/RR<<";n"<<K+2<<"="<<(V/v+MM)/RR<<";\n";
}
for(i=f.no_of_comp()-1;i>=0;i--){
if(iindex==f.exponent(i)){
iindex=0;
continue;
}else{
iindex++;
break;
}
}
if(i<0)break;
}while(1);
}
}
fflush(stdout);
}
void output(int level, int overflow)
{
if(level==K){
solve();
}
if(M>maxM)maxM=M;
c++;
fprintf(stderr,"processed %lld by %ds\n",c,time(NULL)-start_time);
}
void search(int level)
{
if(level>=K){
output(level,0);
return;
}
double ll=(double)dbl(M)/dbl(R);
double uu=ll*(N-level);
if(uu>=4294967295.5){///overflow
output(level,1);
return;
}
unsigned long li=(unsigned long)ceil(ll);
unsigned long ui=(unsigned long)uu;
if(li<=n)li=n+1;
unsigned long i;
for(i=li;i<=ui;i++){
if(gcd(M,i)!=1)
continue;
n=i;
M=M*i;
R=R*i-M;
A++;
search(level+1);
A--;
}
}
void search0()
{
int i;
for(i=0;i<4;i++)A=0;
for(i=2;i<=7;i++){
n=i;
R=i-1;
M=i;
A++;
search(1);
A--;
}
}
int main()
{
start_time=time(NULL);
search0();
int i;
return 0;
}