HDU 4507 吉哥系列故事——恨7不成妻 数位DP--思路分析

本篇博客只涉及我对这个题目的理解和公式方式的推导，没有完整的代码, 本篇一共三个方法分析，第一个是略写，后两个是重点，想看正解思路直接跳第三个方法就好。
题目要求出L，R之间与7无关的数字的平方和，而且范围在1e18之间，很容易想到是数位DP。这种时候就出现了不同的思路。
第一种是直接往爆搜靠拢，pos=0时，说明数据合理，return 这个数字的平方值。这个方法在我看来肯定不行，就算结果算对了，每一个符合条件的数字也至少会计算一次平方值，感觉会T，话说，这就是爆搜+记忆化了啊。

第二种是稍好一点，数位DP的思路已经理解了，不是最后计算平方值，而是一边往下搜索，一边维护当前平方值，每次加一位数字之后，就维护新的平方值，同时保证不同数字的后续状态能重合，这样就能在搜索到不同的数字时，直接return，而不是再搜索一次，有效的节约时间。

我们来分析一下（我的WA做法）这种做法的过程，也便于理解最后的正解。
对于(a+b+c)^2,我们可以写成
a * a
b * a + b * a + b * b
c * a + c * a + c * b +c * b + c * c
注意到了么，每个数字乘以的数字都是都是不超过自己的，对于每次增加的一个新数字，我们的新平方和sqsum=原sqsum + c * c +c * sum(数字a+b之和)
我们把这个结论应用在数位DP维护平方和过程中，那么就可以一边增加新的数字，一边维护平方和。只要把公式变成(a * 100 + b * 10 +c)^2 这样就可以很好契合，加入d之后就变成 (a * 1000 + b * 100 +c * 10 +d) ^ 2 = (a * 100+b * 10+c)^2 * 100 + d * (a * 1000 + b * 100 +c * 10) * 2 + d * d
我们用s_sum维护前缀和，当插入c 时，s_sum=a * 10+b
插入c之后，s_sum=a * 100+b*10+c
这样就能通过s_sum来一直维护平方和
然后用sum维护数位和%7, sta维护搜索的数字%7，这样就能在最后直接判断这个数字是不是0来判断是不是7的倍数。
同时，dp数组就会变成不超过[ len ][ 7 ][ 7 ]的大小，相对于1e18个数字，这个范围很小，也就是非常多的数据都会重合，减少了非常多的重复搜索。
下面的写法也能通过样例，似乎这就是正解。

//pos位置 sum数位之和 sta余数 limit是否达到上限 cal记录当前平方值  s_sum记录当前数字是多少
LL dfs(int pos,int sum,int sta,bool limit, LL cal,LL s_sum)
{
    if(pos==0) {
        if(sum==0||sta==0)
            return 0;
        return cal;
     }
    if(!limit&&dp[pos][sum][sta]!=-1) return dp[pos][sum][sta];
    int up=limit?a[pos]:9;
    LL res=0;
    for(int i=0;i<=up;i++){
        if(i==7) continue;
        res=(  res + dfs(pos-1, (sum+i)%7, (sta*10+i)%7, limit&&i==a[pos] ,(cal*100%mod+i*i+s_sum*10*i*2%mod)%mod ,(s_sum*10+i)%mod)%mod  )%mod;//第一个sum是dp要用的 第二个s_sum是记录当前数字的值
    }
    if(!limit) dp[pos][sum][sta]=res;
    return res;
}

然后它WA了，WA的十分漂亮，原因就在于，这个dp数组并没有真正做到状态合并，dp数组维护的是平方值，但是两个长度相同的数字，在数位和%7相同 和 大小%7相同的情况下，并不能继承相同的平方值，比如12和82。数位DP的这个模板是从高位往低位开始搜索，并且数字从小到大，先搜索00000~ 09999，再搜索10000~19999,依次类推，也就是当我们搜索到82时，这个dp[pos][sum][sta]已经更新过，也就是会直接return ，不会继续搜索，我们明明是82开始，却直接得到了小于等于12的不为7的数字的平方和，没有开展搜索，这肯定不行。
也就是说，我们DP数组维护的状态混淆了，或者，我们DP数组根本就不应该维护平方和。
我们需要想想，12和82到底相同的是什么。

