/*
2002/08 整数解の怪
Hideki ANDOさんの解答にIsaku WADAが手を加えたもの
プログラムの説明は ANDO さんの手によるもの

 ┌────────────────────────────┐
 │ 以下の2つの式を同時に満たす正の整数解を求めたい。  │
 │   x^2+y^2=u^2                  │
 │   x^2−xy+y^2=v^2               │
 │ これを満たす解を1つだけ示しておく。         │
 │   x=8,y=15,u=17,v=13       │
 │ ただし,ここでは「x<y」とする。ほかにも解はいくつも│
 │あるので,この例とは別の解を1つだけお答えいただきたい。│
 │もちろん,この例の数をすべて同じ整数倍したものは,別の解│
 │とはしない。                      │
 └────────────────────────────┘
 ┌────────────────────────┐
 │解答:x=8109409,y=10130640,u=12976609,v=9286489 │
 └────────────────────────┘
 ┌──────────────┬─────┬───────┐
 │コンパイラ          │実行時間 │ファイルサイズ│
 ├──────────────┼─────┼───────┤
 │Microsoft Visual C++ 6.0    │ 0.2312 秒│  36864 byte  │
 │Microsoft Visual C++ .NET   │ 0.2359 秒│  34816 byte  │
 │Borland C++ 5.5 for Win32   │ 0.2281 秒│  53760 byte  │
 │gcc-2.953(djgpp)            │ 0.2198 秒│ 108465 byte  │
 │gcc-2.953(cygwin)           │ 0.2235 秒│  19176 byte  │
 │LSI C-86 Ver 3.30 試食版  │ 1.1920 秒│ 13946 byte │
 │Turbo C Ver 1.5(small model)│ 3.8660 秒│ 11402 byte │
 └──────────────┴─────┴───────┘
    CPU:Pentium4 2.52GHz       OS:Windows XP
┌──────────────────────────────┐
│・プログラムの説明                     │
│                              │
│(x, y, u) のような数の組ををピタゴラス数と呼び、(2kr, k^2 - │
│r^2, k^2 + r^2) の前2つをx<yの大小関係に応じて入れ換えれ│
│ば、全てこれの整数倍で表現できます。ただし、k,rはk>rで、│
│片方が偶数、もう片方が奇数である自然数です。また、k,rが互い│
│に素であるという条件を入れても成り立ちますが、かえって時間が│
│かかるためチェックしていません。              │
│                              │
│vの式をk,rで表すと、( k(k-r), r(k+r), v) が実はピタゴラス│
│数になることがわかります。これもkが偶数なら(2mn, m^2 - n^2,│
│m^2 + n^2)、奇数なら前2つを交換したもので表せます。要するに│
│k,rとm,nの対応関係が成り立つような(k, r, m, n)をみつけれ│
│ば、(x, y, u, v)が求まったことになります。         │
│                              │
│rが偶数のときr↑⇒n↓と仮定すると(↑は増加、↓は減少を意味│
│します)、r↑⇒k(k-r)↓⇒m^2-n^2↓⇒m↓(∵n↓の仮定)⇒2mn↓│
│⇒r(k+r)↓⇒r↓となり矛盾するので、r↑⇒n↑となります。また、r│
│が奇数のときはr↓⇒n↓と仮定すると、r↓⇒r(k+r)↓⇒m^2-n^2↓│
│⇒m↓(∵n↓の仮定)⇒2mn↓⇒k(k-r)↓⇒r↑となり矛盾するので、│
│r↓⇒n↑となります。                     │
│                              │
│kが奇数のとき、rを偶数の範囲で増やしていきながらnを小さい│
│方から増やしていき、m,nが整数になるかどうかをチェックして│
│いきます。ここで、r↑⇒n↑という関係を用いれば、各々のrに│
│おけるnの初期値を直前のrにおけるnの終了値にできます。kが│
│偶数のときは、rを奇数の範囲で減らしながらnを増やしてチェッ│
│クします。こちらもr↓⇒n↑という関係を用いれば同様にnの値│
│を次のrにおいて引き継ぐことができます。          │
│                              │
│なお、安藤解では使われていませんが、k,rの内の偶数の方は4の│
│倍数でなくてはなりません。m,nの両方が奇数であるとすると、│
│2mn と m^2 - n^2 の両方が偶数になり、そうするとk,rの両方が│
│偶数になってしまいます。つまり、m,nの内の少なくとも一方が偶│
│数でなくてはならず、そうすると 2mn は4の倍数になり、k,rの│
│内の偶数の方が4の倍数になります。この事実を用いると実行時間│
│を相当に削減できます。                   │
│                              │
│問題文中の解答例の排除は、元の安藤解ではk≠4rで実現してい│
│ましたが、この和田解においては既出のk/r比を複数記憶しておく│
│ことができるようになっています。              │
└──────────────────────────────┘
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/*┌────────┐
 │定数・変数の宣言│
 └────────┘*/
#define K 3247
double Log[100]; int LogSize,k,r,n;

/*┌───────┐
 │解の確認と表示│
 └───────┘*/
void search(long s,long t)
{   
    long a,m,x,y,u,v; double ratio; int i;
    
    for(;;n++) { m=s/n; a=(m+n)*(m-n); if(a<=t) break; }
    if(a!=t||s%n) return;
    if (k&1) m=(long)r*(k+r); else m=(long)k*(k-r);
    m>>=1; m/=n; x=(long)k*r<<1; y=(long)k*k-(long)r*r;
    if (x>y) { a=x; x=y; y=a; }
    ratio=(double)y/x;
    for (i=0;i<LogSize;i++) if (ratio==Log[i]) break;
    if (i<LogSize) return; else Log[LogSize++]=ratio;
    u=(long)k*k+(long)r*r; v=m*m+(long)n*n;
    printf(" k=%d x=%ld y=%ld u=%ld v=%ld\n",k,x,y,u,v);
}

/*┌────┐
 │解の探索│
 └────┘*/
void Solve(void)
{   
    for(k=2;k<=K;k++)
     if (k&1) for(n=1,r=2;r<k;r+=2)
                  search((long)r*(k+r)>>1,(long)k*(k-r));
     else     for(n=1,r=k-1;r>0;r-=2)
                  search((long)k*(k-r)>>1,(long)r*(k+r));
}

/*┌────────────────────┐
 │MS-DOS コンパイラ用のシンプルな clock() │
 └────────────────────┘*/
#if CLOCKS_PER_SEC == 1
#undef CLOCKS_PER_SEC
#endif
#ifndef CLOCKS_PER_SEC
#define CLOCKS_PER_SEC 1000
#include <dos.h>
static long clock(void) {
    union REGS t; t.h.ah=0x2c; intdos(&t,&t);
    return t.h.dl*10L+t.h.dh*1000L+t.h.cl*60000L+t.h.ch*3600000L;
}
#endif

/*┌───────┐
 │メインルーチン│
 └───────┘*/
int main(int argc,char**argv)
{
    double StartTime=clock(); long i,num;

    if (argc==2) num=atol(argv[1]); else num=1;
    for (i=0;i<num;i++) Solve();
    printf("実行時間 %.3f 秒\n",(clock()-StartTime)/CLOCKS_PER_SEC);
    return 0;
}