数学研发论坛

 找回密码
 欢迎注册
查看: 4971|回复: 12

[转载] 一个多变量丢番图方程

[复制链接]
发表于 2008-11-4 14:10:01 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
精华
http://topic.csdn.net/u/20081102 ... d-971e4ff1e0ea.html 中,northwolves问了一个问题,
i)对于四个不同自然数$a<=b<=c<=d$,满足$a^2+b^2+c^2+d^2=abcd$,求出一亿以内所有的(a,b,c,d)组合
ii)对于n个不同自然数$a_1<=a_2<=...<=a_n$,满足$a_1^2+a_2^2+...+a_n^2=a_1a_2...a_n$,求出一亿以内所有组合$(a_1,a_2,...,a_n)$

推荐
开始northwolves和medie2005试用了很多不错的方法求得了(特别是问题一的)部分解,然后medie2005发现一个很有趣的规律,这些解中很多数字会大量重复出现。在帖子中50楼,litaoye发现了一个重要的规律:
假设我们已经找到$a^2+b^2+c^2+d^2=abcd$的一组解(a,b,cd)
我们知道d是方程$X^2-abcX+a^2+b^2+c^2=0$的一个解,而这个方程还有另外一个解$e=abc-d$.
如果e<d,那么我们就可以将这个解降到另外一个更小一点的解(a,b,c,e)(这里它们顺序可能需要重派)。
所有反过来,如果对于某个解(a,b,c,d),我们分别假设其中a是某个解(b,c,d,x)对应的e;b是某个解(a,c,d,x)对应的e等等就可以通过一个基础解得出很多其他派生解。
所以对于这个问题,我们可以求出所有那些最基本的解,然后其他解都可以通过上面方法派生出来就可以了。
litaoye还分析出来对于4个变量的情况,所有其他解都可以通过基础解(2,2,2,2)直接或间接派生出来。

我们可以进一步分析更加一般的情况,所以如果有某个基础解$(a_1,a_2,...,a_n)$无法通过其他解派生出来(就是一个最基本的解),由于$a_n$是方程
$X^2-a_1a_2...a_{n-1}X+a_1^2+a_2^2+...+a_{n-1}^2=0$的一个解,这个方程有两个解
$X={a_1a_2...a_{n-1}+-sqrt(a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2))}/2$
由于这个解是最基本的,所以$a_n$必须是其中较小的解,由此我们得到
$a_n={a_1a_2...a_{n-1}-sqrt(a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2))}/2>=a_{n-1}$  (1)
于是我们可以得到
$(a_1a_2...a_{n-1}-2a_{n-1})^2>=a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2)$而且$a_1a_2...a_{n-2}>=2$

$2a_{n-1}^2+a_1^2+...+a_{n-2}^2>=a_1a_2...a_{n-2}a_{n-1}^2$   
所以
   $a_{n-1}<=sqrt({a_1^2+a_2^2+...+a_{n-2}^2}/{a_1a_2....a_{n-2}-2})$   (2)
由于$2a_{n-1}^2+a_1^2+...+a_{n-2}^2<=na_{n-1}^2$,所以我们还可以得到
   $2<=a_1a_2...a_{n-2}<=n$    (3)
