题解：P16097 [ICPC 2019 NAIPC] Cost of Living

· · 题解

题意简述

每种商品从某年到下一年的价格，等于原价格乘当年的通货膨胀率，再乘该商品固定的修正系数。已知部分价格和通货膨胀率，回答指定商品在指定年份的价格是否唯一；若唯一，求出其值。

解题思路

对价格和各个乘数取自然对数，乘法关系就变为线性关系。价格是否唯一，等价于其对数是否唯一，因为指数函数是一一对应的。

将年份改为从 0 开始编号。设商品 j 的初始价格对数为 u_j，修正系数对数为 v_j，前 t 次通货膨胀率的对数和为 s_t，其中 s_0=0。连续展开题目的价格递推式，得到第 t 年的价格对数：

\ln p_{j,t}=u_j+tv_j+s_t

同一种商品每经过一年都会多乘一次修正系数，因此 v_j 的系数为 t；通货膨胀率仅与年份有关，因此所有商品共用 s_t。

把 u_j、v_j 和 s_1,\dots,s_{y-1} 作为未知量，共有 n=2c+y-1 个。每个已知价格直接提供上式对应的方程。若从第 t-1 年到第 t 年的通货膨胀率为 r_t，则还有一个方程：

s_t-s_{t-1}=\ln r_t

不能要求整个方程组具有唯一解。部分商品可能完全没有已知价格，但这不妨碍其他商品的某个价格唯一。甚至某个价格对应的线性组合唯一时，其中每个未知量仍可能不唯一。我们要判断的是 询问对应的线性组合能否由已知方程推出。

维护一个按主元列编号的阶梯形方程组。f[i] 存储第一个非零系数在第 i 列、且该系数已经归一化为 1 的方程。插入一个方程时，按列从小到大消去已有主元：

若所有系数都被消掉，原方程没有提供新的独立限制。题目输入的小数可能存在舍入误差，因此不使用这种相关方程的微小常数残差判无解。系数消元使用 long double，以 10^{-12} 判断系数是否为 0。

处理询问时，将目标价格对数的系数写入同样的临时数组，但常数项初始化为 0，然后用已有方程消元。假设消去的方程线性组合右侧为 d，目标表达式就变成「剩余线性组合加 d」。临时数组的常数项记录的是 -d。

若所有未知量系数都变为 0，目标表达式恒等于 d，答案就是 exp(-a[n+1])。反之，若遇到一个非零系数所在列没有主元，则目标表达式不在已知方程的行空间中，无法唯一确定。

后一结论也可以从自由变量理解：阶梯形方程组会把所有主元变量写成自由变量的线性组合。询问不能被完全消掉，说明代入这些表达式后，至少一个自由变量仍有非零系数。改变它并重新确定主元变量，可以保持全部已知方程成立，却改变询问值。所有原始参数都要求为正，而任意实数取指数都为正，所以这种自由变化没有额外的符号限制。

实现中，第 2j-1、2j 列分别保存 u_j、v_j，第 2c+t 列保存 s_t。输入的第 i 年对应 t=i-1；当 i=1 时不需要添加前缀变量，因为 s_0=0。reduce 同时用于插入方程和处理询问，返回首个无法消去的列；返回 n+1 表示系数全部消去。

设已知方程数为 m，总时间复杂度为 O((m+q)n^2)，空间复杂度为 O(n^2)，其中 n\le209、m\le yc+y-1。

参考代码

#include <bits/stdc++.h>
using namespace std;

using ld=long double;
const int N=214;
const ld eps=1e-12;
int n;
ld f[N][N],a[N];
int reduce()
{
    for(int i=1;i<=n;i++)
    {
        if(abs(a[i])<eps)continue;
        if(abs(f[i][i])<eps)return i;
        ld tmp=a[i];
        for(int j=i;j<=n+1;j++)a[j]-=tmp*f[i][j];
    }
    return n+1;
}
void add()
{
    int pos=reduce();
    if(pos<=n)
    {
        ld tmp=a[pos];
        for(int i=pos;i<=n+1;i++)f[pos][i]=a[i]/tmp;
    }
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int y,c,q;
    cin>>y>>c>>q;
    n=2*c+y-1;
    for(int i=1;i<y;i++)
    {
        ld x;
        cin>>x;
        if(x<0)continue;
        fill(a,a+n+2,0);
        a[2*c+i]=1;
        if(i>1)a[2*c+i-1]=-1;
        a[n+1]=log(x);
        add();
    }
    for(int i=1;i<=y;i++)
    {
        for(int j=1;j<=c;j++)
        {
            ld x;
            cin>>x;
            if(x<0)continue;
            fill(a,a+n+2,0);
            a[2*j-1]=1;
            a[2*j]=i-1;
            if(i>1)a[2*c+i-1]=1;
            a[n+1]=log(x);
            add();
        }
    }
    cout<<fixed<<setprecision(10);
    while(q--)
    {
        int x,y;
        cin>>x>>y;
        fill(a,a+n+2,0);
        a[2*x-1]=1;
        a[2*x]=y-1;
        if(y>1)a[2*c+y-1]=1;
        if(reduce()<=n)cout<<-1<<'\n';
        else cout<<exp(-a[n+1])<<'\n';
    }
    return 0;
}