bzoj 1002: [FJOI2007]轮状病毒(生成树计数,高精度)

题目链接

1002: [FJOI2007]轮状病毒

Time Limit: 1 Sec   Memory Limit: 162 MB
Submit: 4816   Solved: 2628
[ Submit][ Status][ Discuss]

Description

　　轮状病毒有很多变种，所有轮状病毒的变种都是从一个轮状基产生的。一个N轮状基由圆环上N个不同的基原子
和圆心处一个核原子构成的，2个原子之间的边表示这2个原子之间的信息通道。如下图所示

　　N轮状病毒的产生规律是在一个N轮状基中删去若干条边，使得各原子之间有唯一的信息通道，例如共有16个不
同的3轮状病毒，如下图所示

　　现给定n(N<=100)，编程计算有多少个不同的n轮状病毒

Input

　　第一行有1个正整数n

Output

　　计算出的不同的n轮状病毒数输出 

题解：

一下是转自VFleaking的博客：(粘贴过来后格式变了，推荐点下面的链接！！！)

博客链接

开篇的矩阵是基尔霍夫矩阵。

用基尔霍夫矩阵（教程：http://www.4ucode.com/Study/Topic/1940063），使用高斯消元解行列式，时间复杂度O(n^3)似乎可以AC。

其实如果你仔细观察矩阵，可以发现它是这样的：（消去了0号点，就是病毒中间的那个圆坨坨）

 3  -1  0  0  0  ...  0  0  0  -1

 -1  3  -1  0  0  ...  0  0  0  0

 0  -1  3  -1  0  ...  0  0  0  0

 0  0  -1  3  -1  ...  0  0  0  0

 0  0  0  -1  3  ...  0  0  0  0

 ...  ...  ...  ...  ...  ...  ...  ...  ...  ...

 0  0  0  0  0  ...  3  -1  0  0

 0  0  0  0  0  ...  -1  3  -1  0

 0  0  0  0  0  ...  0  -1  3  -1

 -1  0  0  0  0  ...  0 0   -1  3

下面就简单了。设F(n)为n轮状病毒的个数，然后自己手算一下上面的矩阵，

可以发现，F(n) = 3F(n - 1) - F(n - 2) + 2。

由此可以得到O（n）的算法。

我无聊，用矩阵乘法+快速幂进行优化……O（logn）

=============================================

update：

吐槽：你说是O(logn)就是O(logn)？敢不敢把高精度计算的时间也算上啊！慢得要死啊！自己弱就不要写题解啊！我当年是sb吗！！！

同学们，只要写高精度加减暴力递推就行了～

我不知道我上面这种含糊不清的题解导致了多少人晕晕地混过去了，深表歉意。时隔一年半，我把详细证明补上 。

/* 

这一段神犇无视啊……

来看这篇题解的同学八成是跟我当年一样是才开始做衡阳八中的新手。

这道题新手可以拿来练练高精度。

估计也没人喜欢看证明的吧。但是对于新手还是建议去看看基尔霍夫矩阵。

这一篇论文挺不错的： http://wenku.baidu.com/view/0c086741be1e650e52ea990e.html

*/

有公式恐惧症的同学，治病的时候到了！！

首先行列式有很多性质，第a行 * k 加到 第b行上去，行列式的值不变。

三角行列式的值等于对角线元素之积。

第a行与第b行互换，行列式的值取反。

常数 * 行列式，可以把常数乘到某一行里去。

如果你行列式不是很熟～建议先google行列式～不然下面会看晕～

那么我们现在对行列式进行变换，我们把第一行与第二行交换，再把第二行与第三行交换...，再把第n - 1行与第n行变换。这些都是行列式初等行列变换，得到新的行列式：

 -1 3 -1 0 0 ... 0 0 0 0

 0 -1 3 -1 0 ... 0 0 0 0

 0 0 -1 3 -1 ... 0 0 0 0

 0 0 0 -1 3 ... 0 0 0 0

 ... ... ... ... ... ... ... ... ... ...

 0 0 0 0 0 ... 3 -1 0 0

 0 0 0 0 0 ... -1 3 -1 0

 0 0 0 0 0 ... 0 -1 3 -1

 -1 0 0 0 0 ... 00  -1 3

 3 -1 0 0 0 ... 0 0 0 -1

