數論

矩陣快速冪

斐波那契第n項和

#include <bits/stdc++.h>
using namespace std;
#define ne 2
#define INF 1000000009

struct node
{
    long long a[ne][ne];
};

node mult(node a, node b)
{
    node c = {0};
    for(long long i = 0; i < ne; i++)
    for(long long j = 0; j < ne; j++)
    for(long long k = 0; k < ne; k++)
    {
        c.a[i][j] += (a.a[i][k] * b.a[k][j]) % INF;
        c.a[i][j] %= INF;
    }
    return
 
 c;
}
node slove(node x,long long n)
{
    node a=x;
    while(n)
    {
        if(n&1) a=mult(a,x);
        x=mult(x,x);
        n>>=1;
    }
    return a;
}
node input(long long n)
{
    node x={0};
    x.a[0][0]=1;
    x.a[0][1]=1;
    x.a[1][0]=1;
    x.a[1][1]=0;
    x=slove(x,n);
    return
 
 x;
}

int main()
{
    long long n;
    while(cin>>n)
    {
        node x=input(n);
        cout<<x.a[1][1]<<endl;
    }
}

乘法逆元

#include <bits/stdc++.h>
using namespace std;
int exgcd(int a,int b,int &x,int &y) 
{
    if (b==0)
    {
        x=1,y=0;
        return
 
 a;
    }
    int q=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return q;
}

int main()
{
    int a,b,x,y;
    cin>>a>>b;
    exgcd(a,b,x,y);
    if(x<0) x+=b;
    cout<<x<<endl;
}

凸包

安德魯演算法 
安德魯演算法判斷凸包需要用到叉積 
關於叉積


首先在二維座標下介紹一些定義：
點：A（x1，y1），B（x2，y2）

向量：向量AB=（ x2 - x1 , y2 - y1 ）= ( x ,  y );

向量的模 |AB| = sqrt ( x*x+y*y );



向量的點積： 結果為 x1*x2 + y1*y2。

點積的結果是一個數值。

點積的集合意義：我們以向量 a 向向量 b 做垂線，則 | a | * cos（a,b）為 a 在向量 b 上的投影，即點積是一個向量在另一個向量上的投影乘以另一個向量。且滿足交換律

應用：可以根據集合意義求兩向量的夾角，
cos（a,b） =（ 向量a * 向量b ） / （| a | * | b |） =  x1*x2 + y1*y2 / （| a | * | b |）



向量的叉積： 結果為 x1*y2-x2*y1

叉積的結果也是一個向量，是垂直於向量a，b所形成的平面，如果看成三維座標的話是在 z 軸上，上面結果是它的模。

方向判定：右手定則，（右手半握，大拇指垂直向上，四指右向量a握向b，大拇指的方向就是叉積的方向）

叉積的集合意義：1：其結果是a和b為相鄰邊形成平行四邊形的面積。

2：結果有正有負，有sin（a，b）可知和其夾角有關，夾角大於180°為負值。

3：叉積不滿足交換律

應用：

1：通過結果的正負判斷兩向量之間的順逆時針關係

若 a x b > 0表示a在b的順時針方向上

若 a x b < 0表示a在b的逆時針方向上

若 a x b == 0表示a在b共線，但不確定方向是否相同

在順時針方向上為凸包 
否則是凹進去的。 
以此依次加點。 
將座標按照x從小到大排序。 
先加上面的凸包 然後加下面的凸包。

凸包
#include <bits/stdc++.h>
using namespace std;
#define mo 100005
struct node
{
    int x,y;
}d[mo],p[mo];

int chaji(node a,node b,node c)//叉積計算
{
    return ((b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x));
}


bool cmp(node a,node b)//座標點排序
{
    if(a.y!=b.y)return a.y<b.y;
    return a.x<b.x;
}
void run(int n)
{
    sort(d,d+n,cmp);
    int tot=0;
    for(int i=0;i<n;i++)
    {
        while(tot>1&&chaji(p[tot-2],p[tot-1],d[i])<0) tot--;
        p[tot++]=d[i];
    }

    for(int i=n-2;i>=0;i--)
    {
        while(tot>1&&chaji(p[tot-2],p[tot-1],d[i])<0) tot--;
        p[tot++]=d[i];
    }
    cout<<tot-1<<endl;
    for(int i=0;i<tot-1;i++)
    {
        cout<<p[i].x<<' '<<p[i].y<<endl;
    }

}

int main()
{
    int n;
    while(cin>>n)
    {
        for(int i=0;i<n;i++)
            cin>>d[i].x>>d[i].y;

        run(n);
    }



}

計算幾何

多邊形面積交

