POJ 1091 扩展欧拉函数

2010年7月20日 05:15

 欧拉函数F(n)=n \prod_{p|n}(1-\frac{1}{p})的概率解释:

1/p看做从1~n中某数有n的因子p的概率,则1-1/p是没有因子p的概率,最前面的n是指1~n的数的个数,由乘法公式可见F的意义

此题设函数F(n,m)为公因子为1的排列的个数,易见当m为质数时 F(n,m)=m^n-1

普遍看,总排列个数为m^n,取p|m,则1/p为1~m的某数有因子p的概率,1/{p^n}为前n个数都有p因子的概率,即这n+1个数有公因子p的概率

所以由乘法公式见:

F(n,m)=m^n\prod_{p|n}(1-\frac{1}{p^n})

恩,我理解的大概就是这个样子了...

评论(0) 阅读(2903)

POJ 2480 积性函数 欧拉函数

2010年7月19日 01:41

因为有性质

 gcd(i,m \times n)=gcd(i,m) \times gcd(i,n) 当m \perp n

所以gcd是积性函数

因为积性函数的和函数也是积性函数(《初等数论及其应用》P182),所以所求函数

F(N)=\sum_{i=1}^{N}gcd(i,N)是积性函数

问题简化到素数阶段,须求F(p^k),此时用到欧拉公式.因为有F(N)=\sum _{d|N}d\cdot \varphi (\frac{N}{d})

(易见因为\varphi (\frac{N}{d})代表的是从1到N中和N最大公约数gcd为d的数的个数)

所以F(p^k)=\sum _{i=0}^k p^i\cdot \varphi (p^{k-i}),又因为\varphi (p^k)=p^k-p^{k-1}则可求出

 

#include <iostream>
#include <cstdio>
using namespace std;
typedef unsigned long long u64;
#define maxn 50000
u64 primet[maxn];
int pnum = 0;
void getprimtable()//筛法求出素数
{
    bool flag[maxn] = {false};
    for (long i = 2; i < maxn; i++)
    {
        if (!flag[i])
        {
            primet[pnum++] = i;
            for (long j = i + i; j < maxn; j += i) flag[j] = true;
        }
    }
}

u64 eularpk(u64 p, int k)//欧拉求p^k
{
    if (k == 0) return 1;
    u64 ans = p - 1;
    while (--k)
    {
        ans *= p;
    }
    return ans;
}

u64 fpk(u64 p, int k)//求f(p^k)
{
    u64 ans = 0, t = 1;
    for (int i = 0; i <= k; i++)
    {
        ans += t * eularpk(p, k - i);
        t *= p;
    }
    return ans;
}

u64 f(u64 n)//求f 因式分解
{
    u64 sum = 1;
    int p, k;
    for (int i = 0; primet[i]*primet[i] <= n; i++)
    {
        p = primet[i];
        //cout<<p<<ends;
        if (n % p == 0)
        {
            k = 0;
            while (n % p == 0) {n /= p; k++;}
            sum *= fpk(p, k);
        }
    }
    if (n > 1) sum *= fpk(n, 1); // 分解到最后剩下一个素因子
    return sum;
}

int main()
{
    //freopen("in.txt", "r", stdin);
    getprimtable();
    u64 n;
    while (cin >> n)
    {
        cout << f(n) << endl;
    }
}

 

评论(0) 阅读(3792)

poj 2891 中国剩余定理

2010年7月16日 18:22

不互质的中国剩余定理,依次合并!

 

