/*
2000/08 三数教室
coded by Isaku WADA
 ┌─────────────────────────────┐
 │ ここに,それぞれ3個の数からなる組が3組ある。      │
 │(12,18,35)                       │
 │(10,27,28)                       │
 │(14,15,36)                       │
 │それぞれの組について,3個の数の和をとるとすべて65になる。 │
 │また積は7560となり,やはり同じになっている。ではこれを4組 │
 │に拡張し,「3数づつが4組あり,各組の和をとっても積をとって│
 │も同じになる」という数の組み合わせをひとつ見つけてほしい。│
 │ここに登場する数はすべて正の整数で,同じ数は使わない。余裕│
 │のある方は5組でも挑戦していただきたい。          │
 └─────────────────────────────┘
 ┌─────────────────────────────┐
 │・プログラムの説明                    │
 │                             │
 │まず、積が M=2520(=8*9*5*7) で、和の等しい組を探します。見│
 │つからなければ積が 2*M の場合を探します。同様に、積として │
 │Mの倍数(x)を次々と与えて、和の等しい組を探索していきます。│
 │                             │
 │より具体的には、積xが与えられたら、xを素因数分解します。次│
 │に、結果を利用して Value1 を x の約数の中から選びます。さ │
 │らに、Value2 を x/Value1 の約数の中から選びます。すると、 │
 │Value3(=x/Value1/Value2)が自動的に決まります。ここで   │
 │Value1 < Value2 < Value3でなければやり直します。うまくいっ│
 │たら、3数の和を求め、リストに登録します。        │
 │                             │
 │同様に、次々と(Value1,Value2,Value3)の組み合わせを探し、和│
 │を求めリストに登録してゆきます。同じ和を持つ組が組数(=N)以│
 │上になったら、リストを表示してプログラムを終了します。  │
 │                             │
 │もし、x の全ての組み合わせを調べても見つからなければ、次の│
 │xを与えなおし、繰り返していきます。M を 1 にすれば、必ず最│
 │小積が見つかりますが、極端に遅くなります。デフォルトのMを │
 │2520 にしたのは、N=4の最小積(=25200)とN=5の最小積(=83160) │
 │の最大公約数だからです。M を大きくすればするほど、高速にな│
 │りますが、解の積が大きくなるリスクが増します。      │
 └─────────────────────────────┘
 ┌─────────────────────────────┐
 │・感想                          │
 │                             │
 │解が見つかりさえすればよいので、最初は乱数を使って解いてい│
 │ました。しかし、膨大な時間がかかるのでやめて、組み合わせを│
 │調べるようにしました。それでも遅いので、素因数分解で2で割│
 │り切れる場合に、割り算ではなくシフト演算を使うようにしたら│
 │5倍も速くなりました。そして M の倍数だけを計算するように │
 │したら、さらに20倍速くなりました。今回の問題は、努力が速│
 │さに反映されるので、やりがいがありました。        │
 └─────────────────────────────┘
*/

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

/*┌─────────────────────────────┐
 │N : 組数                         │
 │                             │
 │LSI C-86 試食版なら最大 N=7                │
 │これはデータ領域の大きさが 64k に制限されているからです。 │
 │32 bit コンパイラなら最大 N=11                            │
 │これは、N>11 だと積が 32 bit 整数の上限を超えるからです。 │
 └─────────────────────────────┘*/
#define N 4

/*┌─────────────────────────────┐
 │M : 高速化のための乗数                  │
 │                             │
 │プログラムは、積が M の倍数の場合だけを計算します。       │
 │                                                          │
 │M を1にすれば、必ず最小積が見つかりますが、遅くなります。 │
 │デフォルトの M を 2520 にしたのは、N=4 の最小積(=25200)と │
 │N=5 の最小積(=83160)の最大公約数だからです。              │
 │M を大きくすればするほど、高速になりますが、              │
 │解の積が大きくなるリスクが増します。                      │
 │                                                          │
 │M=2520でN=2,3の場合はプログラムの結果よりも積が小さい解   │
 │が存在しますが、N>=4 の場合は最も積が小さい解が求まると   │
 │思われます。                                              │
 │N=10の場合2520なら約7分、25200なら約1分で求まります。   │
 │N=11 の場合 M=1437836400(Pentium!!! 550MHz)               │
 └─────────────────────────────┘*/
