编程快速求解丢番图方程

happy886rr
渐入佳境
渐入佳境
帖子: 45
注册时间: 2016年09月27日 16:11
拥有现金: 锁定
储蓄: 锁定
Has thanked: 14 times
Been thanked: 14 times
联系:

编程快速求解丢番图方程

帖子 #1 happy886rr » 2016年09月28日 17:50

形如方程1999a+1299b+1499c+2799d+2499e+3299f=2328450.在约束条件 a∈[0,607]、 b∈[0,534]、 c∈[0,60]、 d∈[0,660]、 e∈[0,100]、 f∈[0,209]下的整数解?


这是一个线性丢番图求解问题。在解线性丢番图方程式,必须学会使用降元法,降此六元丢番图将为两元丢番图方程式,然后利用两元丢番图的通解公式解决。
先化为ti型参数方程:

代码: 全选

{
    1999a+1299b=t2
    t2+1499c=t3
    t3+2799d=t4
    t4+2499e=t5
    t5+3299f=2328450    //其中t2~t5为参数
}

从最后一个ti型参数方程“t5+3299f=2328450”入手也就是t5=2328450-3299f,这是标准的二元线性丢番图,先找出其一个解{f=0,t5=2328450},
直接套二元丢番图公式(易求出)其通解为:

代码: 全选

{
    f=0+θ,               //由于题目还给出f<=209,所以θ∈[0,209]闭区
    t5=2328450-3299θ;    //其中θ∈[0,209]
}

这样f的解就都求出来,直接把求得的t5参值代入到倒数第二个ti型参数方程“t4+2499e=t5”中,继续重复调用二元线性丢番图丢番图的求解公式,得到全部e的解和t4参值。依次类推:
倒数第三个ti型参数方程
倒数第四个ti型参数方程
倒数第五个ti型参数方程
...


这种求解思路极易程序化,那就用C语言来实现一下

代码: 全选

对1999a+1299b+1499c+2799d+2499e+3299f=2328450等式同除以3,降低运算量至30%,
可化为(1999a+1499c+3299f)/3+(433b+933d+833e)=776150,
再化为 433b+933d=776150-833e-(1999a+1499c+3299f)/3

由于b∈[0,534], d∈[0,660]的值域较大,优先对b、d打表,即构建表m[sum=433b+933d]=[b,d]。考虑b、d均小于1000,那么我在打表时用一个int数就可以表示b、d这两个值,即

代码: 全选

m[sum]=b*1000+d;
b=m[sum]/1000,d=m[sum]%1000;  //在拆表时基底分离

当然在构建表数组m时最好使用int *m=(int*)calloc(int size, sizeof(int))自动初始化为0.剩下的几乎没难度,就是丢番图递推了,中间考虑到b的取值x最大才534,因此又省了个while循环,但这只是碰巧。
直接上代码:
Code: [全选] [展开/收缩] [Download] (Untitled.c)
  1. #include <stdio.h>
  2.  
  3. int main(int argc, char *argv[])
  4. {
  5.     int *m=(int*)calloc(776151, sizeof(int));
  6.     int a,b,c,d,e,f,sum,x,y,i;
  7.     for(b=0;b<=534;b++){
  8.         for(d=0;d<=660;d++){
  9.             if(sum=433*b+933*d,sum<28576){continue;}
  10.             else if(sum>776150){break;}
  11.             m[sum]=b*1000+d;
  12.         }
  13.     }
  14.     printf("求解完毕,请按任意键存为solutions.txt");
  15.     getchar();
  16.     FILE *fp=fopen("solutions.txt","w");
  17.     printf("正在保存,这需要一些时间和数GB硬盘空间!");
  18.     for(f=0;f<=209;f++){
  19.         for(c=0;c<=60;c++){
  20.             for(a=0;a<=607;a++){
  21.                 if((1999*a+1499*c+3299*f)%3!=0){continue;}
  22.                 for(e=0;e<=100;e++){
  23.                     if(sum=776150-833*e-(1999*a+1499*c+3299*f)/3,(sum<28576)||(m[sum]==0)){continue;}
  24.                     else if(sum>776150){break;}
  25.                     x=m[sum]/1000,y=m[sum]%1000;
  26.                     /*由于给出的b的取值x最大才534,因此连while循环也省了。
  27.                     while(x>0)
  28.                     {
  29.                         fprintf(fp,"%d,%d,%d,%d,%d,%d\n",a,x,c,y,e,f);
  30.                         x-=933,y+=433;i++;
  31.                     }
  32.                     */
  33.                     fprintf(fp,"%d,%d,%d,%d,%d,%d\n",a,x,c,y,e,f); 
  34.                 }
  35.             }
  36.         }
  37.     }
  38.     fclose(fp);
  39.     return 0;
  40. }

运行无需要时间,但打印需要很长时间,总计3GB求解结果。你可以把打印的fprintf注释掉,你就知道它的求解速度有多惊人。
上次由 happy886rr 在 2016年09月28日 18:40,总共编辑 1 次。

头像
523066680
Administrator
Administrator
帖子: 340
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 30 times
Been thanked: 27 times
联系:

Re: 编程快速求解丢番图方程

帖子 #2 523066680 » 2016年09月28日 18:23

哈哈知乎那个是你

happy886rr
渐入佳境
渐入佳境
帖子: 45
注册时间: 2016年09月27日 16:11
拥有现金: 锁定
储蓄: 锁定
Has thanked: 14 times
Been thanked: 14 times
联系:

Re: 编程快速求解丢番图方程

帖子 #3 happy886rr » 2016年09月28日 18:41

今天刚激活知乎号,否则还不能回复,以前只知道批处理之家。再说,去年我什么编程都不会,以前只是搞数学的,对编程不是很了解。

头像
523066680
Administrator
Administrator
帖子: 340
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 30 times
Been thanked: 27 times
联系:

Re: 编程快速求解丢番图方程

帖子 #4 523066680 » 2016年09月28日 20:16

happy886rr 写了:今天刚激活知乎号,否则还不能回复,以前只知道批处理之家。再说,去年我什么编程都不会,以前只是搞数学的,对编程不是很了解。


有数学这个基础,后面的路就要溜很多了。我现在就是受制于数学水平,工作后即使抽时间学也很难有读书那时候的效果了。

happy886rr
渐入佳境
渐入佳境
帖子: 45
注册时间: 2016年09月27日 16:11
拥有现金: 锁定
储蓄: 锁定
Has thanked: 14 times
Been thanked: 14 times
联系:

Re: 编程快速求解丢番图方程

帖子 #5 happy886rr » 2016年09月28日 20:35

我也好久没学了,都是初中的时候对数学比较热爱,就像你写贪吃蛇那会对批处理的热情。那会买了本砖头厚的奥数,天天研究。兴趣才是最好的teacher。


回到 “Math”

在线用户

用户浏览此论坛: 没有注册用户 和 1 访客