假设C  A1 (mod B1),C  A2 (mod B2)。令C = A1 + X1B,那么X1B1  A2  A1 (mod B2)。用扩展欧几里德算法求出X1,也就求出C。令B = lcm(B1B2),那么上面两条方程就可以被C  C (mod B)代替。迭代直到只剩下一条方程。

 

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;
__int64 x, y, t;
__int64 egcd(__int64 a, __int64 b)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    else
    {
        __int64 e = egcd(b, a % b);
        t = x; x = y; y = t - a / b * y;
        return e;
    }
}
int main()
{
    //freopen("in.txt", "r", stdin);
    __int64 m1, m2, r1, r2, d, c, t;
    bool sb;
    int n;
    while (cin >> n)
    {
        sb = 0;
        cin >> m1 >> r1;
        for (int i = 0; i < n - 1; i++)
        {
            cin >> m2 >> r2;
            if (sb) continue;
            d = egcd(m1, m2);
            c = r2 - r1;
            if (c % d)
            {
                sb = 1;
                continue;
            }
            t = m2 / d;
            x = (c / d * x % t + t) % t;

            r1 = m1 * x + r1;
            m1 = m1 * m2 / d;
        }
        if (sb)
            cout << -1 << endl;
        else
            cout << r1 << endl;
    }

}

 

评论(1) 阅读(6986)

hdoj 1028 1398 母函数

2010年7月11日 23:02

"把组合问题的加法法则和幂级数的t的乘幂的相加对应起来"

"母函数的思想很简单—就是把离散数列和幂级数一一对应起来,把离散数列间的相互结合关系对应成为幂级数间的运算关系,最后由幂级数形式来确定离散数列的构造. "

 

代码基本是这个样子的:

hdoj 1028 整数拆分

#include <iostream>
#include <cstdio>
using namespace std;
#define maxl 300
int ans[maxl+1], t[maxl+1];

int main()
{
    freopen("in.txt", "r", stdin);
    int n, i, j, k;
    while (cin >> n)
    {
        for (i = 0; i <= n; i++)
        {
            ans[i] = 1; t[i] = 0;
        }
        for (i = 2; i <= n; i++)
        {
            for (j = 0; j <= n; j++)
                for (k = 0; k + j <= n; k += i)
                {
                    t[j+k] += ans[j];
                }
            for (j = 0; j <= n; j++)
            {
                ans[j] = t[j]; t[j] = 0;
            }
        }
        cout << ans[n] << endl;
    }
}

评论(0) 阅读(2011)

Miller-Robbin 素数测试

2010年7月02日 18:13

如果一个数是伪素数,他几乎肯定是素数,如果一个数不是伪素数,他一定不是素数。

不断读取2~n-1的基b(s次),计算每次是否都有 exp(b,n-1)mod n == 1

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;
#define maxn 100001
#define s 10
bool isprim[maxn];

void shai()
{
    for(int i=0;i<maxn;i++) isprim[i]=1;
    int t;
    for(int i=2;i<=maxn/2;i++)
    {
        t=i+i;
        while(t<maxn)
        {
            isprim[t]=0;
            t+=i;
        }
    }
}

int modexp(int a,int b,int m)
{
    __int64 ans=1;__int64 t=a%m;
    while(b)
    {
        if(b&1)
        {
            ans=(ans*t)%m;
            b--;
        }
        t=(t*t)%m;
        b>>=1;
    }
    return ans;
}

bool MillerRobbin(int n)
{
    for(int i=1;i<=s;i++)
    {
        int a=rand()%(n-2)+2;
        if (modexp(a,n-1,n)!=1) return false;
    }
    return true;
}

int main()
{
    //freopen("in.txt", "r", stdin);
    shai();
    for(int i=3;i<maxn;i++)
    {
        if(MillerRobbin(i)!=isprim[i])
            cout<<i<<endl;
    }
}

原来的这个是抄黑书上的,貌似只是用费马小定理来检验的,下面的这个应该才是正确的吧,正确率确实很高:

bool MillerRobbin(int n)
{
    if(n==2) return 1;
    else if((n&1)==0) return 0;
    int a,d;
    long long t;
    for(int i=0;i<s;i++)
    {
        a=rand()%(n-1)+1;
        d=n-1;
        while((d&1)==0) d>>=1;
        t=modexp(a,d,n);
        while(d!=n-1&&t!=1&&t!=n-1)
        {
            t=(t*t)% n;
            d<<=1;
        }
        if(t==n-1||(d&1)==1) continue;
        else return 0;
    }
    return 1;
}

评论(4) 阅读(2485)