本文主要探讨数学中的阶乘问题。比如求阶乘最后非零位的位数，计算阶乘的值，等等。
题目：http://acm.hdu.edu.cn/showproblem.php?pid=1042
题目：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1059
分析：题意比较简单，都是给定一个数，求阶乘的值。第二题的数比较大。
首先如果直接相乘肯定是不行的，稍微高级一点的做法就是采用万进位，这样就能过第一道题了。如果要暴
力，可以将阶乘写成幂的形式（先进行素数筛选）。用Java写了一个，能过这两个题，只不过效率太低。
代码：
import java.math.BigInteger;
import java.util.Arrays;
import java.util.Scanner;
public class Main {
public final static int N = 100005;
public static int[] p = new int[N];
public static boolean[] prime = new boolean[N];
public static int cnt;
@SuppressWarnings("resource")
public static void main(String[] args){
isprime();
Scanner cin = new Scanner(System.in);
int n = cin.nextInt();
BigInteger ans = BigInteger.ONE;
for(int i = 0; i < cnt; i++){
if(p[i] > n) break;
ans = ans.multiply(BigInteger.valueOf(p[i]).pow(getFactor(p[i], n)));
}
StringBuffer str = new StringBuffer("");
char[] ch = ans.toString().toCharArray();
for(int i = 0; i < ch.length; i++){
if(i % 1000 == 0 && i > 0){
System.out.println(str.toString());
str.delete(0, str.length());
}
str.append(ch[i]);
}
if(str.length() > 0){
System.out.println(str.toString());
}
}
public static void isprime(){
cnt = 0;
Arrays.fill(prime, true);
for(int i = 2; i < N; i++){
if(prime[i]){
p[cnt++] = i;
for(int j = i + i; j < N; j += i){
prime[j] = false;
}
}
}
}
public static int getFactor(int p, int n){
int ans = 0;
while(n != 0){
ans += n / p;
n /= p;
}
return ans;
}
}
这样的效率确实太低下，必须想更好的方法。有位大神的代码效率比较高，使用了FFT优化，我对其进行了修
改并加上了一些注释。如下
代码：
#include <iostream>
#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <math.h>
using namespace std;
//FFT模板开始
const double PI = acos(-1.0);
struct Virt
{
double r, i;
Virt(double r = 0.0, double i = 0.0)
{
this->r = r;
this->i = i;
}
Virt operator + (const Virt &x)
{
return Virt(r + x.r, i + x.i);
}
Virt operator - (const Virt &x)
{
return Virt(r - x.r, i - x.i);
}
Virt operator * (const Virt &x)
{
return Virt(r * x.r - i * x.i, i * x.r + r * x.i);
}
};
//雷德算法--倒位序
void Rader(Virt F[], int len)
{
int j = len >> 1;
for(int i = 1; i < len - 1; i++)
{
if(i < j) swap(F[i], F[j]);
int k = len >> 1;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}
//FFT实现
void FFT(Virt F[], int len, int on)
{
Rader(F, len);
for(int h = 2; h <= len; h <<= 1) //分治后计算长度为h的DFT
{
Virt wn( cos(-on * 2 * PI / h), sin(-on * 2 * PI / h)); //单位复根e^(2*PI/m)用欧拉公式展开
for(int j = 0; j < len; j += h)
{
Virt w(1, 0); //旋转因子
for(int k = j; k < j + h / 2; k++)
{
Virt u = F[k];
Virt t = w * F[k + h / 2];
F[k] = u + t; //蝴蝶合并操作
F[k + h / 2] = u - t;
w = w * wn; //更新旋转因子
}
}
}
if(on == -1)
for(int i = 0; i < len; i++)
F[i].r /= len;
}
//求卷积
void Conv(Virt a[], Virt b[], int n)
{
FFT(a, n, 1);
FFT(b, n, 1);
for(int i = 0; i < n; i++)
a[i] = a[i] * b[i];
FFT(a, n, -1);
}
//FFT模板结束
/***********************/
//代码正文
typedef long long LL;
const int N = 1000005;
const int maxn = 500000;
const int BITS = 1000;
const int LEN = 3;
Virt x1[N], x2[N];
int num[maxn];
int temp[maxn];
int tmp[20][maxn];
char buf[maxn];
int ans[maxn];
char rel[maxn * LEN];
LL bt[6]= {1, 10, 100, 1000, 10000, 100000};
void Solve(int x[], int &L1, int y[], int L2)
{
int L = 1;
while(L < 2 * L1 || L < 2 * L2) L <<= 1;
for(int i = 0; i < L; i++)
{
if(i < L1) x1[i].r = x[i];
else x1[i].r = 0;
if(i < L2) x2[i].r = y[i];
else x2[i].r = 0;
x1[i].i = x2[i].i = 0;
}
Conv(x1, x2, L);
LL over = 0;
L1 = 0;
for(int i = 0; i < L; i++)
{
over += x1[i].r + 0.5;
x[i] = over % BITS;
if(over) L1 = i;
over /= BITS;
}
if(over) x[++L1] = over;
++L1;
}
//将buf字符串转化为int数组
void CharToInt(char *buf, int x[], int &len)
{
int bLen = strlen(buf);
//求出余项的值
int val = 0;
for(int i = 0; i < bLen % LEN; i++)
val = val * 10 + buf[i] - '0';
len = 0;
if(val) x[len++] = val;
for(int i = bLen % LEN; i < bLen; i += LEN)
{
val = 0;
for (int j = 0; j < LEN; j++)
val = val * 10 + buf[i + j] - '0';
x[len++] = val;
}
//反转字符串
for(int i = 0; i < len / 2; i++)
swap(x[i], x[len - 1 - i]);
}
void PutStr(int x[], int n, char *buf, int &len, int BITS_LEN)
{
int j = 0;
sprintf(buf, "%d", x[n - 1]);
while(buf[j]) j++;
for(int i = n - 2; i >= 0; i--)
{
for(int k = 0; k < BITS_LEN; k++)
{
buf[j + k] = x[i] / bt[BITS_LEN - 1 - k] + '0';
x[i] %= bt[BITS_LEN - 1 - k];
}
j += BITS_LEN;
}
buf[j] = 0;
len = j;
}
//每1000位一千位地输出
void Print(char *buf, int len)
{
for(int i = 0; i < len; i += 1000)
{
char ch = buf[i + 1000];
buf[i + 1000] = 0;
printf("%s\n", buf + i);
buf[i + 1000] = ch;
}
}
void ConvertNum(int *num5, int len5, int *num, int &len)
{
PutStr(num5, len5, buf, len, 5);
CharToInt(buf, num, len);
}
void Run(int l, int r, int *num, int &m)
{
memset(temp, 0, (r - l + 1) * sizeof(int));
memset(num, 0, (r - l + 1) * sizeof(int));
temp[0] = 1;
m = 0;
for(int i = l; i <= r; i++)
{
LL top = 0;
for(int j = 0; j <= m; j++)
{
top += (LL)(temp[j]) * i;
temp[j] = top % 100000;
top /= 100000;
}
if(top)
temp[++m] = top;
}
ConvertNum(temp, m + 1, num, m);
}
//分治法
void Divide(int dept, int ans[], int &n, int l, int r)
{
if (r - l <= 200)
{
Run(l, r, ans, n);
return ;
}
int mid = (l + r) >> 1;
int rlen;
Divide(dept + 1, ans, n, l, mid);
Divide(dept + 1, tmp[dept] + mid * 2 + 2, rlen, mid + 1, r);
Solve(ans, n, tmp[dept] + mid * 2 + 2, rlen);
}
int main()
{
int n, m;
while (scanf("%d", &n) != EOF)
{
Divide(1, ans, m, 1, n);
PutStr(ans, m, rel, n, LEN);
Print(rel, n);
}
return 0;
}