#define M 2520

/*┌─────────────────────────────┐
 │H : 各種配列の大きさ                    │
 │                                                          │
 │16 bit コンパイラなら  2048                               │
 │32 bit コンパイラなら 65536                               │
 │この定数はあまり変更する必要はありません                  │
 └─────────────────────────────┘*/
#define H ((int)(sizeof(int)==2?1<<11:1L<<16))

/*┌─────────────────────────────┐
 │P : 素因数分解したときの、素因数の種類の数の最大値    │
 └─────────────────────────────┘*/
#define P 10

/*┌────────────┐
 │HASH_FUNC : ハッシュ関数│
 └────────────┘*/
#define HASH_FUNCTION(X) ((int)(((((X)+12345)*(X))>>16)&(H-1)))

/*
 * clock()
 *
 * clock() をサポートしない MS-DOS コンパイラ用の自前の clock()
 * LSI C-86 試食版でも 1/1000 秒単位の経過時間を返すようになる
 * Visual C++、Borland C++ Builder、UNIX などでは自動的に無視される
 *
 * 最初に clock() を呼び出すと必ず0を返す。
 * 2回目以降は、最初の呼び出しからの経過時間を 1/1000 秒単位で返す
 * 経過時間が24日以上になると−1を返す
 * 閏年も考慮に入れている
 */
#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)
{
    static int day[12]={0,31,59,90,120,151,181,212,243,273,304,334};
    union REGS t,d; static long Date0=0,Time0=0; long Time1,Date1;

    t.h.ah=0x2c; d.h.ah=0x2a;
    intdos(&t,&t); intdos(&d,&d);
    Time1=t.h.dl*10L+t.h.dh*1000L+t.h.cl*60000L+t.h.ch*3600000L;
    Date1=d.x.cx*365L+d.x.cx/4+day[d.h.dh-1]+d.h.dl;
    if ((d.x.cx%4)==0&&d.h.dh<=2) Date1--;
    if (Date0==0)  { Date0=Date1; Time0=Time1; return 0; }
    if (Date1-Date0>=24) return-1;
    return(Date1-Date0)*86400000L+Time1-Time0;
}
#endif

/*┌────────┐
 │node_t : リスト │
 └────────┘*/
typedef struct {
    unsigned long sum; /* 和 */
    unsigned long value1,value2,value3; /* 3 数 */
    int SameNext; /* 同じ和を持つ次の組 */
    int DiffNext; /* 同じハッシュ値で異なる和を持つ次の組 */
    int count;    /* その和の出現回数 */
} node_t;

/*┌─────┐
 │変数の宣言│
 └─────┘*/
int NumOfPrimes;          /* 分解した時の素因数の種類数 */
unsigned long Primes[P];  /* 分解した時の素因数の値 */
int Max1[P];    /* 分解した時の素因数の個数 */
int Max2[P];    /* Max[P] の中間結果(Result2[P]の範囲)*/
int Result1[P],Result2[P]; /* 2(/3)数に対応する分解結果 */
unsigned long Value1,Value2,Value3;   /* 現在の 3 数 */
node_t*NodeTable;
int*HashTable;
int*UsedMarkers;    /* ハッシュ表の使用記録(クリア時に利用) */
int Index;          /* NodeTable[H] の使用中の領域の大きさ */
int IndexMax;       /* Index の最大値 */
double StartTime;   /* プログラム開始時の clock() の値 */

/*┌────────────────────────────┐
 │x を素因数分解する                                      │
 │素因数が Primes[] に求まり、対応する個数が Max に求まる │
 │素因数の種類数を返す                                    │
 └────────────────────────────┘*/
int Factorization(unsigned long*Primes,int*Max,unsigned long x)
{
    unsigned long divisor=3,prev=0; int NumOfPrimes;

    if ((x&1)==0) { /* x が偶数の場合 */
        NumOfPrimes=1;
        Max[0]=0;
        Primes[0]=2;
        do { x>>=1; Max[0]++; }
        while ((x&1)==0);
    } else NumOfPrimes=0;
    while (x!=1)
     if ((x%divisor)==0) {
         if (prev!=divisor) {
             Max[NumOfPrimes]=1;
             Primes[NumOfPrimes++]=divisor;
         } else Max[NumOfPrimes-1]++;
         x/=divisor;
         prev=divisor;
     } else divisor+=2;
    return NumOfPrimes;
}

