I.什么是生成函数

生成函数（母函数）是一种将无限序列“打包”成一个形式幂级数的工具。比如对于序列 $a_0,a_1,a_2,\ldots$ 其普通生成函数定义为：

$$A(t)=\sum\limits_{k=0}^{\infty}a_kt^k$$

这里的 $t$ 是一个形式变量，我们不关心它的收敛性，只关心它的代数运算。通过研究生成函数的封闭形式（如有理分式），可以快速求出序列的通项或远项系数。

II.什么是 Bostan‑Mori 算法

Bostan - Mori 算法是一种快速求解有理分式远项系数的经典算法。它可以在 $O(k^2\log n)$ 的时间内求出 $[t^n]\frac{P(t)}{Q(t)}$（$k$ 为多项式的度数）。

III.代数代换

题目要求：

$$f(n)=\sum\limits_{x\in S_n}\sum\limits_{y\in S_n\wedge x<y}xy$$

利用对称型可得：

$$f(n)=\frac{1}{2}((\sum\limits_{x\in S_n}x)^2-\sum\limits_{x\in S_n}x^2)$$

因此我们只需要计算两个值：

$$A_n=(\sum\limits_{x\in S_n}x)^2$$ $$B_n=\sum\limits_{x\in S_n}x^2$$

答案为：

$$\frac{1}{2}(A_n-B_n)$$

建立计数生成函数、和生成函数、平方和生成函数：

$$C(t)=\sum\limits_{x\in S}t^{g(x)}$$ $$G(t)=\sum\limits_{x\in S}xt^{g(x)}$$ $$H(t)=\sum\limits_{x\in S}x^2t^{g(x)}$$

定义三个基本多项式，分别代表 数位和的贡献，数值的贡献，平方的贡献：

$$P(t)=\sum\limits_{d=1}^9 t^d$$ $$Q(t)=\sum\limits_{d=1}^9 d t^d$$ $$R(t)=\sum\limits_{d=1}^9 d^2 t^d$$

计数生成函数 $C(t)$

考虑对于一个合法数字，共有两种情况：

• 单独一位 $d$，贡献 $t^d$；
• 在已知的合法数字后添加一位 $d$，贡献 $t^{g(y)+d}$。

所以：

$$C(t)=\sum\limits_{d=1}^9 t^d+\sum\limits_{d=1}^9\sum\limits_{y\in S} t^{g(y)+d} \\ =\sum\limits_{d=1}^9 t^d+\sum\limits_{d=1}^9\sum\limits_{y\in S} t^{g(y)}\cdot t^d \\ =\sum\limits_{d=1}^9 t^d+\sum\limits_{d=1}^9(t^d\cdot\sum\limits_{y\in S} t^{g(y)}) \\ =\sum\limits_{d=1}^9 t^d+\sum\limits_{d=1}^9t^d(\sum\limits_{y\in S} t^{g(y)}) \\ =P(t)+P(t)C(t) \\ C(t)+1=\frac{1}{1-P(t)}$$

和生成函数 $G(t)$

相同的分类：

• 单独一位 $d$，贡献 $dt^d$；
• 在已知的合法数字后添加一位 $d$，$x$ 变为 $10y+d$，贡献 $(10y+d)t^{g(y)+d}=10yt^{g(y)}t^d+dt^{g(y)}t^d$。

所以：

$$G(t)=\sum\limits_{d=1}^9dt^d+\sum\limits_{d=1}^9(\sum\limits_{y\in S}10yt^{g(y)}t^d+\sum\limits_{y\in S}dt^{g(y)}t^d) \\ =\sum\limits_{d=1}^9dt^d+\sum\limits_{d=1}^9(10t^d\sum\limits_{y\in S}yt^{g(y)}+dt^d\sum\limits_{y\in S}t^{g(y)}) \\ =\sum\limits_{d=1}^9dt^d+10\sum\limits_{d=1}^9t^d\sum\limits_{y\in S}yt^{g(y)}+\sum\limits_{d=1}^9dt^d\sum\limits_{y\in S}t^{g(y)} \\ =Q(t)+10P(t)G(t)+Q(t)C(t) \\ G(t)(1-10P(t))=Q(t)(C(t)+1)$$