这个行列式跟一开始的那个行列式的值不一定相等。

因为我们是通过n - 1次交换行的操作得到的，

为了说话方便我们称一开是的那个行列式为A，上面我刚写的是B

那么由行列式性质得：A = (-1)^(n - 1) * B

-1好烦人对吧……

好，那么现在可以正大光明地处理B了～

利用行列式性质，来手算这个行列式。之所以刚才有那么一步，就是为了方便手算。

因为观察B，发现就只剩下左下角的-1、3、-1三个倒霉了。

首先看倒数第二行：

 -1 0 0 0 0 ... 00  -1 3

用第一行的：

 -1   3   -1   0   0   ...   0   0   0   0

乘以-1来消，得：

 0-3 1 0 0 ... 00  -1 3

然后再用第二行的：

 0 -1   3   -1   0   ...   0   0   0   0

乘以-3来消，得：

 00 -8 3 0 ... 00  -1 3

这样就有了初步感觉了～

现在把这个过程一般化，假设现在正在处理的这一行的第k个和第k+1个时是这样的：

0 0 ... F(k) G(k) 0 0 ... -1 3

总能找到上面的某一行：

0 0 ... -1 3 -1 0 ... 0 0

乘以F(k)来消，根据我们设的数，应该得到的是：

0 0 ... 0 F(k + 1) G(k + 1) 0 ... -1 3

于是得到：

F(k + 1) = G(k) + 3 * F(k)

G(k + 1) = -F(k)

整合一下：

F(k + 1) = 3 * F(k) - F(k - 1)

从初始的行和消了一次之后的行中取得边界条件：F(1) = -1，F(2) = -3

经过了漫长的消来消去，我们最终可以把倒数第二行变为：

 0 0  0   0   0   ...  F(n - 3) G (n - 3)     -1   3

然后用倒数第四行：

 0  0  0  0  0  ...  -1  3  -1  0

乘以F(n - 3)来消：

 0 0  0   0   0   ...  0 F (n - 2)    G(n - 2)  - 1   3

还剩下一步……用倒数第三行：

 0  0  0  0  0  ...  0  -1  3  -1

乘以F(n - 2)来消：

 0 0  0   0   0   ...  0   0  F(n - 1)  - 1  G(n - 1) +  3

好现在搞定了倒数第二行，来看看成果：（设f = F(n - 1) - 1, g = G(n - 1) + 3）

 -1 3 -1 0 0 ... 0 0 0 0

 0 -1 3 -1 0 ... 0 0 0 0

 0 0 -1 3 -1 ... 0 0 0 0

 0 0 0 -1 3 ... 0 0 0 0

 ... ... ... ... ... ... ... ... ... ...

 0 0 0 0 0 ... 3 -1 0 0

 0 0 0 0 0 ... -1 3 -1 0

 0 0 0 0 0 ... 0 -1 3 -1

 0 0 0 0 0 ... 00  f g

 3 -1 0 0 0 ... 0 0 0 -1

好，现在来搞倒数第一行。

和倒数第二行的方法是类似的。

再设函数H(k)和I(k)，意义与F(k)、G(k)类似,得：

H(k + 1) = I(k) + 3 * H(k)

I(k + 1) = -H(k)

其实跟F、G的递推式是一样的我会乱说？

H(k + 1) = 3 * H(k) - H(k - 1)

边界条件是： H(1) = 3，H(2) = 8

最后使劲搞一搞，倒数第一行就成了：

 0  0   0   0   0   ...   0   0  H(n - 1)  I(n - 1)  - 1

再来看成果： （设h = H(n - 1), i = I(n - 1) - 1）

 -1 3 -1 0 0 ... 0 0 0 0

 0 -1 3 -1 0 ... 0 0 0 0

 0 0 -1 3 -1 ... 0 0 0 0

 0 0 0 -1 3 ... 0 0 0 0

 ... ... ... ... ... ... ... ... ... ...

 0 0 0 0 0 ... 3 -1 0 0

 0 0 0 0 0 ... -1 3 -1 0

 0 0 0 0 0 ... 0 -1 3 -1

 0 0 0 0 0 ... 00  f g

 0 0 0 0 0 ... 0 0 h i

