51nod 1575 Gcd and Lcm

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#include<queue>
#include<algorithm>
#include<cmath>
#include<stack>
#include<ctime>
using namespace std;

const int maxn = 1E5 + 10;
typedef long long LL;
typedef unsigned int u32;
const u32 Mod = 3000007;
const LL INF = 1000000007;

u32 tot,tp,cnt,pri[maxn],g0[maxn],g1[maxn],g2[maxn],f[maxn],last[maxn]
	,g[maxn],h[maxn],sum[maxn],fp[maxn],fn[maxn],res[maxn],stk[maxn],p[maxn];
bool not_pri[maxn];

stack <int> s;
vector <u32> mi[maxn];
vector <pair<u32,u32> > v[Mod];

inline void Insert(u32 Num,u32 pos)
{
	s.push(Num % Mod);
	v[Num % Mod].push_back(make_pair(Num,pos));
}

inline u32 Search(u32 Num)
{
	int pos = Num % Mod;
	for (int i = 0; ; i++)
		if (v[pos][i].first == Num) return v[pos][i].second;
}

inline u32 Calc_h(u32 p,u32 k)
{
	return mi[p][k] * (2 * k + 1) * (mi[p][k] - mi[p][k - 1]) + mi[p][k - 1];
}

inline u32 Calc_1(const u32 &N)
{
	u32 A = N,B = N + 1;
	if (A & 1) B >>= 1; else A >>= 1;
	return A * B;
}

inline u32 Calc_2(const u32 &N)
{
	u32 A = N,B = N + 1,C = 2 * N + 1;
	if (A % 2 == 0) A /= 2;
	else if (B % 2 == 0) B /= 2; else C /= 2;
	if (A % 3 == 0) A /= 3;
	else if (B % 3 == 0) B /= 3; else C /= 3;
	return A * B * C;
}

void Calc_g(u32 *g,int k)
{
	for (int i = 1; i <= tp; i++)
		last[i] = 0,g[i] = !k ? stk[i] : (k == 1 ? Calc_1(stk[i]) : Calc_2(stk[i]));
	int pos = 1;
	for (int i = 1; i <= cnt; i++)
	{
		u32 now = p[i] * p[i];
		sum[i] = mi[p[i]][k] + sum[i - 1];
		while (stk[pos] < now) ++pos;
		for (int j = tp; j >= pos; j--)
		{
			u32 tmp; int np = Search(stk[j] / p[i]);
			if (p[i] > stk[np]) tmp = 1;
			else tmp = g[np] - (sum[i - 1] - sum[last[np]]);
			g[j] -= mi[p[i]][k] * tmp; last[j] = i;
		}
	}
	for (int i = tp; i; i--)
	{
		if (last[i] == cnt) continue;
		if (p[cnt] >= stk[i]) {g[i] = 1; continue;}
		g[i] -= (sum[cnt] - sum[last[i]]);
	}
}

inline u32 Get_f(u32 Num,u32 p,int np)
{
	if (Num < p) return 1;
	return 1 + sum[Search(Num)] - sum[np];
}

void Clear()
{
	u32 Max = 0; cnt = tp = 0;
	while (!s.empty()) v[s.top()].clear(),s.pop();
}

void Solve()
{
	u32 n; cin >> n; u32 Sqrt = sqrt(n);
	if (n == 1) {cout << 1 << endl; return;}
	for (int i = 1; ; i++)
		if (pri[i] > Sqrt) break; else p[++cnt] = pri[i];
	for (int i = 1,Last; i <= n; i = Last + 1)
		stk[++tp] = n / i,Last = n / stk[tp];
	reverse(stk + 1,stk + tp + 1);
	for (int i = 1; i <= tp; i++) Insert(stk[i],i),f[i] = 1;
	Calc_g(g0,0); Calc_g(g1,1); Calc_g(g2,2);
	//cerr << (double)(clock()) / CLOCKS_PER_SEC << endl;
	for (int i = 1; i <= tp; i++)
	{
		sum[i] = sum[i - 1]; g[i] = (u32)(3) * (g2[i] - g1[i]) + g0[i];
		if (stk[i] <= Sqrt && !not_pri[stk[i]]) sum[i] += Calc_h(stk[i],1);
	}
	int pos = tp; reverse(p + 1,p + cnt + 1);
	for (int i = 1; i <= cnt; i++)
	{
		int np = Search(p[i]); u32 now = p[i] * p[i];
		while (stk[pos] >= now) f[pos] += (sum[pos] - sum[np]),--pos;
		for (int j = tp; j > pos; j--)
			for (LL c = 1,k = p[i]; k <= stk[j]; k *= 1LL * p[i],c++)
			{
				u32 l = stk[j] / k;
				if (l >= now) f[j] += Calc_h(p[i],c) * f[Search(l)];
				else f[j] += Calc_h(p[i],c) * Get_f(l,p[i],np);
			}
	}
	//cerr << (double)(clock()) / CLOCKS_PER_SEC << endl;
	u32 Ans = 0;
	for (int i = 1; i <= Sqrt; i++)
		f[tp] -= h[i],Ans += h[i] * g[Search(n / i)];
	cout << Ans + f[tp] << endl; Clear();
}

int main()
{
	#ifdef DMC
		freopen("DMC.txt","r",stdin);
	#endif
	
	h[1] = res[1] = not_pri[1] = 1;
	for (int i = 2; i < maxn; i++)
	{
		if (!not_pri[i])
		{
			pri[++tot] = i; res[i] = fn[i] = 1; fp[i] = i;
			for (LL j = 1; j <= INF; j *= 1LL * i) mi[i].push_back(j);
		}
		h[i] = Calc_h(fp[i],fn[i]) * res[i];
		for (int j = 1; j <= tot; j++)
		{
			u32 Nex = i * pri[j];
			if (Nex >= maxn) break;
			not_pri[Nex] = 1;
			if (i % pri[j] == 0)
			{
				res[Nex] = res[i]; fp[Nex] = fp[i];
				fn[Nex] = fn[i] + 1; break;
			}
			res[Nex] = h[i]; fp[Nex] = pri[j]; fn[Nex] = 1;
		}
	}
	int T; cin >> T; while (T--) Solve();
	//cerr << (double)(clock()) / CLOCKS_PER_SEC << endl;
	return 0;
}

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值