A sequence of integer{an}$\left\{a_n\right\}$can be expressed as:

an=⎧⎩⎨⎪⎪⎪⎪023an−1−an−22+n+1n=0n=1n>1$$a_n=\{\begin{array}{cc}0& n=0\\ 2& n=1\\ \frac{3a_{n-1}-a_{n-2}}{2}+n+1& n>1\end{array}$$
Now there are two integers nn and mm. I’m a pretty girl. I want to find all b1,b2,b3⋯bp$b_1,b_2,b_3\cdots b_p$that1⩽bi⩽n$1⩽b_i⩽n$ is relatively-prime with the integer m$m$. And then calculate:
∑i

=1pabi$$\sum _{i=1}^p{a}_{b_i}$$
But I have no time to solve this problem because I am going to date my boyfriend soon. So can you help me?
Input
Input contains multiple test cases ( about 15000 ). Each case contains two integers n and m. 1≤n,m≤108$1\le n,m\le {10}^8$
.

Output
For each test case, print the answer of my question(after mod 1,000,000,007).

Hint
In the all integers from 1 to 4, 1 and 3 is relatively-prime with the integer 4. So the answer is a_1+a_3=14
樣例輸入 複製
4 4
樣例輸出 複製
14

題意

已知ai$a_i$序列，給你一個n和m求小於n與m互質的數作為a序列的下標的和

思路

我們通過打表發現ai$a_i$序列就等於 ai=i(i+1)$a_i=i(i+1)$，那麼原問題就可以轉換為求解

∑i=1ni(i+1)[gcd(i,m)==1]$$\sum _{i=1}^{n}i(i+1)[gcd(i,m)==1]$$
那我們可以考慮

[1,n]$[1,n]$中與m有相同因子的數
∑i=1ni(i+1)−∑d|n,d|md(d+1)$$\sum _{i=1}^{n}i(i+1)-\sum _{d|n,d|m}d(d+1)$$
又有i(i+1)=i2+i$i(i+1)=i^2+i$
∑i=1ni2+i=n(n+1)(2n+1)6+n(n+1)2$$\sum _{i=1}^n{i}^2+i=\frac{n(n+1)(2n+1)}{6}+\frac{n(n+1)}{2}$$
我們要是列舉m的素因子，那麼我們可以通過容斥找出[1,n]$[1,n]$中與m不互質的數，例如m的因子中含有2，那麼2,4,6,8,10⋯$2,4,6,8,10\cdots$都是與m不互質的數，這些數的和我們可以提出一個2，然後數的長度就是n2$\frac{n}{2}$，然後套用上面的求和公式即可，然後對不同的因子的貢獻再套用容斥模板即可，記得求和公式要用逆元

#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <string.h>
using namespace std;
const long long mod=1e9+7;
long long a[10005];
long long inv6=166666668;
long long inv2=500000004;
long long clac(long long n,long long i)
{
    n=n/i;
    return (n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod*i%mod*i%mod+n%mod*(n+1)%mod*inv2%mod*i%mod)%mod;
}
int main()
{
    long long n,m;
    while(scanf("%lld%lld",&n,&m)!=EOF)
    {
        int cnt=0;
        for(int i=2;i*i<=m;i++)
            if(m%i==0)
        {
            a[cnt++]=i;
            while(m%i==0)
                m/=i;
        }
        if(m!=1)
            a[cnt++]=m;
        long long ans=clac(n,1);
        long long ans2=0;
        for(int i=1;i<(1<<cnt);i++)
        {
            int flag=0;
            long long temp=1;
            for(int j=0;j<cnt;j++)
                if(i&(1<<j))
            {
                flag++;
                temp=temp*a[j]%mod;
            }
            temp=clac(n,temp);
            if(flag%2==1)
            {
                ans2=(ans2%mod+temp%mod)%mod;
            }
            else
                ans2=(ans2%mod-temp%mod+mod)%mod;
        }
        printf("%lld\n",(ans%mod-ans2%mod+mod)%mod);
    }
    return 0;
}