高精度超几何分布
#include stdio.h #include stdlib.h #include string.h #include stdbool.h #define MAX_DIGITS 2000 // 高精度数最大位数 #define MAX_PRIMES 200 // 最大质因数数量 // 高精度整数结构体低位在前存储 typedef struct { int digits[MAX_DIGITS]; int len; // 有效位数 } BigInt; // 质因数结构体 typedef struct { int prime; int exp; } Factor; // 质因数列表 typedef struct { Factor factors[MAX_PRIMES]; int len; // 质因数个数 } FactorList; // 高精度整数操作 // 初始化高精度整数为1 void bigIntInit(BigInt *num) { memset(num-digits, 0, sizeof(num-digits)); num-len 1; num-digits[0] 1; } // 高精度整数乘以普通整数 void bigIntMultiply(BigInt *num, int factor) { int carry 0; for (int i 0; i num-len; i) { long long product (long long)num-digits[i] * factor carry; num-digits[i] product % 10; carry product / 10; } // 处理剩余进位 while (carry 0 num-len MAX_DIGITS) { num-digits[num-len] carry % 10; carry / 10; num-len; } } // 将高精度整数转换为字符串高位在前 void bigIntToString(BigInt *num, char *str) { int idx 0; if (num-len 1 num-digits[0] 0) { str[idx] 0; str[idx] \0; return; } // 从高位到低位拼接 for (int i num-len - 1; i 0; i--) { str[idx] num-digits[i] 0; } str[idx] \0; } // 质因数相关操作 // Legendre公式计算n!中质因数p的指数 int legendre(int n, int p) { int count 0; while (n 0) { n / p; count n; } return count; } // 分解n的质因数 void primeFactors(int n, FactorList *factors) { factors-len 0; if (n 2) return; // 分解2 if (n % 2 0) { factors-factors[factors-len].prime 2; factors-factors[factors-len].exp 0; while (n % 2 0) { factors-factors[factors-len].exp; n / 2; } factors-len; } // 分解奇数 for (int i 3; i * i n; i 2) { if (n % i 0) { factors-factors[factors-len].prime i; factors-factors[factors-len].exp 0; while (n % i 0) { factors-factors[factors-len].exp; n / i; } factors-len; } } // 剩余的大质数 if (n 1) { factors-factors[factors-len].prime n; factors-factors[factors-len].exp 1; factors-len; } } // 获取组合数C(n,k)的质因数分解 void combPrimeFactors(int n, int k, FactorList *factors) { factors-len 0; if (k 0 || k n) return; if (k 0 || k n) return; // 收集n以内所有质因数 FactorList allPrimes; allPrimes.len 0; for (int i 2; i n; i) { FactorList temp; primeFactors(i, temp); // 去重添加质因数 for (int j 0; j temp.len; j) { bool exists false; for (int m 0; m allPrimes.len; m) { if (allPrimes.factors[m].prime temp.factors[j].prime) { exists true; break; } } if (!exists) { allPrimes.factors[allPrimes.len].prime temp.factors[j].prime; allPrimes.len; } } } // 计算每个质因数的指数 for (int i 0; i allPrimes.len; i) { int p allPrimes.factors[i].prime; int exp legendre(n, p) - legendre(k, p) - legendre(n - k, p); if (exp 0) { factors-factors[factors-len].prime p; factors-factors[factors-len].exp exp; factors-len; } } } // 合并两个质因数列表指数相加 void mergeFactorLists(FactorList *dst, FactorList *a, FactorList *b) { // 复制a到dst memcpy(dst, a, sizeof(FactorList)); // 合并b到dst for (int i 0; i b-len; i) { int p b-factors[i].prime; int exp b-factors[i].exp; bool found false; for (int j 0; j dst-len; j) { if (dst-factors[j].prime p) { dst-factors[j].exp exp; found true; break; } } if (!found) { dst-factors[dst-len].prime p; dst-factors[dst-len].exp exp; dst-len; } } } // 约分抵消分子分母的公共质因数 void reduceFactors(FactorList *numFactors, FactorList *denFactors) { for (int i 0; i numFactors-len; i) { int p numFactors-factors[i].prime; for (int j 0; j denFactors-len; j) { if (denFactors-factors[j].prime p) { int minExp numFactors-factors[i].exp denFactors-factors[j].exp ? numFactors-factors[i].exp : denFactors-factors[j].exp; numFactors-factors[i].exp - minExp; denFactors-factors[j].exp - minExp; // 移除指数为0的质因数 if (numFactors-factors[i].exp 0) { // 移位删除 for (int m i; m numFactors-len - 1; m) { numFactors-factors[m] numFactors-factors[m 1]; } numFactors-len--; i--; // 回退索引 } if (denFactors-factors[j].exp 0) { for (int m j; m denFactors-len - 1; m) { denFactors-factors[m] denFactors-factors[m 1]; } denFactors-len--; j--; } break; } } } } // 超几何分布计算 void hypergeometricDistribution(int N, int K, int n, int k, BigInt *numerator, BigInt *denominator) { // 边界检查 if (k 0 || k K || (n - k) 0 || (n - k) (N - K)) { bigIntInit(numerator); numerator-digits[0] 0; bigIntInit(denominator); denominator-digits[0] 1; return; } // 1. 计算分子质因数C(K,k) * C(N-K, n-k) FactorList c1, c2, numFactors; combPrimeFactors(K, k, c1); combPrimeFactors(N - K, n - k, c2); mergeFactorLists(numFactors, c1, c2); // 2. 计算分母质因数C(N,n) FactorList denFactors; combPrimeFactors(N, n, denFactors); // 3. 约分 reduceFactors(numFactors, denFactors); // 4. 高精度计算分子 bigIntInit(numerator); for (int i 0; i numFactors.len; i) { int p numFactors.factors[i].prime; int exp numFactors.factors[i].exp; for (int j 0; j exp; j) { bigIntMultiply(numerator, p); } } // 5. 高精度计算分母 bigIntInit(denominator); for (int i 0; i denFactors.len; i) { int p denFactors.factors[i].prime; int exp denFactors.factors[i].exp; for (int j 0; j exp; j) { bigIntMultiply(denominator, p); } } } // 主函数测试 int main() { // 测试用例1小参数验证 int N1 10, K1 5, n1 3, k1 2; BigInt num1, den1; char numStr1[MAX_DIGITS], denStr1[MAX_DIGITS]; hypergeometricDistribution(N1, K1, n1, k1, num1, den1); bigIntToString(num1, numStr1); bigIntToString(den1, denStr1); printf(测试用例1 (N%d, K%d, n%d, k%d):\n, N1, K1, n1, k1); printf( 分子: %s\n, numStr1); printf( 分母: %s\n, denStr1); printf( 最简分数: %s/%s\n\n, numStr1, denStr1); // 测试用例2大参数验证 int N2 1000, K2 358, n2 600, k2 240; BigInt num2, den2; char numStr2[MAX_DIGITS], denStr2[MAX_DIGITS]; printf(正在计算大参数请稍候...\n); hypergeometricDistribution(N2, K2, n2, k2, num2, den2); bigIntToString(num2, numStr2); bigIntToString(den2, denStr2); printf(测试用例2 (N%d, K%d, n%d, k%d):\n, N2, K2, n2, k2); printf( 分子: %s\n, numStr2); printf( 分母: %s\n, denStr2); printf( 最简分数: %s/%s\n, numStr2, denStr2); return 0; }