#include<cstdio>  
#include<iostream>  
#include<algorithm>  
#include<cstring>  
#include<cmath>  
using namespace std;  
#define maxn 510  
const double eps=1E-8;  
int sig(double d){  
    return(d>eps)-(d<-eps);  
}  
struct Point{  
    double x,y; Point(){}  
    Point(double x,double y):x(x),y(y){}  
    bool operator==(const Point&p)const{  
        return sig(x-p.x)==0&&sig(y-p.y)==0;  
    }  
};  
double cross(Point o,Point a,Point b){  
    return(a.x-o.x)*(b.y-o.y)-(b.x-o.x)*(a.y-o.y);  
}  
double area(Point* ps,int n){  
    ps[n]=ps[0];  
    double res=0;  
    for(int i=0;i<n;i++){  
        res+=ps[i].x*ps[i+1].y-ps[i].y*ps[i+1].x;  
    }  
    return res/2.0;  
}  
int lineCross(Point a,Point b,Point c,Point d,Point&p){  
    double s1,s2;  
    s1=cross(a,b,c);  
    s2=cross(a,b,d);  
    if(sig(s1)==0&&sig(s2)==0) return 2;  
    if(sig(s2-s1)==0) return 0;  
    p.x=(c.x*s2-d.x*s1)/(s2-s1);  
    p.y=(c.y*s2-d.y*s1)/(s2-s1);  
    return 1;  
}  
//多邊形切割  
//用直線ab切割多邊形p，切割後的在向量(a,b)的左側，並原地儲存切割結果  
//如果退化為一個點，也會返回去,此時n為1  
void polygon_cut(Point*p,int&n,Point a,Point b){  
    static Point pp[maxn];  
    int m=0;p[n]=p[0];  
    for(int i=0;i<n;i++){  
        if(sig(cross(a,b,p[i]))>0) pp[m++]=p[i];  
        if(sig(cross(a,b,p[i]))!=sig(cross(a,b,p[i+1])))  
            lineCross(a,b,p[i],p[i+1],pp[m++]);  
    }  
    n=0;  
    for(int i=0;i<m;i++)  
        if(!i||!(pp[i]==pp[i-1]))  
            p[n++]=pp[i];  
    while(n>1&&p[n-1]==p[0])n--;  
}  
//---------------華麗的分隔線-----------------//  
//返回三角形oab和三角形ocd的有向交面積,o是原點//  
double intersectArea(Point a,Point b,Point c,Point d){  
    Point o(0,0);  
    int s1=sig(cross(o,a,b));  
    int s2=sig(cross(o,c,d));  
    if(s1==0||s2==0)return 0.0;//退化，面積為0  
    if(s1==-1) swap(a,b);  
    if(s2==-1) swap(c,d);  
    Point p[10]={o,a,b};  
    int n=3;  
    polygon_cut(p,n,o,c);  
    polygon_cut(p,n,c,d);  
    polygon_cut(p,n,d,o);  
    double res=fabs(area(p,n));  
    if(s1*s2==-1) res=-res;return res;  
}  
//求兩多邊形的交面積  
double intersectArea(Point*ps1,int n1,Point*ps2,int n2){  
    if(area(ps1,n1)<0) reverse(ps1,ps1+n1);  
    if(area(ps2,n2)<0) reverse(ps2,ps2+n2);  
    ps1[n1]=ps1[0];  
    ps2[n2]=ps2[0];  
    double res=0;  
    for(int i=0;i<n1;i++){  
        for(int j=0;j<n2;j++){  
            res+=intersectArea(ps1[i],ps1[i+1],ps2[j],ps2[j+1]);  
        }  
    }  
    return res;//assumeresispositive!  
}  
//hdu-3060求兩個任意簡單多邊形的並面積  
Point ps1[maxn],ps2[maxn];  
int n1,n2;  
int main(){  
    while(scanf("%d%d",&n1,&n2)!=EOF){  
        for(int i=0;i<n1;i++)  
            scanf("%lf%lf",&ps1[i].x,&ps1[i].y);  
        for(int i=0;i<n2;i++)  
            scanf("%lf%lf",&ps2[i].x,&ps2[i].y);  
        double ans=intersectArea(ps1,n1,ps2,n2);  
        ans=fabs(area(ps1,n1))+fabs(area(ps2,n2))-ans;//容斥  
        printf("%.2f\n",ans);  
    }  
    return 0;  
}  

多邊形面積交

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<queue>
#include<map>
#include<stack>
#include<set>

using namespace std;

const int maxn=555;
const int maxisn=10;
const double eps=1e-8;
const double pi=acos(-1.0);

