【PE 479】Roots on the Rise

题目链接: https://projecteuler.net/problem=479

题意是说令\(a_k,b_k,c_k\)表示方程\(\frac 1 x = (\frac k x)^2(k+x^2)-k x\)的三个解(实数或复数)。

然后定义

$$S(n)=\sum_{p=1}^n\sum_{k=1}^n(a_k+b_k)^p(b_k+c_k)^p(c_k+a_k)^p$$

可以证明\(S(n)\)总是整数

求\(S(10^6)\)模\(1\,000\,000\,007\)的结果。


化简一下方程,即

$$kx^3-k^2x+x-k^3=0$$

其实很容易想到韦达定理,先列出来:

$$a_kb_kc_k=k^2$$

$$a_k+b_k+c_k=k$$

$$a_kb_k+b_kc_k+a_kc_k=\frac{1}{k}$$

那么式子里的\( (a_k+b_k)(b_k+c_k)(c_k+a_k)\)实际上利用上面三行配一下就行了,其实就是:

$$( a_kb_k+b_kc_k+a_kc_k )( a_k+b_k+c_k )- a_kb_kc_k$$

所以

$$ (a_k+b_k)(b_k+c_k)(c_k+a_k) =1-k^2$$

那么用交换求和,等比数列求和公式就可化简\(S(n)\)了:

$$ S(n)=\sum_{p=1}^n\sum_{k=1}^n(a_k+b_k)^p(b_k+c_k)^p(c_k+a_k)^p $$

$$ S(n)=\sum_{k=1}^n \sum_{p=1}^n [(a_k+b_k)(b_k+c_k)(c_k+a_k)]^p $$

$$ S(n)=\sum_{k=1}^n \sum_{p=1}^n (1-k^2)^p $$

$$ S(n)=\sum_{k=1}^n (1-k^2) \frac{1-(1-k^2)^n}{k^2} $$

然后就可以快速幂啦

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int n=1000000;
int power(int x,int n){
    int ans=1;
    while(n){
        if(n&1) ans=1LL*ans*x%mod;
        x=1LL*x*x%mod;
        n>>=1;
    }
    return ans;
}
int main(){
    int ans=0;
    for(int k=2;k<=n;++k){
        int k2=1LL*k*k%mod;
        ans=(ans+1LL*(1-k2)*(1-power(1-k2,n))%mod*power(k2,mod-2)%mod+mod)%mod;
    }
    printf("%d\n",ans);
}

【PE 243】Resilience

题目链接: https://projecteuler.net/problem=243

题意是说求最小的\(n\)满足

$$\frac{\varphi(n)}{n-1}<\frac{15499}{94744}$$

实际上这道题要从欧拉函数的公式入手

$$\frac{\varphi(n)}{n-1}=\frac{n}{n-1}\prod_{i=1}^{k}(1-\frac{1}{p_i})$$