而穷举满足(3)的自然数$a_1,a_2,...,a_{n-2}$非常容易,此后可以根据(2)式穷举对应的$a_{n-1}$,再根据(1)式判断是否存在对应的$a_n$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-4 14:17:24 | 显示全部楼层
再贴一个用gmp程序实现的对应代码,可以用来解出对于给定n的所有基本解或前面若干个解(比如前10000个解)。

  1. #include <stdio.h>
  2. #include <gmp.h>
  3. #include <math.h>
  4. #include <set>
  5. using namespace std;
  6. #ifndef N
  7. #define N 4
  8. #endif
  9. #define MAX_COUNT 10000
  10. typedef struct sDataEntry{
  11.     mpz_t datas[N];
  12.     mpz_t nc;
  13.     sDataEntry(){
  14.         int i;
  15.         for(i=0;i<N;i++)mpz_init(datas[i]);
  16.         mpz_init(nc);
  17.     }
  18.     sDataEntry(const sDataEntry& d){
  19.         int i;
  20.         for(i=0;i<N;i++)mpz_init(datas[i]);
  21.         mpz_init(nc);
  22.         for(i=0;i<N;i++){
  23.             mpz_set(datas[i],d.datas[i]);
  24.         }
  25.         mpz_set(nc,d.nc);
  26.     }
  27.     ~sDataEntry(){
  28.         int i;
  29.         for(i=0;i<N;i++)mpz_clear(datas[i]);
  30.         mpz_clear(nc);
  31.     }
  32.     void set_nc(){
  33.         int i;
  34.         mpz_t tmp;
  35.         mpz_init(tmp);
  36.         mpz_set_ui(nc,0);
  37.         for(i=0;i<N;i++){
  38.             mpz_mul(tmp,datas[i],datas[i]);
  39.             mpz_add(nc,nc,tmp);
  40.         }
  41.         mpz_clear(tmp);
  42.     }
  43.     bool operator<(const struct sDataEntry& s)const{
  44.         int i;
  45.         int c=mpz_cmp(nc,s.nc);
  46.         if(c!=0)return c<0;
  47.         for(i=N-1;i>=0;i--){
  48.             int c=mpz_cmp(datas[i],s.datas[i]);
  49.             if(c!=0)return c<0;
  50.         }
  51.         return false;
  52.     }
  53.     bool operator==(const struct sDataEntry& s)const{
  54.         int i;
  55.         for(i=N-1;i>=0;i--){
  56.             int c=mpz_cmp(datas[i],s.datas[i]);
  57.             if(c!=0)return false;
  58.         }
  59.         return true;
  60.     }
  61. }DataEntry;

  62. set<DataEntry> results;

  63. DataEntry t;
  64. int td[N];
  65. void search_last()
  66. {
  67.     long long mul=1LL;
  68.     long long sum=0LL;
  69.     int i;
  70.     if(td[N-3]==1)return;///No solution available
  71.     for(i=0;i<N-2;i++){
  72.         mul*=td[i];
  73.         sum+=(long long)td[i]*td[i];
  74.     }
  75.     if(mul<=2LL)return;///No solution available;
  76.     long long d=sum/(mul-2LL);
  77.     int up=(int)sqrt(d+0.01);
  78.     for(i=td[N-3];i<=up;i++){
  79.         td[N-2]=i;
  80.         long long a=sum+(long long)i*i;
  81.         long long b=mul*i;
  82.         long long c=b*b-4*a;
  83.         if(c<0)continue;
  84.         int u=(int)sqrt((double)c);
  85.         if(u*(long long)u<c)u=u+1;
  86.         if(u*(long long)u==c){///Found a solution now
  87.             if((b-u)&1)continue;///Invalid
  88.             td[N-1]=(int)((b-u)/2);
  89.             int j;
  90.             for(j=0;j<N;j++){
  91.                 mpz_set_ui(t.datas[j],td[j]);
  92.             }
  93.             t.set_nc();
  94.             results.insert(t);
  95.         }
  96.     }
  97. }

  98. void search(int last)
  99. {
  100.     int i;
  101.     if(last==N-2){
  102.         search_last();
  103.         return;
  104.     }
  105.     long long mul=1LL;
  106.     long long mul2;
  107.     double ub;
  108.     for(i=0;i<last;i++){
  109.         mul*=td[i];
  110.     }
  111.     if(mul>N)return;
  112.     ub=pow((double)N/(double)mul,1.0/(N-2-last))+0.01;
  113.     i=1;if(last>0)i=td[last-1];
  114.     for(;i<=(int)ub;i++){
  115.         td[last]=i;
  116.         search(last+1);
  117.     }
  118. }
  119. int tcount=0;
  120. void dump(const DataEntry& d)
  121. {
  122.     int i;
  123.     for(i=0;i<N;i++){
  124.         mpz_out_str(stdout,10,d.datas[i]);
  125.         printf("\t");
  126.     }
  127.     printf("\n");
  128. }

  129. void dumpall()
  130. {
  131.     set<DataEntry>::iterator it;
  132.     for(it=results.begin();it!=results.end();++it){
  133.         dump(*it);
  134.     }
  135. }

  136. bool gens(const DataEntry& d)///return false when tcount reaches MAX_COUNT
  137. {
  138.     mpz_t sum;
  139.     mpz_t tmp;
  140.     mpz_init(sum);mpz_init(tmp);
  141.     int i,j;
  142.     mpz_set_ui(sum,0);
  143.     for(i=0;i<N;i++){
  144.         mpz_mul(tmp,d.datas[i],d.datas[i]);
  145.         mpz_add(sum,sum,tmp);
  146.     }
  147.     for(i=0;i<N;i++){
  148.         if(i>0&&mpz_cmp(d.datas[i],d.datas[i-1])==0)
  149.             continue;
  150.         mpz_mul(tmp,d.datas[i],d.datas[i]);
  151.         mpz_sub(tmp,sum,tmp);
  152.         mpz_div(tmp,tmp,d.datas[i]);
  153.         if(mpz_cmp(tmp,d.datas[N-1])>0){
  154.             for(j=0;j<i;j++){
  155.                 mpz_set(t.datas[j],d.datas[j]);
  156.             }
  157.             for(j=i;j<N-1;j++){
  158.                 mpz_set(t.datas[j],d.datas[j+1]);
  159.             }
  160.             mpz_set(t.datas[j],tmp);
  161.             t.set_nc();
  162.             set<DataEntry>::iterator lit;
  163.             lit=results.find(t);
  164.             if(lit==results.end()){
  165.                 results.insert(t);
  166.                 tcount++;
  167.             }
  168.         }
  169.     }
  170.     mpz_clear(sum);mpz_clear(tmp);
  171.     return tcount<MAX_COUNT;
  172. }

  173. void gen_more()
  174. {
  175.     set<DataEntry>::iterator it;
  176.     for(it=results.begin();it!=results.end();++it){
  177.         if(!gens(*it))
  178.             break;
  179.     }
  180.     mpz_t maxt;
  181.     mpz_init(maxt);
  182.     mpz_set(maxt,it->nc);
  183.     set<DataEntry>::iterator nit;
  184.     for(nit=results.begin();nit!=results.end();++nit){
  185.         if(mpz_cmp(nit->datas[N-1],maxt)<=0)
  186.             dump(*nit);
  187.     }
  188.     mpz_clear(maxt);
  189. }

  190. int main()
  191. {
  192.     int i;
  193.     search(0);
  194. #ifdef IONLY
  195.     dumpall();
  196. #else
  197.     gen_more();
  198. #endif
  199. }
