在C语言中用蛮力反驳Euler命题

#1 楼

.L6:
        add     rbx, 1
        cmp     rbx, 500000
        je      .L18
.L8:
        mov     rax, QWORD PTR .LC0[rip]
        movsd   xmm0, QWORD PTR [rsp+40]
        movq    xmm1, rax
        call    pow                        ; (1)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+8], xmm0
        movsd   xmm0, QWORD PTR [rsp+48]
        movq    xmm1, rax
        call    pow                        ; (2)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+16], xmm0
        movsd   xmm0, QWORD PTR [rsp+32]
        movq    xmm1, rax
        call    pow                        ; (3)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+24], xmm0
        pxor    xmm0, xmm0
        cvtsi2sd        xmm0, rbx
        movq    xmm1, rax
        call    pow                        ; (4)
        movsd   xmm2, QWORD PTR [rsp+8]
        addsd   xmm2, QWORD PTR [rsp+16]
        movapd  xmm1, xmm0
        movsd   xmm0, QWORD PTR [rsp+24]
        addsd   xmm0, xmm2
        ucomisd xmm0, xmm1
        jp      .L6
        jne     .L6

long int pow4(long int x){
    return x * x * x * x;
}

movq    %rdi, %rax
imulq   %rdi, %rax
imulq   %rax, %rax
ret

long int pow4(long int x){
    const long int t = x * x;
    return t * t;
}

int prop(long int A, long int B, long int C, long int D) {
  return (pow4(A) + pow4(B) + pow4(C) == pow4(D));
}

a = 1, b = 1, c = 1000, time = 114.156000s

a = 1, b = 1, c = 1000, time = 0.296000s
a = 1, b = 1, c = 2000, time = 0.578000s
a = 1, b = 1, c = 3000, time = 0.859000s
a = 1, b = 1, c = 4000, time = 1.140000s
a = 1, b = 1, c = 5000, time = 1.421000s
a = 1, b = 1, c = 6000, time = 1.703000s
a = 1, b = 1, c = 7000, time = 1.984000s
a = 1, b = 1, c = 8000, time = 2.265000s
a = 1, b = 1, c = 9000, time = 2.546000s
a = 1, b = 1, c = 10000, time = 2.828000s
a = 1, b = 1, c = 11000, time = 3.109000s
a = 1, b = 1, c = 12000, time = 3.390000s
a = 1, b = 1, c = 13000, time = 3.687000s
a = 1, b = 1, c = 14000, time = 3.968000s
a = 1, b = 1, c = 15000, time = 4.250000s
a = 1, b = 1, c = 16000, time = 4.531000s

…
a = 1, b = 1, c = 481000, time = 0.031000s
a = 1, b = 1, c = 482000, time = 0.031000s
a = 1, b = 1, c = 483000, time = 0.031000s
a = 1, b = 1, c = 484000, time = 0.031000s
a = 1, b = 1, c = 485000, time = 0.031000s
a = 1, b = 1, c = 486000, time = 0.031000s
a = 1, b = 1, c = 487000, time = 0.031000s
a = 1, b = 1, c = 488000, time = 0.031000s
a = 1, b = 1, c = 489000, time = 0.031000s
a = 1, b = 1, c = 490000, time = 0.031000s
a = 1, b = 1, c = 491000, time = 0.031000s
a = 1, b = 1, c = 492000, time = 0.031000s
a = 1, b = 1, c = 493000, time = 0.031000s
a = 1, b = 1, c = 494000, time = 0.031000s
a = 1, b = 1, c = 495000, time = 0.031000s
a = 1, b = 1, c = 496000, time = 0.031000s
a = 1, b = 1, c = 497000, time = 0.031000s
a = 1, b = 1, c = 498000, time = 0.031000s
a = 1, b = 1, c = 499000, time = 0.031000s
a = 1, b = 1, time = 0.031000s

#include <stdio.h>
#include <time.h>
#include <math.h>