f是不会为0的吧……

那么用倒数第二行来消倒数第一行，得：

 0  0   0   0   0   ...   0   0  0  i - g * (h / f)

现在这个行列式已经是三角行列式了。

那么它的值就是对角线元素之积。

于是乎：B = (-1) * (-1) * (-1) * ... * (-1) * f * (i - g * (h / f))

一共有n - 2个-1

那么如前文所述：A = (-1)^(n - 1) * B

又因为：B = (-1)^(n - 2) * (f * i - g * h)

于是有：A = (-1)^(2n - 3) * (f * i - g * h) = -f * i + g * h

来来……f、g、h、i的值带入：

A = -(F(n - 1) - 1) * (I(n - 1) - 1) + (G(n - 1) + 3) * H(n - 1)

把H、I的值带入：

A = -(F(n - 1) - 1) * (-H(n - 2) - 1) + (-F(n - 2) + 3) * H(n - 1)

然后再展开咯。。回忆下F、H的递推式～

A = F(n - 1) * H(n - 2) + F(n - 1) - H(n - 2) - 1 - F(n - 2) * H(n - 1) + 3 * H(n - 1)

=H(n) + F(n - 1) + F(n - 1) * H(n - 2) - F(n - 2) * H(n - 1) - 1

= H(n) + F(n - 1) + |F(n - 1) H(n - 1)| - 1

|F(n - 2) H(n - 2)|

那啥那个长得丑丑的占了两行的表示的是二阶行列式。

发现不能化简了？

没关系！

在行列式上动动手脚吧！

【FH定理】

对于任意大于2的k有：

|F(k - 1) H(k - 1)| = |F(2) H(2)|

|F(k - 2) H(k - 2)| |F(1) H(1)|

【FH定理的证明】

对于行列式：

|F(k - 1) H(k - 1)|

|F(k - 2) H(k - 2)|

把行列式最下面的行取反，则行列式的值取反：

-|F(k - 1) H(k - 1) |

|-F(k - 2) -H(k - 2)|

把行列式的上面的行乘以3加到下面去：

-|F(k - 1) H(k - 1) |

|F(k - 1) * 3 - F(k - 2) H(k - 1) * 3 - H(k - 2)|

特意构造的递推式出现了：

-|F(k - 1) H(k - 1)|

|F(k) H(k) |

有点眉目了～把第一行与第二行调换位置，行列式的值取反：

|F(k) H(k) |

|F(k - 1) H(k - 1)|

一目了然，这是k++后的行列式的样子。（pascal同学早日转C++）

那么立即推出：

|F(k - 1) H(k - 1)| = |F(2) H(2)|

|F(k - 2) H(k - 2)| |F(1) H(1)|

FH定理得证。

利用FH定理， 把F(1) = -1，F(2) = -3, H(1) = 3, H(2) = 8带入：

|F(n - 1) H(n - 1)| = -1

|F(n - 2) H(n - 2)|

于是就爽了嘛！

A = H(n) + F(n - 1) + (-1) - 1

化简：

A = H(n) + F(n - 1) - 2

进一步我们发现……

设R(n) = H(n) + F(n - 1) - 2

那么立即有：

R(n) = 3 * H(n - 1) - H(n - 2) + 3 * F(n - 2) - F(n - 3) - 2

= 3 * (R(n - 1) + 2) - (R(n - 2) + 2) - 2

R(n) = 3 * R(n - 1) - R(n - 2) + 2

所以，得到如下结论：

【轮状病毒定理】

轮状病毒的方案数满足递推式F(n) = 3 * F(n - 1) - F(n - 2) + 2， 其中F(1) = 1, F(2) = 5

#include<string>
#include<iostream>
#include<iosfwd>
#include<cmath>
#include<cstring>
#include<stdlib.h>
#include<stdio.h>
#include<cstring>
#define MAX_L 2005 //最大长度，可以修改
using namespace std;