这也就引出了我们的第三种思路，也就是正解。
12和82的数位和%7相同，大小%7相同，我们不妨把数字扩展一下，分别记作 1A, 8B, 其中A和B可以是多个数字，比如156， A就代表了56。在数位和 与 大小都相等的情况下， 我们不难想到A与B时等价的，任何一个能够使1A不为7的子集A，都在B中存在，反之亦然。也就是说，对于1A和8B，1和8的后继子集是完全相同的。
那么，我们能不能先把这些后继子集，也就是比较小的数字的平方和维护出来，然后根据每次增加一个高位来维护平方和。毕竟第二种方法是确定高位，每次新增加一个低位然后维护平方和。
可以的。
我们把(a+b+c+d)^2的平方再拆开看一次
a * a
b * a + b * a + b * b
c * a + c * a + c * b +c * b + c * c
d * a + d * a + d * b +d * b + d * c + d * c +d * d
其中a是高位,d是低位，也就是我们会先求出d * d的值，然后开始计算 (c+d)^2
( c + d ) ^ 2 =
d * d
c * c + c * d + c * d
不难发现这和法二的公式相比也就是反了一下，只需要我们用sqsum来表示平方和,sum表示数字的和就也可以维护。

现在开始详细的推导
d^2 =
d * d
(c * 10+ d )^2 =
d * d
c * 10 * d + c * 10 * d +c * 10 * c * 10
(b * 100 + c * 10 +d) ^ 2=
d * d
c * 10 * d + c * 10 * d +c * 10 * c * 10
b * 100 * d + b * 100 * d +b * 100 * c * 10 + b * 100 * c * 10 + b * 100 * b * 100
上下后不难得到，增加高位之后 sqsum = 低位的sqsum + (新增加的数字) b * 权值 * b * 权值 +b * sum（数字之和 10*c+d ）
只要我们能维护 b对应的权值（根据位置直接判断即可，预处理 10的幂次），数字之和就可以完成维护。
*最后，对于一个确定的高位b， 可能有很多个c ,d的组合，对于每一个组合，我们平方和增加的量对应的公式都是之和它们的和有关，所以，不妨直接维护所有数字的前缀和和符合条件的子集个数，得到sqsum=子集sqsum + 所有数字的前缀和 * (新增的数字 * 权值 ） + 子集个数 * (新增数字 * 权值)^2

那么现在和子集相关的三个属性，前缀和，符合条件的个数，平方和 都需要维护然后一并返回，不妨直接用结构体。放代码：代码来源：感谢大佬Orz

struct Node
{
    long long cnt;//与7无关的数的个数
    long long sum;//与7无关的数的和
    long long sqsum;//平方和
}dp[20][10][10];//分别是处理的数位、数字和%7,数%7
int bit[20];
long long p[20];//p[i]=10^i

Node dfs(int pos,int pre1,int pre2,bool flag)
{
    if(pos==-1)
    {
        Node tmp;
        tmp.cnt=(pre1!=0 && pre2!=0);//非1即0
        tmp.sum=tmp.sqsum=0;//下一位没有数字了 也就是子节点的sum总和=0
        return tmp;
    }
    if(!flag && dp[pos][pre1][pre2].cnt!=-1)
        return dp[pos][pre1][pre2];
    int end=flag?bit[pos]:9;
    Node ans;//记录子集总和
    Node tmp;//记录单次的搜索结果
    ans.cnt=ans.sqsum=ans.sum=0;
    for(int i=0;i<=end;i++)
    {
        if(i==7)continue;
        tmp=dfs(pos-1,(pre1+i)%7,(pre2*10+i)%7,flag&&i==end);
        ans.cnt+=tmp.cnt;//子集所有符合要求数字的个数
        ans.cnt%=MOD;
        ans.sum+=(tmp.sum+ ((p[pos]*i)%MOD)*tmp.cnt%MOD )%MOD;
        ans.sum%=MOD;

        ans.sqsum+=(tmp.sqsum + ( (2*p[pos]*i)%MOD )*tmp.sum)%MOD;
        ans.sqsum%=MOD;
        ans.sqsum+=( (tmp.cnt*p[pos])%MOD*p[pos]%MOD*i*i%MOD );
        ans.sqsum%=MOD;
    }
    if(!flag)dp[pos][pre1][pre2]=ans;
    return ans;
}

这样，这个问题就可以圆润的解决了。
比较细小的地方可以自己慢慢理解，数位DP的细节非常多，一个地方写错说不定答案就跑到哪里了，比如up写错啊，litmit的维护啊，sum的维护啊，代码不长但是变量很多，这些都需要慢慢熟悉了，加油加油。