由计数生成函数 $C(t)+1=\frac{1}{1-P(t)}$ 得：

$$G(t)(1-10P(t))=Q(t)(\frac{1}{1-P(t)}) \\ G(t)=\frac{Q(t)}{(1-10P(t))(1-P(t))}$$

平方和生成函数 $H(t)$

同理：

• 单独一位 $d$，贡献 $d^2t^d$；
• 在已知的合法数字后添加一位 $d$，$x^2$ 变为 $(10y+d)^2=100y^2+20dy+d^2$，贡献为 $(100y^2+20dy+d^2)t^{g(y)+d}=100y^2t^{g(y)+d}+20dyt^{g(y)+d}+d^2t^{g(y)+d}$。

$$H(t)=\sum\limits_{d=1}^9d^2t^d+\sum\limits_{d=1}^9(\sum\limits_{y\in S}100y^2t^{g(y)+d}+\sum\limits_{y\in S}20dyt^{g(y)+d}+\sum\limits_{y\in S}d^2t^{g(y)+d}) \\ =\sum\limits_{d=1}^9d^2t^d+\sum\limits_{d=1}^9(100t^d\sum\limits_{y\in S}y^2t^{g(y)}+20dt^d\sum\limits_{y\in S}yt^{g(y)}+d^2t^d\sum\limits_{y\in S}t^{g(y)}) \\ =\sum\limits_{d=1}^9d^2t^d+100\sum\limits_{d=1}^9t^d(\sum\limits_{y\in S}y^2t^{g(y)})+20\sum\limits_{d=1}^9dt^d(\sum\limits_{y\in S}yt^{g(y)})+\sum\limits_{d=1}^9d^2t^d(\sum\limits_{y\in S}t^{g(y)}) \\ =R(t)+100P(t)H(t)+20Q(t)G(t)+R(t)C(t) \\ H(t)(1-100P(t))=R(t)(C(t)+1)+20Q(t)G(t) \\ H(t)=\frac{R(t)(C(t)+1)+20Q(t)G(t)}{1-100P(t)}$$

由计数生成函数 $C(t)+1=\frac{1}{1-P(t)}$ 以及和生成函数 $G(t)=\frac{Q(t)}{(1-10P(t))(1-P(t))}$ 得：

$$H(t)=\frac{R(t)(\frac{1}{1-P(t)})+20Q(t)\frac{Q(t)}{(1-10P(t))(1-P(t))}}{1-100P(t)} \\ =\frac{R(t)(1-10P(t))+20Q(t)^2}{(1-P(t))(1-10P(t))(1-100P(t))}$$

计算前缀

因为：

$$\frac{1}{1-t}\sum\limits_{k\ge0}a_kt^k=\sum\limits_{n\ge0}(\sum\limits_{k=0}^na_k)t^n$$

所以定义 $G_{pre}(t)$ 和 $H_{pre}(t)$ 分别为 $G(t)$ 和 $H(t)$ 的前缀和，则：

$$G_{pre}(t)=\frac{G(t)}{1-t}=\frac{Q(t)}{(1-10P(t))(1-P(t))(1-t)} \\ H_{pre}(t)=\frac{R(t)(1-10P(t))+20Q(t)^2}{(1-P(t))(1-10P(t))(1-100P(t))(1-t)}$$

本题要求 $S_n=\{x\in S\mid g(x)\le n\}$ 的和与平方和，即：

$$A_n=\sum\limits_{k=0}^n[t^n]G(t)=[t^n]\sum\limits_{k=0}^nG(t)=[t^n]G_{pre}(t) \\ B_n=\sum\limits_{k=0}^n[t^n]H(t)=[t^n]\sum\limits_{k=0}^nH(t)=[t^n]H_{pre}(t)$$

所以只要求出这两个有理分式的第 $n$ 项系数，就能得到答案。

而使用 Bostan‑Mori 算法即可求解。

时间复杂度 $O(D^2\cdot\log n)$，其中 $D$ 为多项式最大度数。

IV.代码

::::info[code]

