HDU6037 Expectation Division 期望、高维前缀和

传送门

---

设(f_x)表示答案，那么(f_x = frac{sumlimits_{d mid x} f_d}{sigma_0(x)} + 1 = frac{sigma_0(x) + sumlimits_{d mid x , d eq x} f_d}{sigma_0(x) - 1})

不难发现(f_x)只和(x)每个质因子出现次数构成的可重集合相关，所以考虑将所有质因子出现次数相同的放在一起考虑，经过搜索可以得到总的可重集的个数为(172513)种。

现在相当于设(f_S)表示质因子可重集合为(S)时的答案，那么仍然有(f_S = frac{sigma_0(x) + sumlimits_{T varsubsetneqq S} f_T}{sigma_0(x) - 1})，其中(x)表示任意一个满足质因子可重集合为(S)的数，当然(sigma_0(x))也只和(S)有关，这里这么写只是为了好写

我们需要快速求出(sumlimits_{T varsubsetneqq S} f_T)。不难发现这是一个高维前缀和的形式，但是并不能够先枚举维度、再枚举位置，因为计算是在线的，每一次计算必须要把之前的子集算出来，所以当前计算就一定要把当前的值算出来。我们考虑反过来，我们相当于已经枚举了一个位置，那么设(g_{S,j})表示对于(S)，对于所有满足：(S)和(T)中最大的(j-1)个元素相同，其余元素排序之后对应位置(S)中的元素大于等于(T)中元素的所有(T)的(f_T)的和。转移先枚举第一个满足(S)中元素大于(T)中元素的位置进行贡献，然后计算一下后缀和。那么(sumlimits_{T varsubsetneqq S} f_T = g_{S,1})。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;
#define int __int128

int read(){
	int a = 0; char c = getchar();
	while(!isdigit(c) && c != EOF) c = getchar();
	if(c == EOF) exit(0);
	while(isdigit(c)){
		a = a * 10 + c - 48; c = getchar();
	}
	return a;
}

const int prm[] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79};
int cnt; double ans[300003][20];
map < int , int > Hash;

int poww(int a , int b){
	int times = 1;
	while(b){
		if(b & 1) times = times * a;
		a = a * a; b >>= 1;
	}
	return times;
}

int arr[30];

int find(){
	int tms = 1;
	for(int i = 1 ; arr[i] ; ++i)
		tms = tms * poww(prm[i - 1] , arr[i]);
	return tms;
}

void dfs(int now , int id , int pre , int ys){
	if(now != 1){
		Hash[now] = ++cnt;
		int p = 1;
		for(int i = 1 ; i < id ; ++i)
			if(arr[i] > arr[i + 1]){
				--arr[i]; int id = Hash[find()];
				for(int j = p ; j <= i ; ++j) ans[cnt][j] = ans[id][j];
				p = i + 1; ++arr[i];
			}
		for(int i = id - 2 ; i > 0 ; --i)
			ans[cnt][i] = ans[cnt][i] + ans[cnt][i + 1];
		ans[cnt][0] = (ans[cnt][1] + ys) / (ys - 1);
	}
	for(int i = 1 ; i <= id ; ++i) ans[cnt][i] += ans[cnt][0];
	for(int j = 1 ; now * prm[id - 1] <= 1e24 && j <= pre ; ++j){
		arr[id] = j;
		dfs(now *= prm[id - 1] , id + 1 , j , ys * (j + 1));
	}
	arr[id] = 0;
}

bool cmp(int a , int b){return a > b;}

signed main(){
#ifndef ONLINE_JUDGE
	freopen("in","r",stdin);
	//freopen("out","w",stdout);
#endif
	dfs(1 , 1 , 1e9 , 1);
	int N , M , Case = 0;
	while(N = read()){
		M = read();
		for(int i = 1 ; i <= M ; ++i){
			int P = read() , cnt = 0;
			while(N % P == 0){N /= P; ++cnt;}
			arr[i] = cnt;
		}
		sort(arr + 1 , arr + M + 1 , cmp); arr[M + 1] = 0;
		printf("Case #%d: %.10lf
" , (signed)(++Case) , ans[Hash[find()]][0]);
	}
	return 0;
}