#2 楼

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {

for (a = 1; a < 100000; a++) {
    for (b = a; b < 300000; b++) {
        for (c = b; c < 500000; c++) {

#3 楼

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <queue>
#include <vector>
using std::priority_queue;
using std::vector;

//typedef long long int bigint;
typedef __int128 bigint;

// Look for a^4+b^4+c^4 = d^4 by comparing
// ascending list of a^4+b^4 with ascending list of d^4-c^4

vector<bigint> list;
vector<int> diffptr,sumptr;

bool sumcmp(int a0,int a1){
  int b0=sumptr[a0];
  int b1=sumptr[a1];
  return list[a0]+list[b0]>list[a1]+list[b1];
}

bool diffcmp(int c0,int c1){
  int d0=diffptr[c0];
  int d1=diffptr[c1];
  return list[d0]-list[c0]>list[d1]-list[c1];
}

int main(int ac,char**av){
  int a,c,i,n;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint sumval,diffval;

  list.resize(n);
  diffptr.resize(n);
  sumptr.resize(n);

  priority_queue<int,vector<int>,decltype(&sumcmp)> sumnext(sumcmp);
  priority_queue<int,vector<int>,decltype(&diffcmp)> diffnext(diffcmp);

  for(i=0;i<n;i++)list[i]=bigint(i)*i*i*i;
  for(a=1;a<n;a++){sumptr[a]=a;sumnext.push(a);}
  for(c=0;c<n-1;c++){diffptr[c]=c+1;diffnext.push(c);}

  while(!sumnext.empty()&&!diffnext.empty()){
    a=sumnext.top();sumval=list[sumptr[a]]+list[a];
    c=diffnext.top();diffval=list[diffptr[c]]-list[c];
    if(sumval==diffval){printf("BINGO %d^4+%d^4+%d^4=%d^4\n",a,sumptr[a],c,diffptr[c]);fflush(stdout);}
    if(sumval<=diffval){
      sumnext.pop();
      sumptr[a]++;if(sumptr[a]<n)sumnext.push(a);
    }else{
      diffnext.pop();
      diffptr[c]++;if(diffptr[c]<n)diffnext.push(c);
    }
  }

}

#4 楼

#5 楼

#6 楼

#7 楼

#8 楼

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4).

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

vector<bigint> pow4;
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,s0,s1;
  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=bigint(min(5.*n,MAXMEM/26.));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){
    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)

    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=1;a<n&&pow4[a]<ub;a++){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=min(int(sqrt(sqrt(ub-pow4[a]))),n);
      for(b=bmin;b<bmax;b++)list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    for(c=0;c<n-1;c++){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(d=dmin;d<dmax;d++)list1.push_back(pow4[d]-pow4[c]);
    }
    std::sort(list1.begin(),list1.end());

    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }
  }
}

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

#define MAXMEM 4e9 // Allow the use of up to (approx) this memory in bytes