#include <bits/stdc++.h>
using namespace std;
const long long N=110,MOD=1000003;
long long n,inv2,A_n,B_n,ans;
long long qpow(long long a,long long p){
    long long res=1;
    while(p){
        if(p&1) res=res*a%MOD;
        a=a*a%MOD;
        p>>=1;
    }
    return res;
}
struct Pol{
    int c[N],sz;
    Pol(){
        sz=0;
        memset(c,0,sizeof(c));
    }
};
Pol P,Q,R,den1,den10,den100,dent,A_N,A_D,B_N,B_D;
Pol mul(const Pol &a,const Pol &b){
    Pol res;
    res.sz=a.sz+b.sz-1;
    for(int i=0;i<a.sz;i++){
        if(!a.c[i]) continue;
        for(int j=0;j<b.sz;j++){
            if(!b.c[j]) continue;
            res.c[i+j]=(res.c[i+j]+(long long)a.c[i]*b.c[j])%MOD;
        }
    }
    while(res.sz>1&&res.c[res.sz-1]==0) res.sz--;
    return res;
}
Pol add(const Pol &a,const Pol &b){
    Pol res;
    res.sz=max(a.sz,b.sz);
    for(int i=0;i<res.sz;i++){
        int x=i<a.sz?a.c[i]:0;
        int y=i<b.sz?b.c[i]:0;
        res.c[i]=(x+y)%MOD;
    }
    while(res.sz>1&&res.c[res.sz-1]==0) res.sz--;
    return res;
}
Pol mulc(const Pol &a,int k){
    Pol res=a;
    for(int i=0;i<a.sz;i++) 
        res.c[i]=(long long)a.c[i]*k%MOD;
    return res;
}
Pol even(const Pol &a){
    Pol res;
    for(int i=0;i<a.sz;i+=2) 
        res.c[res.sz++]=a.c[i];
    if(res.sz==0) 
        res.sz=1,res.c[0]=0;
    return res;
}
Pol odd(const Pol &a){
    Pol res;
    for(int i=1;i<a.sz;i+=2) 
        res.c[res.sz++]=a.c[i];
    if(res.sz==0)
        res.sz=1,res.c[0]=0;
    return res;
}
Pol constr(int k){
    Pol res;
    res.sz=10;res.c[0]=1;
    for(int i=1;i<=9;i++) 
        res.c[i]=(MOD-(long long)k*P.c[i]%MOD)%MOD;
    while(res.sz>1&&res.c[res.sz-1]==0) res.sz--;
    return res;
}
long long bostanmori(Pol N,Pol D,long long n){
    while(n){
        Pol D_neg=D;
        for(int i=1;i<D_neg.sz;i+=2) 
            if(D_neg.c[i]) 
                D_neg.c[i]=(MOD-D_neg.c[i])%MOD;
        Pol S=mul(N,D_neg);
        if(n&1) N=odd(S);
        else N=even(S);
        Pol T=mul(D,D_neg);
        D=even(T);
        n>>=1;
    }
    return (long long)N.c[0]*qpow(D.c[0],MOD-2)%MOD;
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    inv2=(MOD+1)/2;
    cin>>n;
    P.sz=Q.sz=R.sz=10;
    for(int i=1;i<=9;i++) P.c[i]=1;
    for(int i=1;i<=9;i++) Q.c[i]=i;
    for(int i=1;i<=9;i++) R.c[i]=(long long)i*i%MOD;
    den1=constr(1);
    den10=constr(10);
    den100=constr(100);
    dent.sz=2;dent.c[0]=1;dent.c[1]=MOD-1;
    A_N=Q;
    A_D=mul(mul(den1,den10),dent);
    Pol tmp1=mul(R,den10);
    Pol tmp2=mul(Q,Q);
    B_N=add(tmp1,mulc(tmp2,20));
    B_D=mul(mul(mul(den1,den10),den100),dent);
    A_n=bostanmori(A_N,A_D,n),B_n=bostanmori(B_N,B_D,n);
    ans=((A_n*A_n%MOD-B_n+MOD)%MOD)*inv2%MOD;
    cout<<ans<<'\n';
    return 0;
}

::::

抢到最优解。