class bign
{
public:
    int len, s[MAX_L];//数的长度，记录数组
    //构造函数
    bign();
    bign(const char*);
    bign(int);
    bool sign;//符号 1正数 0负数
    string toStr() const;//转化为字符串，主要是便于输出
    friend istream& operator>>(istream &,bign &);//重载输入流
    friend ostream& operator<<(ostream &,bign &);//重载输出流
    //重载复制
    bign operator=(const char*);
    bign operator=(int);
    bign operator=(const string);
    //重载各种比较
    bool operator>(const bign &) const;
    bool operator>=(const bign &) const;
    bool operator<(const bign &) const;
    bool operator<=(const bign &) const;
    bool operator==(const bign &) const;
    bool operator!=(const bign &) const;
    //重载四则运算
    bign operator+(const bign &) const;
    bign operator++();
    bign operator++(int);
    bign operator+=(const bign&);
    bign operator-(const bign &) const;
    bign operator--();
    bign operator--(int);
    bign operator-=(const bign&);
    bign operator*(const bign &)const;
    bign operator*(const int num)const;
    bign operator*=(const bign&);
    bign operator/(const bign&)const;
    bign operator/=(const bign&);
    //四则运算的衍生运算
    bign operator%(const bign&)const;//取模（余数）
    bign factorial()const;//阶乘
    bign Sqrt()const;//整数开根（向下取整）
    bign pow(const bign&)const;//次方
    //一些乱乱的函数
    void clean();
    ~bign();
};
#define max(a,b) a>b ? a : b
#define min(a,b) a<b ? a : b

bign::bign()
{
    memset(s, 0, sizeof(s));
    len = 1;
    sign = 1;
}

bign::bign(const char *num)
{
    *this = num;
}

bign::bign(int num)
{
    *this = num;
}

string bign::toStr() const
{
    string res;
    res = "";
    for (int i = 0; i < len; i++)
        res = (char)(s[i] + '0') + res;
    if (res == "")
        res = "0";
    if (!sign&&res != "0")
        res = "-" + res;
    return res;
}

istream &operator>>(istream &in, bign &num)
{
    string str;
    in>>str;
    num=str;
    return in;
}

ostream &operator<<(ostream &out, bign &num)
{
    out<<num.toStr();
    return out;
}

bign bign::operator=(const char *num)
{
    memset(s, 0, sizeof(s));
    char a[MAX_L] = "";
    if (num[0] != '-')
        strcpy(a, num);
    else
        for (int i = 1; i < strlen(num); i++)
            a[i - 1] = num[i];
    sign = !(num[0] == '-');
    len = strlen(a);
    for (int i = 0; i < strlen(a); i++)
        s[i] = a[len - i - 1] - 48;
    return *this;
}

bign bign::operator=(int num)
{
    char temp[MAX_L];
    sprintf(temp, "%d", num);
    *this = temp;
    return *this;
}

bign bign::operator=(const string num)
{
    const char *tmp;
    tmp = num.c_str();
    *this = tmp;
    return *this;
}

bool bign::operator<(const bign &num) const
{
    if (sign^num.sign)
        return num.sign;
    if (len != num.len)
        return len < num.len;
    for (int i = len - 1; i >= 0; i--)
        if (s[i] != num.s[i])
            return sign ? (s[i] < num.s[i]) : (!(s[i] < num.s[i]));
    return !sign;
}

bool bign::operator>(const bign&num)const
{
    return num < *this;
}

bool bign::operator<=(const bign&num)const
{
    return !(*this>num);
}

bool bign::operator>=(const bign&num)const
{
    return !(*this<num);
}

bool bign::operator!=(const bign&num)const
{
    return *this > num || *this < num;
}

bool bign::operator==(const bign&num)const
{
    return !(num != *this);
}

bign bign::operator+(const bign &num) const
{
    if (sign^num.sign)
    {
        bign tmp = sign ? num : *this;
        tmp.sign = 1;
        return sign ? *this - tmp : num - tmp;
    }
    bign result;
    result.len = 0;
    int temp = 0;
    for (int i = 0; temp || i < (max(len, num.len)); i++)
    {
        int t = s[i] + num.s[i] + temp;
        result.s[result.len++] = t % 10;
        temp = t / 10;
    }
    result.sign = sign;
    return result;
}

bign bign::operator++()
{
    *this = *this + 1;
    return *this;
}

