【C++】计算技巧_c++如何计算-优快云博客

本文链接：https://blog.youkuaiyun.com/qq_36393962/article/details/82731564

【C++】计算技巧

文章目录

【C++】计算技巧

取整

特别要注意，如果使用cmath库的函数会返回double等浮点类型而不是整型。所以这里涉及到浮点数强制转整型的问题

double m= -1.3; int a=(int)m;

注意，这样的强制类型转换是向0取整

(int)0.9999999999==0

注意，一个小于1的数，向0取整后，无论原本多么接近1，也得到0

int a=7;double m=(double)a;

注意，浮点数表示的整数大部分只是近似值，xx.999999或xx.00000002 这些情况是很常见的，而eps也主要是应对这种情况

四舍五入（round）：

事实上,#include <cmath>自带lround,llround，可以看源码。
不过还是应该掌握自己写int_round的方法。
```
#define EPS 1e-12
int int_round(double x)
{
    return (x>=0.0) ? (int)(x+0.5+EPS):(int)(x-0.5-EPS);
}

//或使用函数
int int_round(double x)
{
    double tmp = round(x); //有一些老的编译器不支持
    double tmp =(x >= 0.0) ? floor(x + 0.5 + EPS) : ceil(x - 0.5 - EPS);
    return x>=0 ? (int)(tmp + EPS) : (int)(tmp - EPS);
}
//输出四舍五入
//C stdio
printf("%.0lf\n",double_value);

//C++ iostream,iomanip
cout.setf(ios::fixed);//定点数输出 
cout<<setprecision(0)<<double_value<<endl; //iomanip头文件
```
向上(ceil)/下(floor)取整：和int_round的差不多。
计算模

mod>0

如果要使计算出的结果为正数，则应事先加个K*mod
```
//注意先加mod再%mod，这样对负数也能正常处理
(a + mod) % mod;
```
上面这个只是很简单的情况，还有K>=2的（a<(-mod)），不过如果出现这种情况，一般也是没处理好。不过有更好的办法处理
```
(a % mod + mod) % mod
```
这样就完美了，只要小心溢出就行。

计算占用块数

//block_size为整数
//C++的除法会自动取整数
//最后一块没占满也算一块
(a+block_size-1)/block_size

一行gcd和lcm

//不考虑细节时，这个已经足够了
int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
int lcm(int a,int b) {return a/gcd(a,b)*b;}

unsigned截去高位相当于变相取模

unsigned int 加减乘计算时，如果不对mod取膜，会自动对(UINT_MAX+1)取膜，即截去高位。这是，要注意，只有加减乘可以保证取模的正确性。而且其间也不能模其他值，除非(UINT_MAX+1) % mod == 0

a^b%mod 和 a*b%mod

LL mul_mod(LL a,LL b,LL mod)
{
    LL ans=0;
    for(;b;b>>=1)
    {
        if(1&b)ans=(ans+a)%mod;
        a=(a<<1)%mod;
    }
    return ans;
}

LL pow_mod(LL a, LL b,LL mod)
{
    LL ans=1%mod;
    for(;b;b>>=1)
    {
        if(1&b)ans=mul_mod(ans,a,mod);
        a=mul_mod(a,a,mod);
    }
    return ans;
}

拓展欧几里得（拓展欧几里得算法目的是求x,y满足ax+by==gcd(a,b)）

如果gcd(a,b)==1（互质），则x是a膜b的乘法逆元（ax%b==1）

typedef long long LL;
//拓展欧几里得算法目的是求x,y满足ax+by==gcd(a,b)
LL exgcd(LL a,LL b,LL &x, LL &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    else
    {
        LL _gcd=exgcd(b,a%b,y,x);
        y=y-(a/b)*x;
        return _gcd;
    }
}

在这里插入图片描述

乘法逆元（都要求a和mod互质）

乘法逆元、扩展欧几里得算法、二元一次方程、a的n次方取余
 逆元详解

若ax≡1 mod f, 则称a关于1模f的乘法逆元为x。也可表示为ax≡1(mod f)。

当a与f互素时，a关于模f的乘法逆元有解。如果不互素，则无解。如果f为素数，则从1到f-1的任意数都与f互素，即在1到f-1之间都恰好有一个关于模f的乘法逆元。
https://baike.baidu.com/item/%E4%B9%98%E6%B3%95%E9%80%86%E5%85%83/5831857?fr=aladdin
- 拓展欧几里得（要求互质）
```
LL get_inverse_ele(LL a,LL mod)
{
	if(gcd(a,mod)!=1)return -1;
	LL x,y;
	exgcd(a,mod,x,y);
	return (x+mod)%mod;
}
```
- 费马小定理（要求mod为质数，且互质）
  
  a^(p-1)%p==1，逆元为a^(p-2)