复制代码
我们知道n=2时方程没有非平凡解(也就是只有全0的解)
而利用这个程序还可以知道在100以内,在n分别为
2,6,9,11,12,15,16,18,20,21,24,29,32,33,36,41,42,45,48,50,51,56,57,60,66,72,76,77,81,82,84,90,96,99时没有非平凡解。而对于100以内(包括100)的其他n都有非平凡解。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-5 08:04:30 | 显示全部楼层
附件中文件给出了10000以内所有对应方程没有非平凡解的n,总共1360个,毫无规律,是否同素数的分布规律有得一拼,而且可能频率也差不多
sqm.tar.gz (3.14 KB, 下载次数: 5)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-5 08:34:05 | 显示全部楼层
比如n = 2
应该证明确实没非平凡解
但是n = 6是否在整个正整数区域均无非平凡解呢?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-5 08:45:16 | 显示全部楼层
原帖由 无心人 于 2008-11-5 08:34 发表
比如n = 2
应该证明确实没非平凡解
但是n = 6是否在整个正整数区域均无非平凡解呢?

是的,已经用计算机证明了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-5 10:10:55 | 显示全部楼层
我的意思是严格证明
毕竟某些方程的最小整数解是
很大的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-5 11:44:51 | 显示全部楼层
是严格证明。前面已经通过数学方法将最小解规约到一个很小的范围。然后用计算机穷举那个小范围就可以了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-5 14:08:59 | 显示全部楼层


知道了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2010-1-30 13:29:40 | 显示全部楼层
又看到一棵3叉树,也许会于我有用,收下了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2021-10-3 22:42:38 来自手机 | 显示全部楼层
oeis.org/A146974的第一个链接被标记为无法打开,楼主能修复它吗
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2021-10-24 21:28 , Processed in 0.076429 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表