bign bign::operator++(int)
{
    bign old = *this;
    ++(*this);
    return old;
}

bign bign::operator+=(const bign &num)
{
    *this = *this + num;
    return *this;
}

bign bign::operator-(const bign &num) const
{
    bign b=num,a=*this;
    if (!num.sign && !sign)
    {
        b.sign=1;
        a.sign=1;
        return b-a;
    }
    if (!b.sign)
    {
        b.sign=1;
        return a+b;
    }
    if (!a.sign)
    {
        a.sign=1;
        b=bign(0)-(a+b);
        return b;
    }
    if (a<b)
    {
        bign c=(b-a);
        c.sign=false;
        return c;
    }
    bign result;
    result.len = 0;
    for (int i = 0, g = 0; i < a.len; i++)
    {
        int x = a.s[i] - g;
        if (i < b.len) x -= b.s[i];
        if (x >= 0) g = 0;
        else
        {
            g = 1;
            x += 10;
        }
        result.s[result.len++] = x;
    }
    result.clean();
    return result;
}

bign bign::operator * (const bign &num)const
{
    bign result;
    result.len = len + num.len;
    
    for (int i = 0; i < len; i++)
        for (int j = 0; j < num.len; j++)
            result.s[i + j] += s[i] * num.s[j];
    
    for (int i = 0; i < result.len; i++)
    {
        result.s[i + 1] += result.s[i] / 10;
        result.s[i] %= 10;
    }
    result.clean();
    result.sign = !(sign^num.sign);
    return result;
}

bign bign::operator*(const int num)const
{
    bign x = num;
    bign z = *this;
    return x*z;
}
bign bign::operator*=(const bign&num)
{
    *this = *this * num;
    return *this;
}

bign bign::operator /(const bign&num)const
{
    bign ans;
    ans.len = len - num.len + 1;
    if (ans.len < 0)
    {
        ans.len = 1;
        return ans;
    }
    
    bign divisor = *this, divid = num;
    divisor.sign = divid.sign = 1;
    int k = ans.len - 1;
    int j = len - 1;
    while (k >= 0)
    {
        while (divisor.s[j] == 0) j--;
        if (k > j) k = j;
        char z[MAX_L];
        memset(z, 0, sizeof(z));
        for (int i = j; i >= k; i--)
            z[j - i] = divisor.s[i] + '0';
        bign dividend = z;
        if (dividend < divid) { k--; continue; }
        int key = 0;
        while (divid*key <= dividend) key++;
        key--;
        ans.s[k] = key;
        bign temp = divid*key;
        for (int i = 0; i < k; i++)
            temp = temp * 10;
        divisor = divisor - temp;
        k--;
    }
    ans.clean();
    ans.sign = !(sign^num.sign);
    return ans;
}

bign bign::operator/=(const bign&num)
{
    *this = *this / num;
    return *this;
}

bign bign::operator%(const bign& num)const
{
    bign a = *this, b = num;
    a.sign = b.sign = 1;
    bign result, temp = a / b*b;
    result = a - temp;
    result.sign = sign;
    return result;
}

bign bign::pow(const bign& num)const
{
    bign result = 1;
    for (bign i = 0; i < num; i++)
        result = result*(*this);
    return result;
}

bign bign::factorial()const
{
    bign result = 1;
    for (bign i = 1; i <= *this; i++)
        result *= i;
    return result;
}

void bign::clean()
{
    if (len == 0) len++;
    while (len > 1 && s[len - 1] == '\0')
        len--;
}

bign bign::Sqrt()const
{
    if(*this<0)return -1;
    if(*this<=1)return *this;
    bign l=0,r=*this,mid;
    while(r-l>1)
    {
        mid=(l+r)/2;
        if(mid*mid>*this)
            r=mid;
        else
            l=mid;
    }
    return l;
}

bign::~bign()
{
}
/**********************上面是大整数部分*********************/
bign a[110];

int main() {
    a[1]=1,a[2]=5;
    int n;
    cin >> n;
    for(int i=3;i<=n;i++)
    a[i]=a[i-1]*3-a[i-2]+2;
    cout << a[n] << endl;
    return 0;
}