hcnak.blog

posted at 2016-04-10 23:48:03 +0000

C语言证明哥德巴赫猜想在10亿以内都成立（位操作）

这是老师布置的一道小组作业,前前后后一共布置了3次（根据C语言的学习进度，限制代码的使用）：证明10W内，1亿内，10亿内。

第一次请参见：C语言证明哥德巴赫猜想算法以及算法优化相关问题

第二次：暂未贴出

这里来说说10亿内的证明。

这里我们用的筛选法

要是还是用数组来储存质数的话，恐怕申请不到那么大的空间，所以这里就涉及到位操作。我用了unsigned long作为容器，有些同学用的是char作为容器。

用unsigend long作为容器，我们由于要对其每一位操作，所以我用一个数组来储存每次只有1个位为1的情况：（第一个全为0，是为了在后面的循环中表示方便）

unsigned long bits[33]={0x0,0x1
                ,0x2,0x4,0x8,0x10
                ,0x20,0x40,0x80,0x100
                ,0x200,0x400,0x800,0x1000
                ,0x2000,0x4000,0x8000,0x1000
                ,0x20000,0x40000,0x80000,0x10000
                ,0x200000,0x400000,0x800000,0x100000
                ,0x2000000,0x4000000,0x8000000,0x1000000
                ,0x20000000,0x40000000,0x80000000};

生成偶数表

（其实完全不用生成偶数表的，我们可以直接在证明循环里循环偶数来证明。但我们老师要求生成偶数表）

这里我们只需要把偶数数组（利用calloc()生成）里面的每一个元素与 0xAAAAAAAA 做或运算即可

因为0xAAAAAAAA的每一个偶数位都为1

接下来开始利用筛选法筛选出素数

首先我们要排除所有偶数，同上，我们把待筛选数组的每一个元素赋值0x55555555即可

0x55555555的每个奇数位都是1：

最后别忘了，我们还要挖掉1，把待筛选数组的第一个元素和bits[1]做异或运算，另外，2是最小的质数，所以我们还要把2加进去，再把待筛选数组的第一个元素和bits[2]做或运算即可。

接下来就是不断筛选其它素数。

我们只需要依次取出待筛选数组里面的每个数字，然后把他在1—-10亿以内的倍数都挖掉即可。

这里由于是位操作，就涉及到计算数字值的问题，因为我们是用的unsigned long所以我们要处理的数字应该是（待筛选数组为pary，i为该数组中第i个元素的下标，j为该数字在第i个元素所在的第j个位）：

long number = i * sizeof(unsigned long) * 8 + j;

注意sizeof得到的大小是字节，所以我们还要乘以8.

然后在main函数的证明中，我们同样需要定位该数字在pary数组中的位置，以确定其是否为质数。

(ii为pary中的第ii个元素位置，jj为在第ii个元素中第jj个位)

ii = number/(sizeof(unsigned long)*8);
jj = number%(sizeof(unsigned long)*8);

由于我们在筛选中，把是素数的都标记为1，非素数都标记为0，故通过以上简单的定位计算，即可确定（在筛选完成后）当前数字是否为质数。

完整代码：

#include <stdio.h>

#include < time.h >

#include < stdlib.h >

#
define MAX 15625000 * 2
//#define MAX 100
# define SIZE(sizeof(unsigned long) * 8)
//#define DEBUG    //输出A+B=C
//#define DEBUG_A  //输出素数表
/*哥德巴赫猜想之3 检查4 - 1,000,000,000所有偶数是否满足.
# 基于哥德巴赫猜想之2（completed@2015-11-22）修改
用double数组存放：1,000,000,000/(32_2) = 1562500 unsigned long(8) = 12500,000 byte = 12,207.04 kb = 12.21 MB 也就是说我们需要 1562500_2个unsigned long变量(之前以为我的电脑是8字节的long，结果是4字节的，故所需要2倍) 注意：unsigend long在我电脑上是32位的 */

unsigned long _pnum; //save the d-num
unsigned long_ parrary; //save the prime numbers

unsigned long bits[33] = {
  0x0,
  0x1,
  0x2,
  0x4,
  0x8,
  0x10,
  0x20,
  0x40,
  0x80,
  0x100,
  0x200,
  0x400,
  0x800,
  0x1000,
  0x2000,
  0x4000,
  0x8000,
  0x1000,
  0x20000,
  0x40000,
  0x80000,
  0x10000,
  0x200000,
  0x400000,
  0x800000,
  0x100000,
  0x2000000,
  0x4000000,
  0x8000000,
  0x1000000,
  0x20000000,
  0x40000000,
  0x80000000
};

const unsigned long dbits = 0xAAAAAAAA;
const unsigned long delodd = 0x55555555;
//奇数位为1
float InvSqrt(float x);
int isPrimeA(int num); //a very simple version of is prime
unsigned long numGen(int max); //生成偶数
unsigned long_ primeGenA(int max); //get all prime between 1--max

