高精度超几何分布

#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; }