Challenge for Geeks

理系ネタ J

近頃流行のフィボナッチ数列、この1000項目を誤差なく表示せよ。
64 bit int でも100項目でオーバーフローするんで簡単じゃないよ。ちなみに２００桁程度。
64 bit で100項目で足りないから、 double でもそこからは誤差が出るよ。
まあ多倍長演算ライブラリを拾ってくれば終わりだけどね。

ありがとございます：

第二問

では多倍長とか使わず気楽にデカい計算できる問題：
10の10乗項目を10進数で表した時の下６桁を表示せよ。

ちなみに既に周期的になってるらしく、googol (10の100乗)項目も下６桁は同じ。googolplex項目はO(log n)でやっても宇宙が終わる。

第三問

そしてこれはかなり頭使わないとできないよ。
f(1)=3, f(2)=5, f(n+2)=f(n+1)*f(n) の場合に10の10乗項目を10進数で表した時の下６桁を表示せよ。
あ、これは思ってた以上に面倒だ。f(1)=3, f(2)=7の方が簡単なのでそちらを先に。

第四問

第三問に関連して。

  • 1. 数列 f(n) = a^n mod M, についてその周期性を求めよ。
  • 2. f(n+2) = (f(n)+f(n+1))mod M の周期をP(M)と置く。非常に大きなMについてP(M)を計算する効率的な方法はあるか？

ちなみに答えは知りません。2.は変な動きをします。

M 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
P(M) 3 8 6 20 24 16 12 24 60 10 24 28 48 40 24 36 24 18 60

（これに関する考察はid:ita:20071204:p2 に移動)

お、ちゃんと名前がついてたのかhttp://mathworld.wolfram.com/PisanoPeriod.html では天下一数値計算武闘会開催！できるだけ大きいP(M)/Mを求めよ！
追記：ありゃ、上から押えられるんだ。武闘会終了。http://d.hatena.ne.jp/maehara/20071205/1196782700
てことはO(M)個の独立な軌道があるのか。じゃあ別の問題：軌道の個数Q(M)についてQ(M)/Mの最大値を求めよ。

おまけ：上記1.と2.を高速に計算するルーチンがあるとして、それを使って第三問を高速に解くコードを示せ。

以下は1.を数値的に調べるコード

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

#define M 1000000

static int buf[M/2];

// table of (base^n)%M and its period, etc.
typedef struct PowTab
{
    int base;
// if gcd(M, base)!=1, find max n such that M%(base^n)=0
//   and store it to ofs
//  if gcd=1, ofs=0
    int ofs;
//  period = (gcd==1)? min n such that (base^n)%M=1 :
//                     min n such that (base^(ofs+n))%M = (base^ofs)%M
    int prd;
//  tab[i] = (base^i)%M
    int *tab;
} PowTab;

// gcd=1 case
void init_tab_plain(int base, PowTab *ptab)
{
    int n;
    int x=1;
    buf[0]=1;
    for(n=1;n<M;n++)
    {
	x = (x*base)%M;
	buf[n]=x;
	if (x==1)
	{
	    ptab->prd=n;
	    break;
	}
    }
    ptab->base = base;
    ptab->ofs=0;
    ptab->tab = (int *)malloc(ptab->prd*sizeof(int));
    for(n=0;n<ptab->prd;n++)
    {
	ptab->tab[n] = buf[n];
    }
    printf("Tab initialized: Base=%d p=%d ofs=%d\n",base, ptab->prd,ptab->ofs);
}

int gcd(int m, int n)
{
    if (n==0) return m;
    else return gcd(n, m % n);
}

// pre-calculate (base^n)%M table and its period
void init_tab(int base, PowTab *ptab)
{
    int g=gcd(M, base);
    printf("initializing: base=%d gcd=%d\n",base,g);

    if(g==1)
    {
	init_tab_plain(base, ptab);
	return;
    }

    int ofs=0;
    int m0=M;
    int x0=1;
    int n;

    while(m0 % g==0)
    {
	buf[ofs]=x0;
	x0 = (x0*base)%M;
	m0/=g;
	ofs++;
    }

    int x=1;
    for(n=ofs;n<M;n++)
    {
	buf[n]=x0*x;
	x = (x*base)%m0;
	if(x==1) break;
    }
    ptab->base = base;
    ptab->ofs = ofs;
    ptab->prd = n+1-ofs;
    ptab->tab = (int *)malloc((n+1)*sizeof(int));
    for(n=0;n<ptab->prd + ptab->ofs;n++)
    {
	ptab->tab[n] = buf[n];
    }
    printf("Tab initialized: p=%d ofs=%d\n",ptab->prd,ptab->ofs);
}

// return (base^exp)%M
int ModPow(PowTab *ptab, int exp)
{
    int p=ptab->prd;
    int o=ptab->ofs;
    exp -= o;
    if (exp<0) exp+=p;
    return ptab->tab[o+exp%p];
}