要使这个比例尽可能小,就要选取尽可能多的不同素因子。那么就先取\(n=2\times 3\times 5\times 7 \times 11 \times \dots\)试一下,会发现 \(n=6469693230\)时正好满足题意,即前十个素数的乘积是一个解,那么就dfs一下这十个素数分别选几次,搜索的时候用这个上界限制一下即可。注意中间运算可能会爆long long。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll UP=6469693230LL;
const int N=100;
int prime[N+5],phi[N+5],tot;
bool vis[N+5];
int a[20];
ll ans;
void getprime(){
    phi[1]=1;
    for(int i=2;i<=N;++i){
        if(!vis[i]){
            prime[tot++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<tot && prime[j]<=N/i;++j){
            vis[i*prime[j]]=true;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }else {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }
}
void dfs(int x,ll n,int sta){
    if(x==10){
        ll p=n,q=n-1;
        for(int i=0;i<10;++i)
            if((sta>>i)&1)
                p*=prime[i]-1,q*=prime[i];
        if((long double)p/q<(long double)15499/94744 && n<ans){
            ans=n;
        }
        return;
    }
    ll power=1;
    dfs(x+1,n,sta);
    for(int i=1;i<=50 && power<=UP/prime[x] && n<=UP/(power*prime[x]);++i){
        a[x]=i;
        power*=prime[x];
        dfs(x+1,n*power,sta|(1<<x));
    }
}
int main(){
    getprime();
    ans=UP;
    dfs(0,1,0);
    cout<<ans<<endl;
}

【PE 201】Subsets with a unique sum

题目链接:  https://projecteuler.net/problem=201

题意是考虑集合\(S=\{1^2,2^2,\dots,100^2\}\)的所有大小为\(50\)的子集,对每个子集内元素分别求和,求这些和去重后的和是多少。

实际上就是一个二维背包,记\(f[i][j]\)表示物品个数为\(i\),和为\(j\)的方案数。最后累加方案数为\(1\)的状态即可。

之前知道背包和生成函数有一定联系,实际上二维背包也可以看成算多项式乘积。本质上就是求这个二元多项式中\(x\)指数为\(50\)的那些项的系数:

\(g(x,y)=(1+xy)(1+xy^4)(1+xy^9)\dots (1+xy^{10000})\)

#include<bits/stdc++.h>
using namespace std;
const int M=338350;
const long long mod=1e18;
long long f[51][M+5];
int main(){
    f[0][0]=1;
    for(int i=1;i<=100;++i)
        for(int j=50;j>=1;--j)
            for(int k=M;k>=i*i;--k){
                f[j][k]+=f[j-1][k-i*i];
                if(f[j][k]>=mod) f[j][k]-=mod;
            }
    long long ans=0;
    for(int i=1;i<=M;++i)
        if(f[50][i]==1)
            ans+=i;
    cout<<ans<<endl;
}

【PE 686】Powers of Two

题目链接:  https://projecteuler.net/problem=686

题意是定义了\(p(L,n)\)为第 \(n\) 小的 \(j\) ,满足 \(2^j\) 开头前几位是 \(L\) 。

比如 \(p(12,1)=7\),因为 \(2^7=128\)是第一个以 \(12\) 开头的。

现在要求 \(p(123,678910)\)。

因为是 \(n\) 不大,考虑直接从小到大枚举 \(j\) ,验证开头三位。这里用到的技巧和求斐波那契前几位一样,取log。

先用科学计数法表示一下:

\(2^j=s\times 10^t (s<10)\)

然后两边同时取log,以10为底:

\(j\log{2}=\log{s}+t ( \log{s} <1) \)

那么减去整数部分 \(t\)之后再用pow函数还原回去,就得到了 \(s\) 的值。当然只需要求前三位,所以double的有效位数就足够了。 \(s\) 乘 \(100\) 再截断取整即可。

#include<bits/stdc++.h>
using namespace std;
const double L=log10(2);
int main(){
    int j=10,n=0;
    while(true){
        double x=j*L-floor(j*L);
        double y=pow(10,x);
        if((int)(y*100)==123){
            ++n;
            if(n==45) printf("%d\n",j);
            if(n==678910) {
                printf("%d\n",j);
                break;
            }
        }
        ++j;
    }
}

【PE 388】Distinct Lines

题目链接: https://projecteuler.net/problem=388

题意是说,考虑三维空间中的整点 \((a,b,c)\)满足 \(0\le a,b,c\le N\) ,从原点\((0,0,0)\)到所有点连直线,记 \(D(N)\)表示不同的直线的个数,求 \(D(10^{10})\),提交答案的前九位加后九位。

二维的情况就是bzoj 2190,三维做法一样,详情看之前写的这篇数论初级魔法大赏。最后推出来答案为

\(\sum_{d=1}^{n}\mu(d)(\lfloor\frac{n}{d}\rfloor)^3+6\sum_{i=1}^{n}\varphi(i)\)。

直接开int128上杜教筛。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,l,r) for(int i=l;i<=r;++i)
const int N=2000000;
int prime[N+5],tot;
bool vis[N+5];
int mu[N+5],phi[N+5];
ll muS[N+5],phiS[N+5];
ll n;
unordered_map<ll,__int128> mp1,mp2;
void getprime()
{
    mu[1]=1;
    phi[1]=1;
    for(int i=2; i<=N; ++i)
    {
        if(!vis[i])
        {
            prime[tot++]=i;
            mu[i]=-1;
            phi[i]=i-1;
        }
        for(int j=0; j<tot && prime[j]<=N/i; ++j)
        {
            vis[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                mu[i*prime[j]]=0;
                break;
            }
            else
            {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
    for(int i=1; i<=N; ++i)
    {
        phiS[i]=phiS[i-1]+phi[i];
        muS[i]=muS[i-1]+mu[i];
    }
}
__int128 calc_muS(ll n)
{
    if(n<=N)
        return muS[n];
    if(mp1.count(n))
        return mp1[n];
    __int128 ans=1;
    for(ll l=2,r; l<=n; l=r+1)
    {
        r=n/(n/l);
        ans-=calc_muS(n/l)*(r-l+1);
    }
    return mp1[n]=ans;
}
__int128 calc_phiS(ll n)
{
    if(n<=N)
        return phiS[n];
    if(mp2.count(n))
        return mp2[n];
    __int128 ans=(__int128)n*(n+1)/2;
    for(ll l=2,r; l<=n; l=r+1)
    {
        r=n/(n/l);
        ans-=calc_phiS(n/l)*(r-l+1);
    }
    return mp2[n]=ans;
}
__int128 power(ll x,int n)
{
    __int128 res=1;
    while(n--)
        res*=x;
    return res;
}
void show(__int128 x)
{
    if(x==0)
        return;
    show(x/10);
    printf("%d",(int)(x%10));
}
int main()
{
    getprime();
    cin>>n;
    __int128 ans=0,pre=0,now=0;
    for(ll d=1; d<=n; ++d)
    {
        now=calc_muS(n/(n/d));
        ans+=power(n/d,3)*(now-pre);
        d=(n/(n/d));
        pre=now;
    }
    ans+=calc_phiS(n)*6;
    show(ans);
    puts("");
}

【PE 148】Exploring Pascal’s triangle

题目链接: https://projecteuler.net/problem=148

题意就是求杨辉三角前\(10^9\)行中不是7的倍数的数的个数。

因为7是素数,自然想到Lucas定理

\(C_n^{m} \equiv \prod_{i=0}^{k} C_{n_i}^{m_i} \pmod{p} \)

其中\(n=\sum_{i=0}^{k}n_{i}p^{i} ,m=\sum_{i=0}^{k}m_{i}p^{i}\)

所以即统计所有的\(C_n^{m}\),满足\(\forall{0\le i\le k}, m_i\le n_i\),这样把\(10^9\)转化成\(7\)进制,然后扫一下就能算出来了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int P=7;
ll n;
ll q[20];
int d[100],len;
void divide(int n){
    len=0;
    while(n){
        d[++len]=n%P;
        n/=P;
    }
}
ll work(){
    ll ans=0,pre=1;
    for(int i=len;i>=1;--i){
        ans+=pre*(d[i]*(d[i]+1)/2)*q[i-1];
        pre*=(d[i]+1);//0-d[i];
    }
    return ans;
}
int main(){
    q[0]=1;
    for(int i=1;i<20;++i)
        q[i]=q[i-1]*(P*(P+1)/2);
    while(scanf("%lld",&n)!=EOF){
        divide(n);
        printf("%lld\n",work());
    }
}
//input 1000000000
//Lucas

NOIP 2015

D1T1 神奇的幻方

按幻方的构造规则模拟即可。

#include<bits/stdc++.h>
using namespace std;
int n;
int a[50][50];
int main()
{
    scanf("%d",&n);
    int x=1,y=(n+1)/2;
    a[x][y]=1;
    for(int i=2;i<=n*n;i++)
    {
        if(x==1 && y!=n) x=n,y++;
        else if(x!=1 && y==n) x--,y=1;
        else if(x==1 && y==n) x++;
        else if(!a[x-1][y+1]) x--,y++;
        else x++;
        a[x][y]=i;
    }
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            printf("%d%c",a[i][j],(j<n)?' ':'\n');
}

D1T2 信息传递

考察图论知识,n个点每个点一条出边构成了一个基环树森林,答案就是求所有基环树环长的最小值。

#include<bits/stdc++.h>
using namespace std;
#define N 200010
#define rep(i,l,r) for(int i=l;i<=r;++i)
int n;
int a[N];
int fa[N];
int getfa(int x)
{
    return x==fa[x]?x:fa[x]=getfa(fa[x]);
}
bool Union(int x,int y)
{
    x=getfa(x);
    y=getfa(y);
    if(x==y) return false;
    fa[x]=y;
    return true;
}
int walk(int x)
{
    int ans=0,y=x;
    while(true)
    {
        y=a[y];
        ++ans;
        if(y==x) break;
    }
    return ans;
}
int main()
{
    scanf("%d",&n);
    rep(i,1,n) scanf("%d",a+i),fa[i]=i;
    int ans=~0U>>1;
    rep(i,1,n)
        if(!Union(i,a[i]))
            ans=min(ans,walk(i));
    printf("%d\n",ans);
}

D1T3 斗地主

搜索,首先要用最优性剪枝,然后就是考虑搜索顺序,肯定是要先走顺子这种出牌多的,最后考虑单张牌,对子这些。

#include<bits/stdc++.h>
using namespace std;
int T,n,ans;
int a[16];
int calc()
{
    int now=0;
    for(int i=0;i<=14;++i)
        if(a[i])
            ++now;
    return now;
}
void dfs(int step,int x)
{
    if(step>=ans) return;
    ans=min(ans,step+calc());
    for(int i=3;i<=13;++i)//san shun
    {
        int j=i;
        while(j<=14 && a[j]>=3) ++j;
        if(j-i<2) continue;
        a[i]-=3;
        for(int k=i+1;k<j;++k)
        {
            a[k]-=3;
            dfs(step+1,x-(k-i+1)*3);
        }
        for(int k=i;k<j;++k)
            a[k]+=3;
    }

    for(int i=3;i<=12;++i)// shuang shun
    {
        int j=i;
        while(j<=14 && a[j]>=2) ++j;
        if(j-i<3) continue;
        a[i]-=2;a[i+1]-=2;
        for(int k=i+2;k<j;++k)
        {
            a[k]-=2;
            dfs(step+1,x-(k-i+1)*2);
        }
        for(int k=i;k<j;++k)
            a[k]+=2;
    }
    for(int i=3;i<=10;++i) //dan shun
    {
        int j=i;
        while(j<=14 && a[j]>=1) ++j;
        if(j-i<5) continue;
        a[i]--;a[i+1]--;a[i+2]--;a[i+3]--;
        for(int k=i+4;k<j;++k)
        {
            a[k]--;
            dfs(step+1,x-(k-i+1));
        }
        for(int k=i;k<j;++k)
            a[k]++;
    }
    for(int i=2;i<=14;++i)
    if(a[i]==4)
    {
        a[i]=0;
        dfs(step+1,x-4);
        for(int j=0;j<=14;++j)
        if(a[j]>0)
        {
            a[j]--;
            for(int k=j;k<=14;++k)
            if(a[k]>0)
            {
                a[k]--;
                dfs(step+1,x-6);
                a[k]++;
            }
            a[j]++;
        }
        for(int j=2;j<=14;++j)
        if(a[j]>1)
        {
            a[j]-=2;
            for(int k=j;k<=14;++k)
            if(a[k]>1)
            {
                a[k]-=2;
                dfs(step+1,x-8);
                a[k]+=2;
            }
            a[j]+=2;
        }
        a[i]=4;
    }

    for(int i=2;i<=14;++i)
    if(a[i]>=3)
    {
        a[i]-=3;
        dfs(step+1,x-3);
        for(int j=0;j<=14;++j)
        if(a[j]>=1)
        {
            a[j]--;
            dfs(step+1,x-4);
            a[j]++;
        }
        for(int j=2;j<=14;++j)
        if(a[j]>=2)
        {
            a[j]-=2;
            dfs(step+1,x-5);
            a[j]+=2;
        }
        a[i]+=3;
    }
}
int main()
{
    scanf("%d%d",&T,&n);
    while(T--)
    {
        memset(a,0,sizeof(a));
        for(int i=1;i<=n;++i)
        {
            int x;
            scanf("%d%*d",&x);
            if(x==1) x=14;
            a[x]++;
        }
        ans=n;
        dfs(0,n);
        printf("%d\n",ans);
    }
}

D2T1 跳石头

比较裸的二分答案+贪心验证

#include<bits/stdc++.h>
using namespace std;
#define N 50010
#define rep(i,l,r) for(int i=l;i<=r;++i)
int a[N];
int n,m,H;
bool ok(int len)
{
    int last=0,i=0,ans=0;
    while(true)
    {
        while(i<=n&& a[i]-last<len) ++i;
        if(i>n) break;
        ++ans;
        last=a[i];
    }
    if(ans>0 && H-last<len) --ans;
    return n-ans<=m;
}

int main()
{
    scanf("%d%d%d",&H,&n,&m);
    rep(i,1,n) scanf("%d",a+i);
    int l=1,r=H;
    while(l<r)
    {
        int mid=l+(r-l+1)/2;
        if(ok(mid)) l=mid;
        else r=mid-1;
    }
    printf("%d\n",l);
}

D2T2 子串

一看就是dp,写了个暴力居然跑的飞快,改天再优化一下。

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
int n,m,K;
int dp[2][1005][205];
char a[1005],b[205];
void add(int &x,int y){
    x+=y;
    if(x>=mod) x-=mod;
}
int main(){
    scanf("%d%d%d",&n,&m,&K);
    scanf("%s",a+1);
    scanf("%s",b+1);
    for(int i=0;i<=n;++i)
        dp[0][i][0]=1;

    for(int k=1;k<=K;++k){
        memset(dp[k&1],0,sizeof(dp[k&1]));
        for(int i=1;i<=n;++i)
        for(int j=1;j<=m;++j)
        {
            add(dp[k&1][i][j],dp[k&1][i-1][j]);
            int l=0;
            while(l<j && l<i && a[i-l]==b[j-l]) {
                ++l;
                add(dp[k&1][i][j],dp[(k&1)^1][i-l][j-l]);
            }
            //cout<<i<<" "<<j<<" "<<k<<" "<<dp[i][j][k]<<endl;
        }
    }
    printf("%d\n",dp[K&1][n][m]);
}

D2T3 运输计划

给一棵带边权的树,再给m条路径,你可以选择一条边使其边权变为0,从而使得m条路径的长度的最大值最小,求最小为多少。

做法是二分答案,然后树上差分。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define N 300010
#define rep(i,l,r) for(int i=l;i<=r;++i)
#define dow(i,l,r) for(int i=l;i>=r;--i)
int e,head[N<<1],w[N<<1],last[N<<1],p[N];
int deep[N],fa[N],num[N],son[N];
int top[N];
int a[N],dis[N];
int u[N],v[N],lca[N],val[N];
int n,m;
void add(int x,int y,int c)
{
    head[++e]=y;w[e]=c;
    last[e]=p[x];
    p[x]=e;
}
void dfs(int x,int pre,int dep)
{
    deep[x]=dep;
    fa[x]=pre;
    num[x]=1;
    for(int j=p[x];j;j=last[j])
    {
        int y=head[j];
        if(y==pre) continue;
        dis[y]=dis[x]+w[j];
        dfs(y,x,dep+1);
        num[x]+=num[y];
        if(son[x]==-1||num[y]>num[son[x]])
            son[x]=y;
    }
}
void dfs_son(int x,int rt)
{
    top[x]=rt;
    if(son[x]==-1) return;
    else dfs_son(son[x],rt);
    for(int j=p[x];j;j=last[j])
        if(head[j]!=fa[x] && head[j]!=son[x])
            dfs_son(head[j],head[j]);
}
void dfs(int x)
{
    for(int j=p[x];j;j=last[j])
    {
        int y=head[j];
        if(y==fa[x]) continue;
        dfs(y);
        a[x]+=a[y];
    }
}
int calc(int x,int y)
{
    int f1=top[x],f2=top[y];
    while(f1!=f2)
    {
        if(deep[f1]<deep[f2])
        {
            swap(x,y);
            swap(f1,f2);
        }
        x=fa[f1];
        f1=top[x];
    }
    if(deep[x]>deep[y]) swap(x,y);
    return x;
}

bool ok(int mid)
{
    int tot=0,mx=0,ans=0;
    memset(a,0,sizeof(a));
    rep(i,1,m)
    if(val[i]>mid)
    {
        a[u[i]]++;
        a[v[i]]++;
        a[lca[i]]-=2;
        ++tot;
        ans=max(ans,val[i]);
    }
    dfs(1);
    rep(i,1,n)
        if(a[i]>=tot)
            mx=max(mx,dis[i]-dis[fa[i]]);
    return ans-mx<=mid;
}
int main()
{
    scanf("%d%d",&n,&m);
    rep(i,1,n-1)
    {
        int x,y,c;
        scanf("%d%d%d",&x,&y,&c);
        add(x,y,c);
        add(y,x,c);
    }
    memset(son,-1,sizeof(son));
    dfs(1,0,1);
    dfs_son(1,1);
    rep(i,1,m)
    {
        scanf("%d%d",u+i,v+i);
        lca[i]=calc(u[i],v[i]);
        val[i]=dis[u[i]]+dis[v[i]]-2*dis[lca[i]];
    }
    int l=0,r=4e8;
    while(l<r)
    {
        int mid=(l+r)>>1;
        if(ok(mid))  r=mid;
        else l=mid+1;
    }
    printf("%d\n",l);
    return 0;
}

还是改日吧

这几个月陆陆续续做了一些课件,然而都没时间写博客了。

昨天刚换了个服务器,为了方便,这次把wordpress放到docker里了。在迁移的时候,发现个坑就是有的插件可能跟不上wordpress更新的版本了,不再维护了这种。之前用的一个代码高亮的插件用不了了,于是今天把所有文章的代码块全部换了一遍,简直累死。。。

Burnside引理与Polya计数公式

Burnside引理

设\(A\)和\(B\)为有限集合,\(X\)表示所有从\(A\)到\(B\)的映射集合(着色方案集合)。\(X/G\)表示\(X\)所有元素的轨道的集合(去掉重复),即\(G\)作用在\(X\)上产生的所有等价类的集合,\(X^{g}\)表示\(X\)中在\(g\)作用下的不动元素,则有

\(|X/G|={\frac {1}{|G|}}\sum _{g\in G}|X^{g}|\)

\(X\)中非等价的着色数等于在\(G\)中的置换作用下保持不变的着色的平均数!

Polya计数公式

考虑置换\(g\)可以分解成若干不相交的循环置换的乘积,那么每个循环内的颜色必须相同,才能在\(g\)的作用下染色不变,设颜色一共\(m\)种,置换\(g\)的循环分解中的循环个数为\(c(g)\),那么在\(g\)作用下着色不变的着色数为:

\(|X^{g}|=m^{c(g)}\)

替换一下,就得到了:

\(|X/G|={\frac {1}{|G|}}\sum _{g\in G}m^{c(g)}\)

 

咕咕咕

转眼暑假就要过去,马上开学了,今年这个时间点格外的忙,本来这两月计划的总结也只能咕了,真棒!