```
LL get_inverse_ele(LL a, LL mod)
{
	return pow_mod(a,  mod-2, mod);
}
```
- 另一个mod为质数时的逆元求法
  
  （要求：函数传入参数a<p, a和p互质，p为质数）
  
  事实上这个公式p不为质数也有可能生效。比如p=4和a=3。
  
  其实应该这么讲：p为素数时，一定生效。
  
  看下面的证明可知，这是一个递归的的证明，即要求inv(a)则应该先知道inv(p%a)，换句话说，inv(a)存在时，要求inv(p%a)存在，由此可知p和a，p%a, p%(p%a)，p%(p%(p%a)), …，1（到1停止）互质。一个合数无法保证这一点。但如果p为质数，a<p，则a和p一定互质（当然是在正整数范围内讨论，0单独讨论）
  
  除法取模与逆元/费马小定理
```
#include<cstdio>
typedef long long LL;
LL inv(LL t, LL p) 
{//求t关于p的逆元，注意:t要小于p，最好传参前先把t%p一下 
    return t == 1 ? 1 : (p - p / t) * inv(p % t, p) % p;
}
int main()
{
    LL a, p;
    while(~scanf("%lld%lld", &a, &p))
    {
        printf("%lld\n", inv(a%p, p));
    }
}
```
  逆元的几种求法（扩展欧几里得，费马小定理或欧拉定理，特例，打表等）
- O(n)求多个数逆元
  除法取模与逆元/费马小定理
```
#include <cstdio>
const int N = 200000 + 5;
const int MOD = (int)1e9 + 7;
int inv[N];
int init()
{
    inv[1] = 1;
    for(int i = 2; i < N; i ++)
        inv[i] = (MOD - MOD / i) * 1ll * inv[MOD % i] % MOD;
}
int main()
{
    init();
}
```
(a/b)%mod

由于取模运算只对加，减，乘有分配率，除法取模就显得麻烦了。

一般来说，乘法和加法之所以要取模，是因为结果很大，除法不存在这个但问题，一般来说要除法取模的情况会出现在中间步骤中，比如等比数列求和使用通项公式时，使用了除法，这破坏了分配率。
- b和mod互质
  可以用逆元解决。