vector<bigint> pow4;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,r,s0,s1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  pow4.resize(n);

  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint((MAXMEM-2*n*sizeof(int))*2.8));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));
  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(chunk/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.
// Parallelise over [A_i,B_i).

// Compile with -pthread flag

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
#include <thread>
#include <mutex>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

int n;// Maximum value of d
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
const int nthreads=12;
bigint chunk,chunksize;
vector<bigint> pow4;
std::mutex chunkserver,printer;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  std::lock_guard<std::mutex> lock(printer);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

bigint nextchunk(){// Returns next chunk, or -1 if none left
  std::lock_guard<std::mutex> lock(chunkserver);
  if(chunk>=bigint(n)*n)return -1;
  bigint chunkcopy=chunk;
  chunk+=chunksize;
  return chunkcopy;
}

void searcher(){
  int a,b,c,d,i,j,r,s0,s1;
  bigint ch;
  vector<bigint> list0,list1;
  while((ch=nextchunk())>=0){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=ch*ch;
    bigint ub=min(ch+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(ch/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

int main(int ac,char**av){
  int i,t;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint(MAXMEM*2.8/nthreads));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  vector<std::thread> thr(nthreads);
  for(t=0;t<nthreads;t++)thr[t]=std::thread(searcher);
  for(t=0;t<nthreads;t++)thr[t].join();
}

#9 楼

unsigned sq = 1;
for (unsigned x=1;; x++) {
  sq += x + x + 1;
  printf("%u\n", sq);
}

unsigned sq = 1;
unsigned cube = 1;
for (unsigned x=1;; x++) {
  cube += (sq + sq + sq) + (x + x + x) + 1;
  sq += x + x + 1;
  printf("%u\n", cube);
}

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

// Adjust to some 128 bit type
// The below is a place-holder for a true 128-bit type
typedef uintmax_t uint_fast128t;

typedef struct {
  uint_fast128t x4;  // 4th power of the value
  uint_fast128t x3;  // cube of the value
  uint_fast64_t x2;  // square of the value
  uint_fast32_t x1;  // location of OP's value a,b,c or d
                     // Assume a,b,c,d are 20 bits or less.
}quad;

#define quad_next(x) { \
  x.x4 += (uint_fast128t) x.x3*4 + x.x2*6 + 4*x.x1 + 1; \
  x.x3 += (uint_fast128t) x.x2*3 + x.x1*3 + 1; \
  x.x2 += (uint_fast64_t) x.x1*2 + 1; \
  x.x1++; \
  }

#define factor 20
uint_fast32_t a_max = 95800 / factor;
uint_fast32_t b_max = 217519 / factor;
uint_fast32_t c_max = 414560 / factor;
uint_fast32_t d_max = 422481 / factor;

void Elkies_B(quad a, quad b) {
  quad c = b;
  quad d = c;
  c.x4 += b.x4;

  // Could use a more efficient step here
  while (d.x4 < c.x4) {
    quad_next(d);
  }
  while (c.x1 <= c_max) {
    while (d.x4 < c.x4) {
      quad_next(d);
    }
    if (d.x4 == c.x4) {
      // Success! Print a,b,c,d
      printf("+ a:%lu b:%lu c:%lu d:%lu\n", //
              (unsigned long) a.x1, (unsigned long) b.x1, //
              (unsigned long) c.x1, (unsigned long) d.x1);
      exit(0);
    }
    quad_next(c);
  }
}

void Elkies_A(quad a) {
  quad b = a;
  b.x4 += a.x4;
  while (b.x1 <= b_max) {
    Elkies_B(a, b);
    quad_next(b);
  }
}

void Elkies(void) {
  quad a = { 1, 1, 1, 1 };
  while (a.x1 <= a_max) {
    printf("- a:%lu %llu %ju %ju\n", //
            (unsigned long) a.x1, (unsigned long long) a.x2, //
            a.x3, a.x4);
    fflush( stdout);
    Elkies_A(a);
    quad_next(a);
  }
  printf("Failed\n");
}

int main() {
  Elkies();
}

#10 楼

#11 楼

C:&gt euler413.exe
I haven't found unfinished work!
Please give R parameter, the Range will be R*10240000: 10
Do you want to start an exhaustive search in this Range or
only to search for special solutions?
( y=full search, n=special search ) n
Give the start value of a0. It should be 0&lt=start_a0&lt=16384: 0
Give the end value of a0. It should be start_a0&lt=end_a0&lt=16384: 16384
Building up some tables
Done
Testing: a0=1
Finished: a0=1,Range=102400000,type=0,Time: 0h0m0s,Date: Fri Nov 04 18:20:33 2016
Testing: a0=9
Finished: a0=9,Range=102400000,type=0,Time: 0h0m1s,Date: Fri Nov 04 18:20:34 2016
....

euler413work.txt
stat_euler(4,3,1).txt
results_euler(4,1,3).txt

Solution found! 98438073^4=50681927^4+96592480^4+22321400^4
Solution found! 17321721^4=8918279^4+16996960^4+3927800^4
Solution found! 40980657^4=21099343^4+40212320^4+9292600^4
Solution found! 20615673^4=15365639^4+18796760^4+2682440^4  (primitive)
Solution found! 20615673^4=15365639^4+2682440^4+18796760^4  (primitive)
Solution found! 7182177^4=3697823^4+7047520^4+1628600^4
Solution found! 30841113^4=15878887^4+30262880^4+6993400^4
Solution found! 101817921^4=52422079^4+99908960^4+23087800^4
Solution found! 12197457^4=8282543^4+11289040^4+5870000^4  (primitive)
Solution found! 34220961^4=17619039^4+33579360^4+7759800^4
Solution found! 57879897^4=29800103^4+56794720^4+13124600^4
Solution found! 81538833^4=41981167^4+80010080^4+18489400^4
Solution found! 16003017^4=4479031^4+14173720^4+12552200^4  (primitive)
Solution found! 422481^4=217519^4+414560^4+95800^4  (primitive)
Solution found! 16430513^4=16281009^4+3642840^4+7028600^4  (primitive)
Solution found! 47740353^4=24579647^4+46845280^4+10825400^4
Solution found! 37600809^4=19359191^4+36895840^4+8526200^4

#12 楼

#13 楼

#14 楼

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {
            for (d = 1; d < 500000; d++) {
                if (prop(a, b, c, d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

for (a = 1; a < 100000; a++) {
    long int a4 = a * a * a * a;
    for (b = 1; b < 300000; b++) {
        long int a4b4 = a4 + b * b * b * b;
        for (c = 1; c < 500000; c++) {
            long int a4b4c4 = a4b4 + c * c * c * c;
            for (d = 1; d < 500000; d++) {
                if (a4b4c4 == d * d * d * d)
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

#15 楼

// Ref 
a = 95800; b = 217519; c = 414560; d = 422481


#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define a (95800 + 1)
#define b (217519 + 1)
#define c (414560 + 1)
#define d (422481 + 1)

// Adjust to some 128 bit type on your platform
// Can leave it at 64-bit if `num/den < 65535/d` or about 0.155 for initial testing
typedef uintmax_t uint128t;

#define pow2(x) ((uint128t)(x)*(x))
#define pow4(x) (pow2((uint_fast64_t)(x)*(x)))
void doit_times(unsigned num, unsigned den) {
  time_t t0;
  time(&t0);
  printf("Portion:%6g%%\t", 100.0 * num / den);
  uint_fast32_t a_max = (uint_fast32_t) round(1.0 * a * num / den);
  uint_fast32_t b_max = (uint_fast32_t) round(1.0 * b * num / den);
  uint_fast32_t c_max = (uint_fast32_t) round(1.0 * c * num / den);
  uint_fast32_t d_max = (uint_fast32_t) round(1.0 * d * num / den);
  printf("a_max:%6" PRIuFAST32 " b_max:%6" PRIuFAST32 //
      " c_max:%6" PRIuFAST32 " d_max:%6" PRIuFAST32 "\t", //
      a_max, b_max, c_max, d_max);
  uint_fast32_t q0 = 1;
  for (uint_fast32_t qa = q0; qa <= a_max; qa++) {
    assert(qa <= 0xFFFF || sizeof(uint128t) > 4);
    uint128t qa4 = pow4(qa);
    for (uint_fast32_t qb = qa; qb <= b_max; qb++) {
      uint128t sum_ab4 = qa4 + pow4(qb);
      uint_fast32_t qc = qb;
      uint128t sum_abc4 = sum_ab4 + pow4(qc);
      uint_fast32_t qd = qc + 1;
      uint128t qd4 = pow4(qd);

      // Loop until `c` to too large.
      for (;;) {
        if (sum_abc4 < qd4) {
          if (qc >= c_max) break;
          qc++;
          sum_abc4 = sum_ab4 + pow4(qc);
        } else if (sum_abc4 > qd4) {
          //if (qd >= d_max) break;
          qd++;
          qd4 = pow4(qd);
        } else {
          puts("Success!");
          printf("a %" PRIuFAST32 "\n", qa);
          printf("b %" PRIuFAST32 "\n", qb);
          printf("c %" PRIuFAST32 "\n", qc);
          printf("d %" PRIuFAST32 "\n", qd);
          fflush(stdout);
          exit(0);
        }
      }

    }
  }
  time_t t1;
  time(&t1);
  double dt = difftime(t1, t0);
  printf("t:%g\n", dt);
  fflush(stdout);
}

int main(void) {
  for (unsigned num = 1; num <= 1000; num++) {
    doit_times(num, 1000);
  }
  puts("Done");
  return -1;
}

Portion:   0.1% a_max:    96 b_max:   218 c_max:   415 d_max:   422 t:0
Portion:   0.2% a_max:   192 b_max:   435 c_max:   829 d_max:   845 t:0
Portion:   0.3% a_max:   287 b_max:   653 c_max:  1244 d_max:  1267 t:1
Portion:   0.4% a_max:   383 b_max:   870 c_max:  1658 d_max:  1690 t:1
Portion:   0.5% a_max:   479 b_max:  1088 c_max:  2073 d_max:  2112 t:2
Portion:   0.6% a_max:   575 b_max:  1305 c_max:  2487 d_max:  2535 t:3
Portion:   0.7% a_max:   671 b_max:  1523 c_max:  2902 d_max:  2957 t:5
Portion:   0.8% a_max:   766 b_max:  1740 c_max:  3316 d_max:  3380 t:6
Portion:   0.9% a_max:   862 b_max:  1958 c_max:  3731 d_max:  3802 t:10
Portion:     1% a_max:   958 b_max:  2175 c_max:  4146 d_max:  4225 t:12
Portion:   1.1% a_max:  1054 b_max:  2393 c_max:  4560 d_max:  4647 t:16
Portion:   1.2% a_max:  1150 b_max:  2610 c_max:  4975 d_max:  5070 t:21
Portion:   1.3% a_max:  1245 b_max:  2828 c_max:  5389 d_max:  5492 t:25
Portion:   1.4% a_max:  1341 b_max:  3045 c_max:  5804 d_max:  5915 t:32
Portion:   1.5% a_max:  1437 b_max:  3263 c_max:  6218 d_max:  6337 t:39
Portion:   1.6% a_max:  1533 b_max:  3480 c_max:  6633 d_max:  6760 t:48
Portion:   1.7% a_max:  1629 b_max:  3698 c_max:  7048 d_max:  7182 t:56
Portion:   1.8% a_max:  1724 b_max:  3915 c_max:  7462 d_max:  7605 t:67
Portion:   1.9% a_max:  1820 b_max:  4133 c_max:  7877 d_max:  8027 t:77
Portion:     2% a_max:  1916 b_max:  4350 c_max:  8291 d_max:  8450 t:91
Portion:   2.1% a_max:  2012 b_max:  4568 c_max:  8706 d_max:  8872 t:104
Portion:   2.2% a_max:  2108 b_max:  4785 c_max:  9120 d_max:  9295 t:120

#16 楼

#17 楼