//struct nx{char x[64]}p;

int main() {
  clock_t s, f;
  unsigned int m, n, t = 2, T = 0;
  register int x = 6;
  int loops = MAX;
  //store how many loops do we need to demonstrate
  unsigned long tlp = 0, tt = 0, txx = 0, sxx = 0, pn = 0;
  s = clock();
  /////////////////////////////////////////////////////
  numGen(MAX);
  //generate d-numbers between min to max
  primeGenA(MAX);
  //generate prime
  for (x = 0; x < loops; x++) {
    for (t = 1; t <= SIZE; t++) {
      if (_(pnum + x) & bits[t])
      //是偶数
      {
        if (x == 0 && t < 3)
        //从4开始证明
        {
          continue;
        }
        tlp = x_ SIZE + t;
        //当前正在进行判断的偶数
        for (tt = 0; tt < tlp; ++tt, ++m) {
          txx = tt / SIZE;
          //第几个
          long sxx = tt % SIZE;
          //第几位5
          if ((parrary + txx) & bits[sxx])
          //找到第一个质数
          {
            pn = txx_ SIZE + sxx;
            //第一个质数的值
            n = tlp - pn;
            //得到下一个质数的值及其位置
            txx = n / SIZE;
            sxx = n % SIZE;
            if ( * (parrary + txx) & bits[sxx]) {

              #
              ifdef DEBUG
              printf("%d+%d=%d\t", pn, n, tlp);#
              endif
              break;
            }
          }
        }
      }
    }
  }
  //////////////////////////////////////////////////////
  ///////////////////////////这里的代码用来测试生成的素数表是否正确////////////////

  #
  ifdef DEBUG_A

  for (x = 0; x < MAX; x++) {
    for (t = 1; t <= SIZE; t++) {
      if (parrary[x] & amp; bits[t]) {
        printf("%ld\t", x * SIZE + t);
        //getchar();
      }
    }
  }#
  endif

  //////////////////////////////////////////////
  f
    = clock();
  free(pnum);
  free(parrary);
  printf("%.6lf seconds\n", (double)(f - s) / CLOCKS_PER_SEC);
  getchar();
  return 0;
}

float InvSqrt(float x) {
  float xhalf = 0.5 f_x;
  int i = _(int_) & x; // get bits for floating VALUE
  i = 0x5f375a86 - (i >> 1);
  // gives initial guess y0
  x = (float_) & i;
  // convert bits BACK to float
  x = x(1.5 f - xhalf);
  // Newton step, repeating increases accuracy return 1 / x; }

  unsigned long _numGen(int max)
  //标记偶数=1.COMP
  {
    int i = 0;
    int a = 0, b = 0, c = 1;
    pnum = (unsigned long_) calloc(max, sizeof(unsigned long));
    for (i = 0; i < max; i++) {
      pnum[i] = pnum[i] | dbits;
    }
    return pnum;
  }

  unsigned long _primeGenA(int max)
  //标记质数=1
  {
    int i, j;
    unsigned long ltp = 0, temp = 0, ii = 0, jj = 0, TOOOP = 0;
    //T = InvSqrt(MAX);
    parrary = (unsigned long_) calloc(max + 1, sizeof(unsigned long)); //把所有位设置为1 for(i=0;i<=max;i++) { _(parrary+i) = delodd;
    //首先就排除了偶数
  } //printf("bf=%ld\n",_parrary);
  _parrary ^= bits[1];
  //然后我们挖掉1_
  parrary |= bits[2];
  //然后我们加入2
  //printf("af=%ld\n",_parrary);
  //筛选法找出质数
  TOOOP = InvSqrt(max_SIZE);
  for (i = 0; i <= max; i++) {
    for (j = 1; j <= SIZE; j++) {
      //printf("_(parrary+i) & bits[j] = %ld , i=%d,j=%d,_(parrary+i)=%d,bits[j]=%d\n",_(parrary+i) & bits[j],i,j,_(parrary+i),bits[j]);
      if ((_(parrary + i) & bits[j])) {
        ltp = SIZE_ i + j;
        temp = 2 _ltp;
        //printf("--ltp=%ld,temp=%ld\n",ltp,temp);
        //TOOOP = InvSqrt(max_SIZE);
        while (temp <= TOOOP && ltp != 1) {
          ii = temp / SIZE;
          jj = temp % SIZE;
          if (_(parrary + ii) & bits[jj]) {
            (parrary + ii) ^= bits[jj];
          }

          // printf("------ii=%d,jj=%d,temp=%ld\t*(parrary+ii)=%d,bits[jj]=%d\n",ii,jj,temp,*(parrary+ii),bits[jj]);
          //getchar();

          temp += ltp;
        }
      }
    }
  }
  return parrary;
}

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