int dcmp(double x){
    if(x>eps) return 1;
    return x<-eps ? -1 : 0;
}
inline double Sqr(double x){
    return x*x;
}
struct Point{
    double x,y;
    Point(){x=y=0;}
    Point(double x,double y):x(x),y(y){};
    friend Point operator + (const Point &a,const Point &b) {
        return Point(a.x+b.x,a.y+b.y);
    }
    friend Point operator - (const Point &a,const Point &b) {
        return Point(a.x-b.x,a.y-b.y);
    }
    friend bool operator == (const Point &a,const Point &b) {
        return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
    }
    friend Point operator * (const Point &a,const double &b) {
        return Point(a.x*b,a.y*b);
    }
    friend Point operator * (const double &a,const Point &b) {
        return Point(a*b.x,a*b.y);
    }
    friend Point operator / (const Point &a,const double &b) {
        return Point(a.x/b,a.y/b);
    }
    friend bool operator < (const Point &a, const Point &b) {
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    }
    inline double dot(const Point &b)const{
        return x*b.x+y*b.y;
    }
    inline double cross(const Point &b,const Point &c)const{
        return (b.x-x)*(c.y-y)-(c.x-x)*(b.y-y);
    }

};

Point LineCross(const Point &a,const Point &b,const Point &c,const Point &d){
    double u=a.cross(b,c),v=b.cross(a,d);
    return Point((c.x*v+d.x*u)/(u+v),(c.y*v+d.y*u)/(u+v));
}
double PolygonArea(Point p[],int n){
     if(n<3) return 0.0;
     double s=p[0].y*(p[n-1].x-p[1].x);
     p[n]=p[0];
     for(int i=1;i<n;i++){
        s+=p[i].y*(p[i-1].x-p[i+1].x);
     }
     return fabs(s*0.5);
}
double CPIA(Point a[],Point b[],int na,int nb){
    Point p[maxisn],temp[maxisn];
    int i,j,tn,sflag,eflag;
    a[na]=a[0],b[nb]=b[0];
    memcpy(p,b,sizeof(Point)*(nb+1));
    for(i=0;i<na&&nb>2;++i){
        sflag=dcmp(a[i].cross(a[i+1],p[0]));
        for(j=tn=0;j<nb;++j,sflag=eflag){
            if(sflag>=0) temp[tn++]=p[j];
            eflag=dcmp(a[i].cross(a[i+1],p[j+1]));
            if((sflag^eflag)==-2)
                temp[tn++]=LineCross(a[i],a[i+1],p[j],p[j+1]);
        }
        memcpy(p,temp,sizeof(Point)*tn);
        nb=tn,p[nb]=p[0];
    }
    if(nb<3) return 0.0;
    return PolygonArea(p,nb);
}
double SPIA(Point a[],Point b[],int na,int nb){
    int i,j;
    Point t1[4],t2[4];
    double res=0.0,if_clock_t1,if_clock_t2;
    a[na]=t1[0]=a[0];
    b[nb]=t2[0]=b[0];
    for(i=2;i<na;i++){
        t1[1]=a[i-1],t1[2]=a[i];
        if_clock_t1=dcmp(t1[0].cross(t1[1],t1[2]));
        if(if_clock_t1<0) swap(t1[1],t1[2]);
        for(j=2;j<nb;j++){
            t2[1]=b[j-1],t2[2]=b[j];
            if_clock_t2=dcmp(t2[0].cross(t2[1],t2[2]));
            if(if_clock_t2<0) swap(t2[1],t2[2]);
            res+=CPIA(t1,t2,3,3)*if_clock_t1*if_clock_t2;
        }
    }
    return res;//面積交
    //return PolygonArea(a,na)+PolygonArea(b,nb)-res;//面積並
}

Point a[222],b[222];
Point aa[222],bb[222];

int main(){


    double x1,y1,x2,y2;
    double x3,y3,x4,y4;
    while(scanf("%lf %lf %lf %lf",&x1,&y1,&x2,&y2)!=EOF){
        scanf("%lf %lf %lf %lf",&x3,&y3,&x4,&y4);
        a[0]=Point(x1,y1);   // 逆時針？？
        a[1]=Point(x2,y1);
        a[2]=Point(x1,y2);

        b[0]=Point(x3,y3);
        b[1]=Point(x4,y3);
        b[2]=Point(x4,y4);
        b[3]=Point(x3,y4);

        printf("%.8f\n",fabs(SPIA(a,b,3,4)));
        //printf("%.8f\n",ConvexPolygonArea(out,m));
    }
    return 0;
}

圓與多邊形相交

#include <iostream>
#include <cstdio>
#include <cmath>
#include <vector>
#include <cstring>
#include <algorithm>
#include <string>
#include <set>
#include <functional>
#include <numeric>
#include <sstream>
#include <stack>