- b和mod不互质
  我曾经对这个问题挣扎过
  
  （下面这个是错的）
  首先，先让分子除以分母为整数，即先a = a - a % b
  这时，a,b,mod必有一个最大公约数gcd，它同时也是b和mod的最大公约数。
  那么应该想到
  ans = ( (a/gcd) / (b/gcd) ) % mod
  但是这样写是错的，因为(b/gcd)和(mod/gcd）互质，不意味着(b/gcd)和(mod)互质（举个例子 b=8，mod=20）。
  不过理清楚后，思路也简单了，只要递归求当前值与mod的gcd直到gcd为1，把之前所有gcd累乘，得到super_gcd，就能ans = ( (a/super_gcd) / (b/super_gcd) ) % mod，此时按照互质来解。
  
  实际上这样又回到了原点。
  
  (b/super_gcd)可以求出逆元，但又多了一个(a/super_gcd)这个分母还是与mod不互质。
  
  实际上换一个思路就可以得到一种通用的解法（b可以整除a时）：
  
  (转自博客：逆元详解)
  
  这个已经不算什么逆元了。
  其实是巧妙地回避了除法取模。
${\frac{a}{b}}\bmod{m}=({a}\bmod{bm})/b$

至于b不整除a时，其实还是一样的（C/C++除法是向0取整，对正数而言是向下取整），只要在证明过程中改几步

$a\bmod b \\ a \bmod bm = bx+ a \bmod b \\ a \bmod bm = bx+ (a \bmod bm) \bmod b \\ 记 g=(a \bmod bm) \\ 则 x=(g-g \bmod b)/b= \lfloor \frac{a \bmod bm}{b} \rfloor$

这说明，即使不整除，在C++里也是一个写法。

不过要注意的是bm可能会很大，如果出现溢出的情况，还是不能这么做。

欧拉函数

https://blog.youkuaiyun.com/qq_36409190/article/details/53173777

在数论，对正整数n，欧拉函数是小于或等于n的正整数中与n互质的数的数目（φ(1)=1）。此函数以其首名研究者欧拉命名(Euler’s totient function)，它又称为Euler’s totient function、φ函数、欧拉商数等。例如φ(8)=4，因为1,3,5,7均和8互质。从欧拉函数引伸出来在环论方面的事实和拉格朗日定理构成了欧拉定理的证明。

//直接求解欧拉函数  
#include<cstdio>
int euler(int n){ //返回euler(n)   
    int res=n,a=n;  
    for(int i=2;i*i<=a;i++){//从小到大尝试n的质因数 
        if(a%i==0){//如果i是n的质因数 
            res=res/i*(i-1);//提了一个1/i出来，先进行除法是为了防止中间数据的溢出   
            while(a%i==0) a/=i;//欧拉函数只记算一种质因数 
        }  
    }  
    if(a>1) res=res/a*(a-1);//如果最后还剩因子 
    return res;  
}
int main(){
    int x;
    scanf("%d",&x);
    printf("%d",euler(x));
    return 0;
}

//筛选法打欧拉函数表   
#include<cstdio>
#define Max 1000001  
int euler[Max];  
void Init(){
    euler[1]=1;
    for(int i=2;i<Max;i++)
        euler[i]=i;
    for(int i=2;i<Max;i++)
        if(euler[i]==i)//如果i是质数 
            for(int j=i;j<Max;j+=i) 
                euler[j]=euler[j]/i*(i-1);//提一个1/i,先进行除法是为了防止中间数据的溢出
    return ;
}
int main(){
    Init();
    for(int i=1;i<=100;i++)
        printf("%d\n",euler[i]);
    return 0;
}

素数打表

根据素数定理，N以下的素数约有N/lnN个

正常的打表法用的多，也很简单，就不多说了，现在来看几种高效的打表法：

1和2来自李煜东的《算法竞赛进阶指南》

Eratosthenes筛法（最常用，O(nloglogn)）
复杂度计算看这里：Sieve of Eratosthenes

bool v[MAX_N];//合数标记，true为合数
void primes(int n){
	memset(v,false,sizeof(v));
	for(int i=2;i<=n;i++){
		if(v[i]) continue;
		//从i*i开始筛，可以证明小于这个数的i的倍数已经被筛掉了
		for(int j=i;j*i<=n;j++) v[i*j]=true;
	}
}

欧拉线性筛（时间复杂度为O(n)）
每个数只会被它最小的质因数筛掉，记每个数最小质因数为第ppi个质数，prime[j]*i>n时第二层循环数为qqi，则第二层循环数为mi=min(ppi, qqi)，则时间复杂度为m2+m3+…+mn=O(n)

const int N = 1000000;
int v[N], prime[80000]; //要注意保证pirme[]数组足够存储[2,N]之间的素数（不要低于N/lnN）
int m;
int primes(){
	memset(v,0,sizeof(v)); //最小质因子
   	m = 0;//质数数量
   	for(int i=2;i<=n;i++){
   		if(v[i]==0) {v[i]=i; prime[++m]=i;} //i是质数
   		//给当前i乘上一个质因子
   		for(int j=1;j<=m;j++){
   			//有比prime[j]更小的质因子，或者prime[j]*i超出n范围，停止循环
   			//还有一种 if(i%prime[j]==0)break 的写法
   			if(prime[j]>v[i] || prime[j]*i>n) break;
   			v[i*prime[j]]=prime[j];
   		}
   	}
}

只存素数

/*=================================*\
  素数打表
  该函数执行后在prim[]数组中存入
  从2开始的从小到大的numOfPrim个素数
\*=================================*/
const int numOfPrim = 1000;
int prim[numOfPrim] = {2,3};
void prime(){
    int tally=2;
    bool flag;
    for(int i=5 ; tally<numOfPrim ; i+=2){
        flag = true;
        for(int j=0 ; prim[j]*prim[j]<=i ; j++)
            if( i%prim[j]==0 ){ flag = false; break; }
        if( flag ){
            prim[tally]=i;
            tally++;
        }
    }
}

二分模板

下面是整数二分模板
注意>>1是向下取整，/2是向0取整

//找>=x最小的一个，搜索[min, max]时不会搜到max
while(l < r)
{
	int mid= (l + r) >> 1;
	if(a[mid] >= x) r = mid;
	else l=mid + 1;
}
return a[l];

//找<=x最小的一个，搜索[min, max]时不会搜到min
while(l < r)
{
	int mid = (l + r + 1) >> 1;
	if(a[mid] <= x) l = mid;
	else r = mid - 1;
}
return a[l];

质因数分解

#include <iostream>
#include <vector>
#include <utility>

typedef long long LL;
using namespace std;

//质因数分解
//x>=2
vector<pair<LL,LL> > get_prime_factors(LL x)
{
    vector<pair<LL, LL> > prime_factors;
    for(LL i = 2; i*i<=x; i++)
    {
        if(x % i)continue;
        LL cnt = 0;
        while(x % i == 0)
        {
            x /= i;
            cnt++;
        }
        prime_factors.push_back(make_pair(i, cnt));
    }
    //注意，最后可能会留一个大于n^(1/2)的质数，n为x原本的大小
    //例如14=2*7
    if(x!=1)prime_factors.push_back(make_pair(x,1));
    return prime_factors;
}

位运算

下面是一些比较基本的

//以下是在int（有符号整数）的情况下讨论的，如果要unsigned，有的要有略微变动
n>>1 //除以2并且向下取整，注意/2是向0取整
n<<1 //乘以2
1&n //判断奇偶，1奇0偶
(n>>k)&1 //取第k位
n & ((1<<k)-1) //取0~(k-1)位，或者说对2^k-1取模
n ^ (1<<k) //第k位取反
n | (1<<k) //第k位置1
n & (~(1<<k)) //第k位置0
n & (-n) //n>0，lowbit运算，lowbit(01101000)=00001000 （即取最右边的1构成的数）
(x ^ y)>=0 //同号时，为true（即x^y>=0），不过这里是把0当正数看的。如果不想这样，可以额外加特判

位运算其实有很多讲的，其实可以看下面这两个链接（都是斯坦福的那一篇）。里面有很多神奇的技巧，比如位交错(interleave)，位反转(reverse)，打表求整数log，判断整数是否为2的幂等等。我觉得这些很多是应该掌握的。

http://graphics.stanford.edu/~seander/bithacks.html
https://www.cnblogs.com/Atela/archive/2011/04/27/2030393.html

要以2^k为索引时的打表法（哈希）

因为2^k可能很大，不能直接用数组。但是有个hash法可以非常巧妙地解决这个问题。对32位的，以37为模；对64位的，以67为模（其实都可以以67为模；也不一定是这几个数，只要大于最大值且满足一定条件就可以）
```
//举个例子：H[2^k]=k----->H[(1<<k)%67]=k
int H[67] = {0};
for(long long i = 0;i < 60;i++)
	H[ (1ll << i) % 67] = i;
```

计算整数log_2

如果要计算一个正整数的log_2向下取整结果，如果数大（如64位的）且要求准确，显然不能直接用cmath的int ans=log(x)/log(2)，首先如果数大了，会有精度损失的问题，其次，log效率其实不是很高（相对来说不高，但是一般可以忽略不计）。
实际上，整数计算时，求取log_2(x)其实就是求二进制下x最高位的1所在的位数。遗憾的是，非汇编写法下，没有很方便的方法求取（最低位的可以通过 (x^(-x))配合H[(1ll<<k)%67]=k表在O(1)时间求出来））。

//朴素方法
inline int log2_int(register int x){
    register int ans = 0;
    while((x >>= 1)) ++ ans;
    return ans;
}

//打表法
Log[0]=-1;
for(int i=1;i<=n;i++)Log[i]=Log[i>>1]+1;

//据说这是lua原码
//https://blog.youkuaiyun.com/w1964332/article/details/16923337
//http://www.bkjia.com/Linuxjc/774385.html
//道理很简单，就是打表法结合朴素法。
//lua这段代码的作者认为超过256的就不常用了，平常用就是一次索引获得结果。
//其实下面代码可以不用256而使用其他值（如65536）
//这样，就得掌握log2的打表方法
int luaO_log2 (unsigned int x) { 
  static const lu_byte log_2[256] = { 
    0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, 
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, 
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 
    8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, 
    8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, 
    8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, 
    8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8 
  }; 
  int l = -1; 
  while (x >= 256) { l += 8; x >>= 8; } 
  return l + log_2[x]; 
}

//据说这个是汇编写法，复杂度O(1)。没研究过
//基本没什么兼容性
//https://blog.youkuaiyun.com/qq_35703773/article/details/79167210
unsigned int LOG2(unsigned int x){
    unsigned int ret;
    __asm__ __volatile__ ("bsrl %1, %%eax":"=a"(ret):"m"(x));
    return ret;
}