设\(f_x\)表示答案,那么\(f_x = \frac{\sum\limits_{d \mid x} f_d}{\sigma_0(x)} + 1 = \frac{\sigma_0(x) + \sum\limits_{d \mid x , d \neq x} f_d}{\sigma_0(x) - 1}\)
不难发现\(f_x\)只和\(x\)每个质因子出现次数构成的可重集合相关,所以考虑将所有质因子出现次数相同的放在一起考虑,经过搜索可以得到总的可重集的个数为\(172513\)种。
现在相当于设\(f_S\)表示质因子可重集合为\(S\)时的答案,那么仍然有\(f_S = \frac{\sigma_0(x) + \sum\limits_{T \varsubsetneqq S} f_T}{\sigma_0(x) - 1}\),其中\(x\)表示任意一个满足质因子可重集合为\(S\)的数,当然\(\sigma_0(x)\)也只和\(S\)有关,这里这么写只是为了好写
我们需要快速求出\(\sum\limits_{T \varsubsetneqq S} f_T\)。不难发现这是一个高维前缀和的形式,但是并不能够先枚举维度、再枚举位置,因为计算是在线的,每一次计算必须要把之前的子集算出来,所以当前计算就一定要把当前的值算出来。我们考虑反过来,我们相当于已经枚举了一个位置,那么设\(g_{S,j}\)表示对于\(S\),对于所有满足:\(S\)和\(T\)中最大的\(j-1\)个元素相同,其余元素排序之后对应位置\(S\)中的元素大于等于\(T\)中元素的所有\(T\)的\(f_T\)的和。转移先枚举第一个满足\(S\)中元素大于\(T\)中元素的位置进行贡献,然后计算一下后缀和。那么\(\sum\limits_{T \varsubsetneqq S} f_T = g_{S,1}\)。
#include//this code is written by Itstusing namespace std;#define int __int128int 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\n" , (signed)(++Case) , ans[Hash[find()]][0]); } return 0;}