/*┌─────────────────┐
 │結果を表示してプログラムを終了する│
 └─────────────────┘*/
void PrintResult(node_t*p)
{
    for (;p!=NodeTable;p=NodeTable+p->SameNext)
        printf("(%4lu,%4lu,%4lu)\n",p->value1,p->value2,p->value3);
    printf("使用した大きさ:%d/%d  ",IndexMax,H);
    printf("実行時間:%f 秒\n",(clock()-StartTime)/CLOCKS_PER_SEC);
    exit(0);
}

/*┌───────────────────────┐
 │次の組み合わせを求める。全部調べたら 0 を返す │
 └───────────────────────┘*/
int NextCombination(int*result,const int*Max)
{
    int i,j;
    
    for (i=0;i<NumOfPrimes;i++) if (result[i]!=Max[i]) break;
    if (i==NumOfPrimes) return 0;
    result[i]++;
    for (j=0;j<i;j++) result[j]=0;
    return 1;
}

/*┌─────────────────────────┐
 │Value1 < Value2 < Value3 ならば、記録し処理をする │
 └─────────────────────────┘*/
void SubProcedure(void)
{
    node_t*p,*q; int i,j,hash; unsigned long sum;

    for (Value2=1,i=0;i<NumOfPrimes;i++)
      for (j=Result2[i];j>0;j--) Value2*=Primes[i];
    if (Value2<=Value1) return;
    for (Value3=1,i=0;i<NumOfPrimes;i++)
      for (j=Max2[i]-Result2[i];j>0;j--) Value3*=Primes[i];
    if (Value3<=Value2) return;
    p=NodeTable+Index;
    sum=Value1+Value2+Value3;
    p->value1=Value1;
    p->value2=Value2;
    p->value3=Value3;
    hash=HASH_FUNCTION(sum);
    q=NodeTable+HashTable[hash];
    UsedMarkers[Index]=hash;
    while (q!=NodeTable)
     if (q->sum==sum) {
        /* 同じ和を持つ組が、既に登録済みの場合の処理 */
        p->SameNext=q->SameNext;
        q->SameNext=Index;
        if (++q->count>=N) PrintResult(q);
        break;
     } else q=NodeTable+q->DiffNext;
    if (q==NodeTable) {
        /* 同じ和を持つ組が無かった場合の処理 */
        p->SameNext=0;
        p->sum=sum;
        p->count=1;
        p->DiffNext=HashTable[hash];
        HashTable[hash]=Index;
    }
    if (++Index>=H) {
        printf("N を小さくするか H を大きくしてください\n");
        exit(1);
    }
    if (Index>IndexMax) IndexMax=Index;
}

/*┌───────┐
 │メインルーチン│
 └───────┘*/
void main(void)
{
    unsigned long x; int i,j,k=0;

    StartTime=clock();
    NodeTable=(node_t*)calloc(sizeof*NodeTable,H);
    HashTable=(int*)calloc(sizeof*HashTable,H);
    UsedMarkers=(int*)calloc(sizeof*UsedMarkers,H);
    printf("組数N:%d  配列の大きさH:%d  乗数M:%d\n",N,H,M);
    for (x=M;x-M<x;x+=M) { /* M の倍数を積として調べてゆく */
        if ((k++&255)==0)
        { fprintf(stderr,"%lu\r",x); fflush(stderr); }
        for (i=1;i<Index;i++) HashTable[UsedMarkers[i]]=0;
        Index=1; NumOfPrimes=Factorization(Primes,Max1,x);
        for (i=0;i<NumOfPrimes;i++) Result1[i]=0;
        do {
            for (Value1=1,i=0;i<NumOfPrimes;i++) {
                for (j=Result1[i];j>0;j--) Value1*=Primes[i];
                Max2[i]=Max1[i]-Result1[i];
                Result2[i]=0;
            }
            do SubProcedure();
            while (NextCombination(Result2,Max2));
        } while (NextCombination(Result1,Max1));
    }
    printf("積が unsigned long の扱える範囲を超えました\n");
}