POJ 1091（数论）

题目大意是给定两个整数n和m，求出长度为n+1满足条件的数列data的个数，数列的要求下：
1）1<=data[i]<=m,for1<=i<=n
2)data[n+1]=m;
3）这个n+1个数满足：存在x1,x2,...,xn,xn+1,满足x1*data[1]+x2*data[2]+...+x(n+1)*data[n+1]=1；
根据数论的知识，若存在这样的x1,x2...xn+1,则data[1],data[2]...data[n+1]的最大公约数为1

证明：若data[1],data[2]...data[n+1]满足题意，并且存在最大公约数d（为整数）;则x1*data[1]+x2*data[2]+...+x(n+1)*data[n+1]的和是d的整数倍，必不等于1

其实举个例子就明白了，例如:n=2,m=360
360=3^2*2^3*5  所有不满足条件的数列，最大公约数是360质因子的乘积，只要将这些组合去掉，就是要求的答案

具体解题步骤如下：
1、求出满m的所有质因子,存入数组num
2、求出总的序列个数吗m^n
3、设t(k)表示数列最大公约数为（k个质因子乘积）的数列的个数

f=m^n-t(1)+t(2)-t(3)+..(-1)^k*t(k);
答案 = (m ^ n) - （有公因数2的n元组）- （有公因数3的n元组）- （有公因数5的n元组）+ （有公因数2，3的n元组） +（有公因数2，5的n元组） + （有公因数3，5的n元组）- （有公因数2，3，5的n元组）。这个比公式形象些
有公因数d的n元组，每个位置上有 (m/d)个选择（1 ~ m里面有m/d个d的倍数），根据乘法原理，可以得出有公因数d的n元组有 (m/d)^n 个。

  1 #include<iostream>
 2 #include<cstdio>
 3 #include<cmath>
 4 using namespace std;
 5 __int64 n,m,per,total;
 6 __int64 s[130000],num[130000];
 7 
 8 void totalnum(__int64 x)//求质因子
 9 {
10     __int64 i;
11     total=0;
12     for(i=2;i*i<=x;i++)
13     {
14         if(x%i==0)
15         {
16             while(x%i==0) x=x/i;
17             num[total++]=i;
18         }
19     }
20     if(x!=1) num[total++]=x;
21 }
22 
23 __int64 por(__int64 x,__int64 y)//总的序列个数
24 {
25     __int64 i,k;
26     k=x;
27     for(i=1;i<y;i++)
28         x=k*x;
29     return x;
30 }
31 
32 void get(__int64 a,__int64 b,__int64 c)
33 {//a:序列起始位置 b个质因子乘积，c:公共质因子个数
34     __int64 i;
35     if(b==c)
36     {
37         __int64 t=m;
38         for(i=0;i<c;i++)
39             t=t/s[i];
40         per+=por(t,n);
41     }
42     else 
43     {
44         for(i=a;i<total;i++)
45         {
46             s[b]=num[i];
47             get(i+1,b+1,c);
48         }
49     }
50 }
51 
52 int main()
53 {
54     __int64 sum,i;
55     while(scanf("%I64d%I64d",&n,&m)!=EOF)
56     {
57         totalnum(m);
58         sum=por(m,n);
59         for(i=1;i<=total;i++)
60         {
61             per=0;
62             get(0,0,i);
63             if(i%2==1)
64                 sum-=per;
65             else sum+=per;
66         }
67         printf("%I64d\n",sum);
68     }
69     return 0;
70 }