Problem 362 「無平方因数」 †
54という数字について考えよう.
54 は 1 より大きい因数を1つ以上使って7つの異なる方法で因数分解することができる.
54, 2×27, 3×18, 6×9, 3×3×6, 2×3×9 そして 2×3×3×3.
無平方(squarefree)の因数のみを使う場合は, 2つが残る: 3×3×6 と 2×3×3×3.
1 よりも大きい無平方の1つ以上の因数で n を因数分解する方法の個数を Fsf(n) と呼ぼう, つまり Fsf(54)=2 となる.
k=2 から n までの ΣFsf(k) を S(n) としよう.
S(100)=193 である.
S(10 000 000 000) を求めよ.
計算時間１時間１０分の解法を思いついてとりあえず答えを見ました。
愚直に全パターンを計算しましたが、条件を満たさない数の補集合から計算したほうが速そうです。
中学校レベルの教育を受けずにそだった身としてはこの辺で満足してます。
競技プログラムの問題なんてのはプロの技術屋から見たらレベルが低すぎるんです。
といわれたことがあるので技術屋ではない身としては技術屋ってすげえなと思うことしきり。
解法解説
つたない幼い解法ですが解説を。
Ｓ（100）の場合を考えます。
計算に必要なのは10までの素数と10から100までの素数です。
10までの素数は低い計算量で出ます。
10以上の素数は10以下の素数から求めます。
まず100倍できる数として１を考えます。
まず２を考えます。
1を2倍2倍、、、としていくと2倍でできる数は
(後ｎ倍できる、その条件を満たす数の個数)を考えます。
(100,1),(50,1),(25,1),(12,1),(6,1),(3,1),(1,1)
と集合が得られます。
次はこの集合に3倍3倍、、、を掛けていく。
5倍5倍と掛けていく。
7倍7倍と掛けていく。
２，３、5,7を素因数とした数が
(10,a1),(8,a2),(7,a3),,,(1,an)
と個数が求まります。
これの穂集合をとれば10以上の素数の個数があと何倍できるかという条件で集計できます。
ここで(9,b)がないので簡単な計算で補完します。
後はこの発想を繰り返し適用することで
10以上の素数*10以下の無平方数
10以下の無平方数の積
等を計算できます。
#include<stdio.h>
#include<map>
#include<vector>
#include<iostream>
#include<set>
const __int64 L=100000;
const __int64 LIMIT=L*L;
std::vector<__int64> ps,vec;
std::map<__int64,__int64> memo2;
std::set<__int64> sets;
//L=100の正しい答えans4=18049 ans2+ans3=33223 51272
void setPrime(){
for(int i=2;i<=L;i++){
bool isP=true;
for(int j=2;j*j<=i;j++){
if(i%j==0){
isP=false;
break;
}
}
if(isP==true){
ps.push_back(i);
sets.insert(i);
}
}
}
void set2(){
for(int i=2;i<=L;i++){
bool inOK=true;
for(int j=2;j*j<=i;j++){
__int64 a=j*j;
if(i<a)break;
if(i%a==0){
inOK=false;
break;
}
}
if(inOK==true){
vec.push_back(i);
}
}
}
int main(){
setPrime();
set2();
__int64 ans4=0,ans1=0,ans2=0,ans3=0;
std::map<__int64,__int64>::iterator it,it2;
std::map<__int64,__int64> memo,memo3;
memo[LIMIT]=1;
printf("a\n");
for(int i=0;i<ps.size();i++){
__int64 p=ps[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,ps.size());
std::map<__int64,__int64> next;
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
memo[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
}
printf("b\n");
printf("\naaaaaaaaaa\n");
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
it=memo.begin();
it2=it;
it2++;
for(;(*it).first<L;it++,it2++){
(*it).second=LIMIT/(*it).first-LIMIT/(*it2).first-(*it).second;
}
it=memo.begin();
it2=it;
it2++;
for(;(*it).first<L;it++,it2++){
while(((*it).first<L)&&(*it).first+1<(*it2).first){
__int64 a=(*it).first;
__int64 d=LIMIT/(a+1)-LIMIT/(a+2);
(*it).second-=d;
memo[a+1]=d;
it=memo.find(a+1);
it2=it;
it2++;
}
}
it=memo.lower_bound(L);
it2=memo.end();
memo.erase(it,it2);
printf("\nbbbbbbbb\n");
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
for(int i=0;i<ps.size();i++){
__int64 p=ps[i];
it=memo.end();
it--;
std::map<__int64,__int64> next;
if(i%500==0)printf("(%d/%d)",i,ps.size());
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
next[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
for(it=next.begin();it!=next.end();it++){
memo[(*it).first]-=(*it).second;
}
}
printf("\ncccccccccccc\n");
//ここでmemoにはL以上の素数の個数が入っている
printf("\nu");
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
//L以上の素数にvec[i]を掛けた一個の数を全て作る
memo3.clear();
for(int i=0;i<vec.size();i++){
__int64 p=vec[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,vec.size());
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
memo3[n/p]+=add;
if(it==memo.begin())break;//念のためのおまじない
it--;
}
}
std::cout<<"\ndddddddddddddd";
for(it=memo3.begin();it!=memo3.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
std::cout<<"\n";
for(it=memo3.begin();it!=memo3.end();it++){
memo[(*it).first]+=(*it).second;
}
printf("\nu");
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
memo3.clear();
//(L以上の素数*(1 or L以下の無平方数))*0個以上の無平方数を計算
for(int i=0;i<vec.size();i++){
__int64 p=vec[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,vec.size());
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
memo[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
}
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
printf("\ne\n");
for(it=memo.begin();it!=memo.end();it++){
ans2+=(*it).second;
}
memo.clear();
memo3.clear();
//L以下の素因数で作られたＬ以上の合成数を検討
memo[LIMIT]=1;
for(int i=0;i<ps.size();i++){
__int64 p=ps[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,vec.size());
std::map<__int64,__int64> next;
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
next[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
for(it=next.begin();it!=next.end();it++){
memo[(*it).first]+=(*it).second;
}
}
it=memo.lower_bound(L);
it2=memo.end();
memo.erase(it,it2);
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
}
for(int i=0;i<vec.size();i++){
__int64 p=vec[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,vec.size());
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
memo[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
}
printf("\nans3 sums\n");
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
ans3+=(*it).second;
}
//L以下の無平方合成数の因数分解を検討
//たぶんここのパートは正しい処理
memo.clear();
memo[LIMIT]=1;
printf("\ng\n");
for(int i=0;i<vec.size();i++){
__int64 p=vec[i];
it=memo.end();
it--;
if(i%500==0)printf("(%d/%d)",i,vec.size());
while(1){
__int64 n=(*it).first;
__int64 add=(*it).second;
if(n<p)break;
memo[n/p]+=add;
it=memo.find(n);
if(it==memo.begin())break;//念のためのおまじない
it--;
}
}
memo.erase(LIMIT);
for(it=memo.begin();it!=memo.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
ans4+=(*it).second;
}
std::cout<<"\nans="<<ans1<<" ans2="<<ans2<<" ans3="<<ans3<<" "<<ans2+ans3<<" " <<ans4<<"\n";
std::cout<<"ans="<<ans1+ans2+ans3+ans4;
}
答え
457895958010
最終更新：2015年11月27日 19:11
