プロジェクトオイラー
Problem 521 「最小の素因数」 †
解法
ある数 n の最小の素因数を smpf(n) と定義しよう.
smpf(91)=7 となる, なぜなら 91=7×13, そして smpf(45)=3, なぜなら 45=3×3×5.
2 ≤ i ≤ n の範囲における smpf(i) の和を S(n) としよう.
たとえば, S(100)=1257.
S(10^12) mod 10^9 を求めよ.
とりあえず計算時間４５分で解けた。
本当は一分切らないといけないのだが
数学知識のない身としてはまあまあだと思う。
使った数学知識は素因数分解の一意性だけ。
以下数学知識０プログラムテクニックだけに頼った我流の解法
10^12は値が大きすぎるのでS(10000)で説明する。
まず10000を
Ｓ（10000）で説明すると
a=100以下の素数
b=aを素因数とする合成数
c=100以上の素数
d=(a xor b)c
と定義する。
a,bは1から初めてあと何倍できるかで10000から割っていくと
10000～5001までは倍数化できない。
5000～3334までは２倍はできる
、、、
101～101までは９９倍できる。
で後最大何倍できるかをキーに集約できる。
この手法により10000以下のa,b合成数を取り除くことができる。
101～101までは素数のみである。
99倍までできる。
後何倍できるかをキーに割り算で計算するとaとｂは求まる。
10000～101までの数からaとｂを引いていくとｃ∨ｄが得られる。
ｃは
101～101を99倍できるまでで101を素因数とする合成数をc∨ｄからｄを除外できる。
98倍９７倍、、、２倍までとキーをさげていきながら除算すると答えが出る。
わかりにくい説明ですいません。
頭の中ではわかってるのです。
ホワイトボードを後ろに口頭で説明できたら説明する自信はあります。
#include <iostream>
using namespace std;
#include<stdio.h>
#include<map>
#include<vector>
#include<math.h>
#include<set>
//プロジェクトオイラー問521
//a=10^6以下の素数
//b=aの範囲を素因数とする合成数
//c=10^6以上の素数
//d=(a xor b)*c
//A=(a∨b∨ｃ∨d)
//ヒント 10^6の答え37568404989
//計算時間 最新のパソコンで４５分。
const long long int LIMIT=pow(10,12);
const long long int MOD=pow(10,9);
const int UP=1000000;
int minPs[UP+1];
long long int minPSum[UP+1];
bool isPrime[UP+1];
std::vector<long long int> primes;
std::map<long long int,long long int> sumAB,countAB,sumAllAB;
std::map<long long int,long long int>::iterator it,it2,it3;
long long int calc2(long long int a){
long long int b;
if(a%2==1){
b=(a+1)/2;
}else{
b=a+1;
a/=2;
}
a=a % MOD;
b=b % MOD;
return (a*b)%MOD;
}
long long int calc(long long int a,long long int b){
//if(a==b)return a;
a=calc2(a);
b=calc2(b);
long long int c= (a-b)%MOD;
if(c<0)c+=MOD;
return c;
}
int minP(int n){
for(int i=2;i*i<=n;i++){
if(n%i==0)return i;
}
return n;
}
int main(){
for(int i=0;i<=UP;i++)isPrime[i]=true;
for(int i=2;i<=UP;i++){
if(isPrime[i]==false)continue;
for(int j=i*2;j<=UP;j+=i){
isPrime[j]=false;
}
}
for(int i=2;i<=UP;i++){
if(isPrime[i]){
primes.push_back(i);
}
}
std::cout<<LIMIT;
sumAB[LIMIT]=0;
sumAllAB[LIMIT]=1;
countAB[LIMIT]=1;
for(int i=primes.size()-1;i>=0;i--){
//割り算で考える
std::map<long long int,long long int> sumAB2,countAB2,sumAllAB2;
std::cout<<"第一ステップ残り"<<i<<"\n";
for(it=sumAB.begin();it!=sumAB.end();it++){
long long int p=primes[i];
long long int p2=p;
long long int f1,f2,s2,f3;
f1=(*it).first;
f2=f1;
s2=countAB[f2];
while(p2<=f1){
f3=f1/p2;
sumAB2[f3]=(sumAB2[f3]+s2*p)%MOD;
countAB2[f3]=(countAB2[f3]+s2)%MOD;
long long int temp=(sumAllAB[f1]*p2)%MOD;
sumAllAB2[f3]=(sumAllAB2[f3]+temp)%MOD;
p2*=p;
}
}
for(it=sumAB2.begin();it!=sumAB2.end();it++){
long long int f1,s1,f2,s2,f3,s3;
f1=(*it).first;
s1=(*it).second;
f3=f2=f1;
s2=countAB2[f2];
s3=sumAllAB2[f3];
sumAB[f1]=(sumAB[f1]+s1)%MOD;
countAB[f2]=(countAB[f2]+s2)%MOD;
sumAllAB[f1]=(sumAllAB[f1]+s3)%MOD;
}
}
long long int ansAB=0;
for(it=sumAB.begin();it!=sumAB.end();it++){
//std::cout<<"("<<(*it).first<<" "<<(*it).second<<")";
if((*it).first!=LIMIT){
ansAB=(ansAB+(*it).second)%MOD;
}
}
//ここで区切る
//Ｓ（c）を求める
long long int ansC=0;
long long int ansD=0;
for(it=sumAllAB.begin();(*it).first<=UP;it++){
long long int f=(*it).first;
long long int areaSum=calc(LIMIT/f,LIMIT/(f+1));
(*it).second=areaSum-(*it).second;
long long int areaSize=((LIMIT/f-LIMIT/(f+1)))%MOD;
countAB[f]=(areaSize-countAB[f])%MOD;
if(countAB[f]<0)countAB[f]+=MOD;
}
std::cout<<"\n\n";
it=sumAllAB.upper_bound(UP-1);
it--;
long long int sum=0;
minPSum[0]=0;
minPSum[1]=0;
for(int i=2;i<UP;i++){
minPs[i]=minP(i);
sum=(sum+minPs[i])%MOD;
minPSum[i]=sum;
}
long long int ansCount=0;
for(int i=UP-1;i>1;i--){
int f=i;
long long int s,m;
if(sumAllAB.find(f)==sumAllAB.end()){
s=LIMIT/f;//ここちょっと怪しい
m=1;
}else{
s=sumAllAB[f];
m=countAB[f];
}
std::cout<<f<<"\t"<<s<<"\n";
int p2=2;
std::set<int> memo;
while(p2*p2<=f){
int point=f/p2;
sumAllAB[point]=(sumAllAB[point]-s*p2)%MOD;
countAB[point]=(countAB[point]-m)%MOD;
if(sumAllAB[point]<0)sumAllAB[point]+=MOD;
if(countAB[point]<0)countAB[point]+=MOD;
ansD=(ansD+m*minPs[p2])%MOD;
if(ansD<0)ansD+=MOD;
p2++;
}
for(int p2=1;p2*p2<f;p2++){
int pR=f/(p2+1);
int pL=f/p2;
long long int w=pL-pR;
long long int dell=calc(pL,pR);
if((p2+1)*(p2+1)>=f){
pR=(int)sqrt(f);
w=pL-pR;
dell=calc(pL,pR);
}
sumAllAB[p2]=(sumAllAB[p2]-dell*s)%MOD;
countAB[p2]=(countAB[p2]-w*m)%MOD;
if(sumAllAB[p2]<0)sumAllAB[p2]+=MOD;
if(countAB[p2]<0)countAB[p2]+=MOD;
long long int a=minPSum[pL]-minPSum[pR];
ansD=(ansD+a*m)%MOD;
if(ansD<0)ansD+=MOD;
}
ansCount+=countAB[f];
ansC=(ansC+s)%MOD;
}
std::cout<<1<<"\t"<<sumAllAB[1];
ansC=ansC+sumAllAB[1];
std::cout<<"\nS(a∨b)="<<ansAB;
std::cout<<"\nS(c)="<<ansC;
std::cout<<"\nS(d)="<<ansD;
std::cout<<"\n ans=S(a∨b∨c∨d)="<<(ansAB+ansC+ansD)%MOD;
std::cout<<"\ncount="<<ansCount;
}
最終更新：2015年07月19日 01:34
