hcnak.blog

posted at 2015-11-03 22:32:39 +0000

C语言证明哥德巴赫猜想算法以及算法优化相关问题

上周，我们C语言老师给我们布置了一道小组作业：用C语言编写程序，证明哥德巴赫猜想在4—1000,000之间成立。最后成绩以各组算法在他机器上运行时间为准。特别注意：只能使用《C语言程序设计》前五章的知识，即不能用指针，不能用数组，不能用自定义函数等等…

我们才学C语言没多久，就布置这样的题目？跪了。不过大体想了想，要证明哥德巴赫猜在4—1000,000之间成立，其实还是简单的，关键是如何优化算法，让它的时间最短。

首先我们先来熟悉下哥德巴赫猜想：哥德巴赫猜想主张每个大于等于4的偶数都是哥德巴赫数——可表示成两个素数之和的数。

基本思路

来说说基本解题思路：对于任意给定偶数a，我们要先把它拆分成x和y，并满足x+y=a。然后循环判断x与y都为质数即可。

所以就有了如下代码：

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

int main() {
  clock_t s, f;
  int x = 4, m, n, t = 2, r = 0, T = 0, R = 0;
  s = clock();
  //process start 
  //第一个循环，用来取出4——1000000之间的，所有偶数x 
  for (; x < 1000000; x += 2) {
    m = x / 2;
    n = m;
    R = m;
    //这个循环，用来对取出的偶数x的分解出的2个数做x/2次判断因为m和n变到0和x都需x/2次自减自加 
    for (r = 0; r < R; r++, m--, n++) {
      T = sqrt(n);
      //这个循环用来判断偶数x分解出的2个数m，n是否为质数，因为 n>=m 恒成立，所以我们做sqrt(n)次判断 
      for (t = 2; t < T; t++) {
        //如果n和m中有一个数能够被大于2的整数i整除，则跳出这个for循环，对m，n变换后继续判断 
        if (m % t == 0 || n % t == 0) {
          break;
        }
      }
    }
    //如果循环结束的时候，n！>=x或者m!<=0则说明，这个偶数不满祝哥德巴赫的猜想 
    if ((!(n >= x)) || (!(m <= 0))) {
      printf("%d不满足n", x);
    }
  }
  //process compltet 
  f = clock();
  printf("%.3lf secondsn", (double)(f - s) / CLOCKS_PER_SEC);
  return 0;
}

但是，这个算法在我的电脑上大概需要跑7–9秒才可以证明完毕。所以就引出了如何优化的问题。

优化

我们看到，上面代码对于分解出来的每个因子x和y，判断质数都是从2开始，一直到√x或者√y，所以，这里就有了第一个优化方法：

只判断奇数，即只用奇数去除待判断数字m或n

//第一个循环，用来取出4——1000000之间的，所有偶数x 
for (; x < 1000000; x += 2) {
  m = x / 2;
  n = m;
  R = m;
  //这个循环，用来对取出的偶数x的分解出的2个数做x/2次判断因为m和n变到0和x都需x/2次自减自加 
  for (r = 0; r < R; r++, m--, n++) {
    T = sqrt(n);
    //这个循环用来判断偶数x分解出的2个数m，n是否为质数，因为 n>=m 恒成立，所以我们做sqrt(n)次判断 
    for (t = 3; t < T; t += 2) {
      //如果n和m中有一个数能够被大于2的整数i整除，则跳出这个for循环，对m，n变换后继续判断 
      if (m % t == 0 || n % t == 0) {
        break;
      }
    }
  } //如果循环结束的时候，n！>=x或者m!<=0则说明，这个偶数不满祝哥德巴赫的猜想 
  if ((!(n >= x)) || (!(m <= 0))) {
    printf("%d不满足n", x);
  }
}

其次呢，我们可以把判断

x和y是否为质数中，大的那个数字，放到最内层循环，来减少CPU跨切循环层次数

 T = sqrt(m);
 //首先判断m是否为质数，如果是质数，就进入到下个if块，判断n是否为质数 
 //如果m不是质数，跳出#2，开始新一次#1，并对m--，n++，因为一个不满足，这一对就不满足 
 for (t = 3; t <= T; t += 2) //#2 
 {
   if (m % t == 0) {
     break;
   }
 } //得到了第一个质数m，我们再判断n是否为质数，是的话就跳出#1，不是的话，回到#2继续寻找下一个质数m 
 if (t > T) {
   n = x - m;
   T = sqrt(n); //判断n是否为质数。如果是，就执行下个if块，并跳出#1，回到#0，判断下一个偶数是否满足
   for (t = 3, T += 1; t <= T; t += 2) //#3
   {
     if (n % t == 0) {
       break;
     }
   }
   if (t > T) {
     TT = 0;
     //printf(" %d+%d=%d",n,m,x); 
     continue;
   }
 }
 ++m, ++m;
 goto label;

再者，我们每次循环的过程中都要判断x%i==0，即x能否被i整除，在以上代码中，我们可以看到：每次判断x能被i整除后，我们用的是break语句，跳出这个循环，继续执行循环后的语句（即if块）。但是由于

我们已经知道这个数不是质数。所以我们可以用一个goto语句，直接跳到顶部，改变x和y的值，继续判断。这样便去除了冗余的if判断

做完以上优化后，最终版本如下：

//第一组作品，by Kanch 算法优化参考了：http://m.blog.csdn.net/blog/sunjiajiang/7887724
#include < stdio.h > #include < time.h > #include < math.h >

  int main() {
    clock_t s, f;
    int x = 6, m, n, t = 2, T = 0, TT = 0;
    s = clock();
    //第一个循环，用来取出4——1000000之间的，所有偶数x 
    for (; x < 1000000; ++x, ++x) //#0 { 
      //这个循环，用来对取出的偶数x的分解出的2个数做x/2次判断因为m和n变到0和x都需x/2次自减自加 
      //由于每次判断后，都要对m++,n--，所以，m+n=x 
      m = 3; 
      label_dm:
      T = sqrt(m);
    //首先判断m是否为质数，如果是质数，就进入到下个if块，判断n是否为质数
    //如果m不是质数，跳出#2，开始新一次#1，并对m--，n++，因为一个不满足，这一对就不满足
    for (t = 3; t <= T; t += 2) //#2
    {
      if (m % t == 0) {
        //break;
        ++m, ++m;
        goto label_dm;
      }
    }
    //得到了第一个质数m，我们再判断n是否为质数，是的话就跳出#1，不是的话，回到#2继续寻找下一个质数m
    if (t > T) {
      n = x - m;
      T = sqrt(n);
      //判断n是否为质数。如果是，就执行下个if块，并跳出#1，回到#0，判断下一个偶数是否满足
      for (t = 3, ++T; t <= T; t += 2) //#3
      {
        if (n % t == 0) {
          break;
        }
      }
      if (t > T) {
        TT = 0;
        //printf(" %d+%d=%d",n,m,x);
        continue;
      }
    }
    ++m, ++m;
    goto label_dm;

    //因为每次最外层循环开始前我都会把TT设置为1，这样，只要无法执行到51行，就说明这个偶数不满足哥德巴赫猜想啊
    if (TT) {
      printf("%d不满足n", x);
      TT = 1;
    }
  }
f = clock();
printf("%.6lf secondsn", (double)(f - s) / CLOCKS_PER_SEC);
return 0;

}

<2015-11-5日更新，这里再补上一张算法解释图>

kanch@20151103

last updated on 2022-07-27 01:57:54 +0000