#define CL(arr, val) memset(arr, val, sizeof(arr))
#define REP(i, n)for((i) = 0; (i) < (n); ++(i))
#define FOR(i, l, h)for((i) = (l); (i) <= (h); ++(i))
#define FORD(i, h, l)for((i) = (h); (i) >= (l); --(i))
#define L(x)(x) << 1
#define R(x)(x) << 1 | 1
#define MID(l, r)(l + r) >> 1
#define Min(x, y)x < y ? x : y
#define Max(x, y)x < y ? y : x
#define E(x)(1 << (x))
#define iabs(x) (x) < 0 ? -(x) : (x)
#define OUT(x)printf("%I64d\n", x)
#define lowbit(x)(x)&(-x)
#define Read()freopen("data.in", "r", stdin)
#define Write()freopen("data.out", "w", stdout);

const double eps = 1e-10;
typedef long long LL;
const int inf = ~0u>>2;

using namespace std;

inline double max (double a, double b) { if (a > b) return a; else return b; }
inline double min (double a, double b) { if (a < b) return a; else return b; }
inline int fi (double a)
{
    if (a > eps) return 1;
    else if (a >= -eps) return 0;
    else return -1;
}
class Vector
{
public:
    double x, y;
    Vector (void) {}
    Vector (double x0, double y0) : x(x0), y(y0) {}
    double operator * (const Vector& a) const { return x * a.y - y * a.x; }
    double operator % (const Vector& a) const { return x * a.x + y * a.y; }
    Vector verti (void) const { return Vector(-y, x); }
    double length (void) const { return sqrt(x * x + y * y); }
    Vector adjust (double len)
    {
        double ol = len / length();
        return Vector(x * ol, y * ol);
    }
    Vector oppose (void) { return Vector(-x, -y); }
};
class point
{
public:
    double x, y;
    point (void) {}
    point (double x0, double y0) : x(x0), y(y0) {}
    Vector operator - (const point& a) const { return Vector(x - a.x, y - a.y); }
    point operator + (const Vector& a) const { return point(x + a.x, y + a.y); }
};
class segment
{
public:
    point a, b;
    segment (void) {}
    segment (point a0, point b0) : a(a0), b(b0) {}
    point intersect (const segment& s) const
    {
        Vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b;
        double s1 = v1 * v2, s2 = v3 * v4;
        double se = s1 + s2;
        s1 /= se, s2 /= se;
        return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1);
    }
    point pverti (const point& p) const
    {
        Vector t = (b - a).verti();
        segment uv(p, p + t);
        return intersect(uv);
    }
    bool on_segment (const point& p) const
    {
        if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 &&
            fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true;
        else return false;
    }
};

double radius;
point polygon[110];

double kuras_area (point a, point b, double cx, double cy)
{
    point ori(cx, cy);
    int sgn = fi((b - ori) * (a - ori));
    double da = (a - ori).length(), db = (b - ori).length();
    int ra = fi(da - radius), rb = fi(db - radius);
    double angle = acos(((b - ori) % (a - ori)) / (da * db));
    segment t(a, b); point h, u; Vector seg;
    double ans, dlt, mov, tangle;

    if (fi(da) == 0 || fi(db) == 0) return 0;
    else if (sgn == 0) return 0;
    else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn;
    else if (ra >= 0 && rb >= 0)
    {
        h = t.pverti(ori);
        dlt = (h - ori).length();
        if (!t.on_segment(h) || fi(dlt - radius) >= 0)
            return radius * radius * (angle / 2) * sgn;
        else
        {
            ans = radius * radius * (angle / 2);
            tangle = acos(dlt / radius);
            ans -= radius * radius * tangle;
            ans += radius * sin(tangle) * dlt;
            return ans * sgn;
        }
    }
    else
    {
        h = t.pverti(ori);
        dlt = (h - ori).length();
        seg = b - a;
        mov = sqrt(radius * radius - dlt * dlt);
        seg = seg.adjust(mov);
        if (t.on_segment(h + seg)) u = h + seg;
        else u = h + seg.oppose();
        if (ra == 1) swap(a, b);
        ans = fabs((a - ori) * (u - ori)) / 2;
        tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length()));
        ans += radius * radius * (tangle / 2);
        return ans * sgn;
    }
}

const double pi = acos(-1.0);

int main ()
{
    //freopen("data.in", "r", stdin);
    //radius 為圓的半徑
    int n;
    double area, x, y, cx, cy;//cx ，cy 為圓的座標
    double x0, y0, v0, th, t, g;
    double vx, vy;
    while(cin>>cx>>cy>>radius)
    {
        int n;
            cin>>n;
        for(int i = 0; i < n; ++i) {
            scanf("%lf%lf", &x, &y);
            polygon[i] = point(x, y);
        }
        area = 0;
        for (int i = 0; i < n; i++)
            area += kuras_area(polygon[i], polygon[(i + 1) % n], cx, cy);
        printf("%.2f\n", fabs(area));
    }
    return 0;
}