【简●解】POJ 1845 【Sumdiv】

【题目大意】#

给定AB,求AB的所有约数之和,对9901取模。
(对于全部数据,0<=A<=B<=50,000,000

【算法关键词】#

  • 数论
  • 综合模板
  • 二分,乘法逆元

【分析】#

不管什么题首先思考的肯定是暴力解法。起码可以骗分啊,当然,如果能一眼标算,那再好不过了。
这道题暴力做法就不说了,其实仔细思考也不会真的打暴力吧。。。

看见约数,首先想到的应该就是数论的有关知识。那么这道题实际上就是在求AB的约数之和,问题的难点就在于求解答案的时间优化。

首先,可以思考唯一分解定理,即任意一个自然数都可分解且只能分解成以下形式:

n=p1k1p2k2p3k3...pmkm

其中,p为质因数,k为自然数。

那如何实现呢?

首先通过线性筛把题目范围内的质数筛出来存在一个数组里,然后枚举数组里的质数是否能被A整除,即枚举p,如果能,就以p为底枚举k,并存入相对应的数组。

于是就可分解为:

AB=p1k1Bp2k2Bp3k3B...pmkmB

a所有约数和为ans(a)
对于两个互质的数a,b,必定有ans(ab)=ans(a)ans(b):

a约数为x1,x2...xmb约数为y1,y2,...yn
a,b互质,a1>amb2>bn无公共元素
那么

ans(ab)=i=1,j=1i<=m,j<=nxiyj=i=1mxij=1nyj=ans(a)ans(b)

得证 这不是积性函数的性质吗

同理,pikiBpjkjB互质
所以:

ans(AB)=i=1msumi

其中:

sumi=j=0kiBpij

即:

ans(AB)=i=1mj=0kiBpij

(理解不透彻的可以把这个拆开看看)

公式推到这里,其实就可以写代码了,在这里,还要注意sumi的计算。

sumi=pi0+pi1+pi2+...+piki

这里有两种方法。

  1. 二分+递归思想
    ki为奇数:

sum(p,k)=sum(p,k/2)(1+pn/2+1)

ki为偶数:

sum(p,k)=sum(p,k/21)(1+pn/2+1)+pn/2

证明的话可以把上面的式子拆开来看,这里不再赘述。

  1. 等比数列通项公式

sum(p,k)=pk+11p1

当然,如果在模的意义下这样做就需要乘法逆元的相关知识,这里就只给出法一的代码。法二的代码可以自己尝试下。其实是因为我懒

Copy
#include<cmath> #include<ctime> #include<queue> #include<cstdio> #include<cstdlib> #include<cstring> #include<iostream> #include<algorithm> #define debug() puts("FBI WARNING!") #define ll long long #define R register using namespace std; /*常量申明*/ const int MAX = 500000; const int BIG_SEVEN = 777777; const int MEIZI = 100000; const int mod = 9901; /*函数申明*/ void ola(); void test(); void divide(); void cls();void solve(); ll binary_pow(ll, ll); ll sum(ll, ll); ll read(); /*变量申明*/ ll a, b, s, prim[BIG_SEVEN], p[MEIZI], k[MEIZI]; bool vis[MAX + 5]; int cnt, tot; int main(){ ola(); while (~scanf("%lld %lld", &a, &b)) { cls(); divide(); solve(); //test(); } return 0; } inline ll read(){ // 快读 (貌似没用上) ll f = 1, x = 0;char ch; do { ch = getchar(); if (ch == '-') f = -1; } while (ch < '0'||ch>'9'); do {x = x*10+ch-'0'; ch = getchar(); } while (ch >= '0' && ch <= '9'); return f*x; } inline void ola() { // 欧拉筛 vis[1] = 1; for (R int i = 2;i < MAX; ++i) { if (!vis[i]) { prim[cnt++] = i; } for (R int j = 0;j < cnt ; ++j) { if (i*prim[j] > MAX) break; vis[i*prim[j]] = 1; if (i%prim[j] == 0) break; } } } inline void cls() {tot = 0, cnt = 0;} // 多组数据初始化 inline void divide() { // 唯一分解定理 memset(k, 0, sizeof(k));ll rose = a; for (R int i = 0; prim[i] <= rose/prim[i]; ++i) { if (rose%prim[i] == 0) { p[tot] = prim[i]; while (rose%prim[i] == 0) { // debug(); k[tot]++; rose /= prim[i]; } tot++; } } if (rose != 1) { p[tot] = rose; k[tot++] = 1; } } inline void solve() { // 累加答案 ll ans = 1; for (R int i = 0; i < tot; ++i) { ans *= (sum(p[i], b*k[i]) %mod); ans %= mod; } printf("%lld\n", ans); } inline void test() { // 测试函数 /*check ola()*/ for (R int i = 0;i < cnt; ++i) { cout << prim[i] << " "; } printf("%.2lf\n",(double)clock()/CLOCKS_PER_SEC); /*check divide()*/ for (R int i = 0;i < tot; ++i) { printf("%lld^%lld\n", p[i], k[i]); } cout << tot; } inline ll binary_pow(ll a, ll b) { // 二进制实现的快速模幂 ll ans = 1, tmp = a%mod; while (b) { if (b&1) {ans *= tmp; ans %= mod;} b >>= 1; tmp *= tmp; tmp %= mod; } return ans; } inline ll sum(ll p, ll k) { // 约数和二分实现(p^0+p^1+...+p^k) if (p == 0) return 0; if (k == 0) return 1; if (k & 1) return ((1+binary_pow(p, k/2+1)) %mod * sum(p, k/2) %mod) %mod; // 奇数 else return ((1+binary_pow(p, k/2+1)) %mod*sum(p, k/2-1) + binary_pow(p, k/2) %mod) %mod; }
posted @   SilentEAG  阅读(252)  评论(0编辑  收藏  举报
编辑推荐:
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
· 没有源码,如何修改代码逻辑?
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示
CONTENTS