24 の約数は 8 個ある, 1, 2, 3, 4, 6, 8, 12, 24.
100 以下で, ちょうど8個の約数を持っている数は 10 個ある, 24, 30, 40, 42, 54, 56, 66, 70, 78, 88.
ちょうど 8 個の約数を持つ n 以下の数の個数を f(n) としよう.
f(100) = 10, f(1000) = 180, そして f(10^6) = 224427.
f(10^12) を求めよ.
2013年購入の廉価パソコンでメモリを食わず計算時間30分で解。
(10^12の性質を悪用して少しコードを手抜きしているので他の数字では正しい答えが出ない可能性あり)
1兆は考えにくいので100でたとえます
考え方としては
c++言語的に
std::map<long long int,long long int> memo;
意味的には
memo[主キー]=個数
とし
主キーの意味を100を越えないであと最大何倍できる数が何個あるかと定義します。
1倍できるのは100~51
2倍できるのは50~34
3倍できるのは33~26
、、、
まず10までの素数を求めます、これは計算量が低いので力技で済みます。
memo[100]=1
と設定し、主キーmを整数演算で10以下の素数Pで割り、結果が1以上なら。
memo[m/p]+=memo[m]
とします。
するとmemoには10以下の素数でできた合成数がでます。
主キーの範囲でこの合成数の補集合をとると
memo[1~10]には、”10以上の素数”か”10以上の素数*10以下の合成数”の個数がでます。
ここからひと手間かかるのですが大事な点だけを言えば、”10以下の合成数*10以上の素数”でできた合成数を10以下の主キーから引いていけば10以上の素数の個数が主キーごとに分類されて出ます。
あと主キーには歯抜けができて歯抜けに分類される数は必ず素数ですのでその辺を調整します。
主キーとそれに対応した素数の個数が出れば後は簡単な計算で答えとなり”正答”しました。
#include<stdio.h>
#include<map>
#include<vector>
#include<iostream>
#include<algorithm>
const __int64 LIMITH=1000000;
const __int64 LIMIT=LIMITH*LIMITH;
bool isPrime(__int64 n){
for(__int64 i=2;i*i<=n;i+=(i&1)+1){
if(n%i==0)return false;
}
return true;
}
int main(){
std::vector<__int64> as;
//10^6までの素数を求める
for(__int64 i=2;i*i<=LIMIT;i++){
if(isPrime(i)==true)as.push_back(i);
}
printf("a1\n");
std::map<__int64,__int64> memo,memo2,memo3;
__int64 n=LIMIT;
memo[n]=1;
std::map<__int64,__int64>::iterator it,it2;
for(int i=as.size()-1;i>=0;i--){
it=memo.end();
it--;
__int64 m=as[i];
if(i%300==0)std::cout<<"("<<i<<"/"<<as.size()<<")";
while(1){
n=(*it).first;
if(n<m)break;
__int64 c=(*it).second;
memo[n/m]+=c;
it=memo.find(n);
if(it==memo.begin())break;
it--;
}
}
//for(it=memo.begin();it!=memo.end();it++){
// std::cout<<(*it).first<<" "<<(*it).second<<"\n";
//}
printf("a2\n");
it=memo.begin();
it2=it;
it2++;
for(;it2!=memo.end();it2++,it++){
if((*it).first>LIMITH)break;
memo2[(*it).first]=(LIMIT/(*it).first)-(LIMIT/(*it2).first)-(*it).second;
}
for(int i=1;i<LIMITH;i++){
if(memo2.find(i)==memo2.end()){
__int64 add=LIMIT/i-LIMIT/(i+1);
memo2[i]=add;
memo2[i-1]-=add;
}
}
//for(it=memo2.begin();it!=memo2.end();it++){
// if((*it).first*(*it).first>LIMIT)break;
// std::cout<<(*it).first<<" "<<(*it).second<<"\n";
//}
printf("a3\n");
for(int i=as.size()-1;i>=0;i--){
__int64 m=as[i];
std::map<__int64,__int64> next;
std::map<__int64,__int64>::reverse_iterator rIt;
if(i%300==0)std::cout<<"("<<i<<","<<as.size()<<")";
it=memo2.end();
while(it--,1){
__int64 n=(*it).first;
__int64 c=(*it).second;
if(n<m)break;
next[n/m]+=c;
if(it==memo2.begin())break;
}
for(it=next.begin();it!=next.end();it++){
__int64 f=(*it).first;
__int64 c=(*it).second;
memo2[f]-=c;
}
}
memo3.insert(memo2.begin(),memo2.end());
//for(it=memo2.begin();it!=memo2.end();it++){
// if((*it).first*(*it).first>LIMIT)break;
// std::cout<<(*it).first<<" "<<(*it).second<<"\n";
//}
__int64 sum=0,sum1=0;
std::map<__int64,__int64>::reverse_iterator rIt;
for(rIt=memo3.rbegin();rIt!=memo3.rend();rIt++){
sum+=(*rIt).second;
memo3[(*rIt).first]=sum;
}
//for(it=memo3.begin();it!=memo3.end();it++){
// if((*it).first>100)break;
// std::cout<<(*it).first<<" "<<(*it).second<<"\n";
//}
std::map<__int64,__int64> primeNo;
for(int i=0;i<as.size();i++){
primeNo[as[i]]=i+1;
}
__int64 maxPrime=as[as.size()-1];
//a*b*c
__int64 ans=0;
for(int i=0;i<as.size();i++){
if(i%300==0)std::cout<<"("<<i<<" "<<as.size()<<")";
for(int j=i+1;j<as.size();j++){
__int64 c=as[i]*as[j];
if(LIMIT<=c)break;
__int64 c1=LIMIT/c;
if(c1>as[j]){
it=primeNo.upper_bound(c1);
it--;
if((*it).first>as[j]){
ans+=(*it).second-primeNo[as[j]];
}
if(LIMIT>=c){
ans+=memo3[c];
}
}else{
break;
}
}
}
__int64 ans2=0;
for(int i=0;i<as.size();i++){
__int64 c=as[i]*as[i]*as[i];
if(LIMIT<=c)break;
__int64 c1=LIMIT/c;
it=primeNo.upper_bound(c1);
it--;
__int64 add=(*it).second;
ans2+=add-((*it).first>=as[i]);
if(LIMIT>=c){
ans2+=memo3[c];
}
}
__int64 ans3=0;
for(int i=0;i<as.size();i++){
__int64 a=as[i]*as[i]*as[i];
a*=a;
a*=as[i];
if(LIMIT<a)break;
ans3++;
}
std::cout<<"\nans="<<ans<<" "<<ans2<<" "<<ans3<<" "<<ans+ans2+ans3<<"\n";
}
最終更新:2016年02月26日 12:45