poj 3361 Gaussian Prime Factors 高斯素数约数

poj 3361 Gaussian Prime Factors

题意：

在复数 a+bj (a,b为整数)范围内，约数只有 1, -1, a+bj, -(a+bj)的称为高斯素数。求任给正整数N的所有素因子(|b|>a>0)。默认N<2e9

解法：

定理有云：
n≡3(mod 4)者，无因式分解。
n≡1(mod 4)者，可唯一分解为两个共轭高斯整数乘积。

当然做的时候并不知道！
我的思路是，既然可以分解因子，那么可以先分解成素因子，再把素因子分解出来。发现3,7等数是不能分解的，其他的素数(忽略2)都是奇数，要分解成乘积，最后不能有虚部存在，故得分解成共轭高斯整数相乘，即(a+bj)(a-bj)=n，即 a^2+b^2 = n,则a,b一奇一偶。且因为若a+bj不为素数，则可分解为 (c+dj)(e+fj)，那么(a+bj)(a-bj)=(c+dj)(c-dj)(e+fj)(e-fj)=n,此时前两项与后两项都为一个实整数，与n为素数矛盾，故只要分解出来，a+bj与a-bj则一定为n的素因子。
这样就好办了，对所有的N，都先进行素因子分解，然后对每个素因子再分解其高斯素因子。我选择遍历a,b来预处理出所有素数的标准因子a+bj（另一个则为a-bj）。

自己的坑：
低估了45000内素数数量，大约为十分之一个。
使用了set,但是貌似不是这么用的。(我用来验证m-i^2是否为完全平方数，多了一个log(num_pri)的复杂度)。应该是直接检验比较快。

so,学到了STL set:

myset.count(n)  ->得到数量
myset.insert(n) ->插入
myset.find(n)   ->得到迭代器，找不到则是末尾。 .erase(.find(n))

140+行代码，500+ms过。可能是预处理和数据结构比较费时间。看别人的当场暴力分解质因数，暴力分解素数的0ms过。只能说数据太水太水啦！~

我的代码：

#include<iostream>
#include<cstdio>
#include<string>
#include<cstring>
#include<vector>
#include<cmath>
#include<set>
#include<algorithm>
using namespace std;
typedef long long ll;
const int maxn = 45000;
const int maxq = 212;
struct stt{
    int a,b;
}q[maxn],tmpst;
vector <stt> ans;
set <int> myset;
int prim[5000];
bool isnotpri[maxn];
int prinum = 0;
int prime,sise;
int sq[maxq+1];//square

void ini(){
    memset(prim,0,sizeof prim);
    memset(isnotpri,false,sizeof isnotpri);
    memset(q,0,sizeof q);
    myset.clear();
    //memset(ans,0,sizeof ans);
    ans.clear();
    for(int i = 2; i <= maxq; i++){
        if(!isnotpri[i]){
            for(int j = i*i; j < maxn; j += i){
                isnotpri[j] = 1;
            }
        }
    }
    for(int i = 2; i <= maxn; i++){
        if(!isnotpri[i]) prim[prinum++]=i;
    }
    q[2].a=1;
    q[2].b=1;
    for(int i = 1; i <= maxq; i++){
        sq[i] = i*i;
    }
    int tmp;
    for(int i = 1; i <= maxq; i++){
        for(int j = i+1; j <= maxq; j += 2){
            tmp = sq[i]+sq[j];
            if(tmp >= maxn) break;
            if(!isnotpri[tmp]){
                q[tmp].a = i;
                q[tmp].b = j;
            }
        }
    }
    for(int i = 2; i <= maxn; i++){
        if(q[i].b == 0) q[i].a = i;
        myset.insert(i*i);
    }
}
void gener(int m){
    ans.clear();
    if(m == 0 || m == 1) return;
    for(int i = 0; i < prinum; i++){
        prime = prim[i];
        if(prime > m) break;
        if(m % prime == 0){
            ans.push_back(q[prime]);
            while(m % prime == 0){
                m /= prime;
            }
        }
    }
    if(m != 1){
        if(m%4==3){
            tmpst.a = m;
            tmpst.b = 0;
            ans.push_back(tmpst);
        }
        else{
            int tmp = (int)sqrt(m+0.1);
            int tmp2;
            for(int i = 1; i <= tmp; i++){
                if(myset.count(m-i*i)!=0){
                    tmpst.a = i;
                    tmpst.b = (int)sqrt(m-i*i+0.1);
                    ans.push_back(tmpst);
                    break;
                }
            }
        }
    }
}
bool cmp(stt i, stt j){
    if(i.a == j.a){
        return i.b < j.b;
    }
    return i.a < j.a;
}
void printout(){
    sise = ans.size();
    if(sise == 0){
        printf("\n");
        return;
    }
    sort(ans.begin(),ans.end(),cmp);
    int aa,bb;
    for(int i = 0; i < sise; i++){
        aa = ans[i].a, bb = ans[i].b;
        if(bb==0){
            printf(" %d",aa);
        }
        else{
            if(bb==1){
                printf(" %d+j, %d-j",aa,aa);
            }
            else{
                printf(" %d+%dj, %d-%dj",aa,bb,aa,bb);
            }

        }
        printf("%c",",\n"[i==sise-1]);
    }
}
int main(){
    freopen("in.txt","r",stdin);
    ini();
    int k = 1;
    int m;
    //return 0;
    while(scanf("%d",&m)!=EOF){
        printf("Case #%d:",k++);
        gener(m);
        printout();
    }
    return 0;
}

别人的代码：

from yjCola

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define MAXN 1005
struct ele
{
    int a,b;
    bool operator < (const ele &c) const
    {
        if(a==c.a)
        {
            if(abs(b)==abs(c.b))
                return b>c.b;
            else
                return abs(b)<abs(c.b);
        }
        else
            return a<c.a;
    }
}e[MAXN];
int up;
void get(int p)
{
    if((p-3)%4==0) e[up].a=p,e[up++].b=0;
    else
    {
        for(int i=1;i*i<=p;i++)
        {
            int j=sqrt(1.0*(p-i*i));
            if(i*i+j*j==p)
            {
                e[up].a=i,e[up++].b=j;
                e[up].a=i,e[up++].b=-j;
                break;
            }
        }
    }
}
int main()
{
    int cas=1,n;
    while(scanf("%d",&n)!=EOF)
    {
        up=0;
        for(int i=2;i*i<=n;i++)
            if(n%i==0)
            {
                get(i);
                while(n%i==0) n/=i;
            }
        if(n>1) get(n);
        sort(e,e+up);
        printf("Case #%d:",cas++);
        for(int i=0;i<up;i++)
        {
            printf(" %d",e[i].a);
            if(e[i].b<0)
            {
                if(e[i].b==-1) printf("-j");
                else printf("%dj",e[i].b);
            }
            else if(e[i].b>0)
            {
                if(e[i].b==1) printf("+j");
                else printf("+%dj",e[i].b);
            }
            if(i<up-1) printf(",");
        }
        printf("\n");
    }
    return 0;
}

疑问：

看到有个人说 1 = j·(-j)，为何1不是素因子，想想如果酱紫的话，j = 1·j, 所有的单位元都可以认为是素因子了。思考了下，可能是因为，由于a+bj和 j*(a+bj)是等价的因子，所以1 = j·(-j)又可以认为是1=1·(-1)，1与本身重合，所以 1 仍是作为单元，而非素因子。