'''
下面学习算法导论29章线性规划问题.'''
'''
解决方法:1.单纯型法:是一个最坏情况指数型的算法
2.内点法.:是一个多项式级数的算法
3.整数线性规划:nphard问题.没有找到多项式的时间算法
标准型和松弛型:标准型中,所有的约束都是不等式,而在松弛型中所有的约束都是等式,
并且要求变量非负.
非常有用的转化思想:
1.我们可以吧单源最短路径问题转化为一个线性规划.
'''
'''
单纯型法:就是化成松弛型矩阵之后
的等解变换.然后逐渐把目标函数的变量前面的系数,都修改成负数,这样,最后的前面的
常数项就是我们要的答案.'''
#算法很复杂就不写代码了
'''
30章,多项式乘法和傅里叶变换.
1.多项式的表达:1.系数的表达:比较trivial
2.点值表达:一个次数界为n的多项式A(x)的点值表达是由
{(x0,y0),...,(xn-1,yn-1)}来表达的集合
其中yk=A(xk)'''
'''
通过系数表达和点值表达来进行多项式乘法,和在复平面的单位圆里面取点值表达.来操作的'''
'''
28章:矩阵计算.
1.线性方程组:用LUP分解,也就是高斯消元法.
31章:经典算法:
1.密码里面有一个问题:求解模线性方程:ax≡b(mod n)
2.素数定理:素数跟自然数的个数是logn:n 经典问题证明非常难.
3.工具:欧拉函数φ(n)=n*∏(1-1/p)其中p是n的素因子:
这个函数表示不大于n的整数有多少个跟n互素.
4.对求模运算有一个经典的问题:如果gcd(a,n)=1,那么方程ax≡1(mod n)对模n有唯一的解,
叫做这个解释a对模n的乘法逆元.如果gcd(a,n)!=1,那么方程无解.
这就是ras加密的原理.
5.拓展欧几里得算法:设d=gcd(a,n):那么d=ax+ny.需要这个算法给a,n输出d,x,y:因为这些量
对模的求逆运算很重要.
'''
def Extended_Euclid(a,n):#假设a大n小.
if a<n:
x,y,z=Extended_Euclid(n,a)
return (x,y,z)
if n==0:
return (a,1,0)
else:#下面的步骤需要一个反,也就是辗转一下.
(x,y,z)=Extended_Euclid(a % n,n)
return (x,y,z-(a//n)*y)
print (Extended_Euclid(1423423403242342300,13477088832423432423432000))
'''
书中定理:
1.令d=gcd(a,n),假设对某些整数x1,y1有d=a*x1+n*y1
那么有如果d|b那么方程ax≡b(mod n)才有解.并且
解是:xi=x0+i(n/d)这里i=0,1,...,d-1
x0=x1(b/d)mod n'''
#证明带入验证即可.所以下面利用这个定理可以直接给出模方程的算法.
def modequation(a,b,n):#方程ax≡b(mod n)
d,x1,y1=Extended_Euclid(a,n)
if b%d==0:
c=[]
c0=x1*(b//d) %n
c.append(c0)
for i in range(1,d):
c0+=(n//d)
c0=c0%n
c.append(c0)
return c
else:
return None
print (modequation(14,30,100))
#利用这个modequation可以高效的求出模的逆元.比如下面的例子.
print (modequation(7,1,10)) #输出答案3
'''
元素的幂的模操作:
1.欧拉定理:对于任意整数n>1,有a的φ(n)次幂≡1(mod n)对所有的a小于n都成立.
2.费马定理:如果p是素数,则a的p-1次幂≡1(mod p)对所有的a小于p都成立.
'''
'''
希望有一个高效的求a的b次幂mod n的值.下面给出.利用二进制来实现.
'''
def modular_expotion(a,b,n):#高效的求a的b次幂mod n的值
#本质也就是把每一次分解成平方和平方乘a运算,每一次用mod来加速运算.
b=bin(b)
b=b[2:]
d=1
for i in range(len(b)):
if b[i]=='0':
d**=2
d%=n
if b[i]=='1':
d**=2
d%=n
d*=a
d%=n
return d
print (modular_expotion(751221326,5604675868673545300,5634231)) #秒算出3496792
'''
然后直接算幂再mod的话,直接电脑卡死
'''
'''
rsa秘钥加密系统总结:
在rsa秘钥加密系统中,一个参与者按下列过程创建他的公钥和秘钥:
1.随机选取两个大的素数p和q,使得p≠q,例如素数p和q都取1024位.
2.计算n=pq
3.选取一个与φ(n)互质的小的奇数e.(其中φ(n)=(p-1)(q-1)). 欧几里得算最大公约数很快.
另外从书中我们知道自然数n和φ(n) 比是lnln n :1 所以跟他互素的数非常多.
从小到大遍历我们轻松算出是14:1.所以基本不超过尝试100次左右就能找到这个小奇数e
并且因为是lnln 所以就算p,q取多大比如100万位也都秒算.
4.对模φ(n),计算出e的乘法逆元d的值.
5.将对P=(e,n)公开,并作为参与者的RSA公钥.
6.将对S=(d,n)保密,并作为参与者的RSA秘钥.
对于上面的方案,消息M是一个数字.
把P(M)定义为M**e mod n.
把S(C)定义为C**d mod n.
所以容易知道Alice可以通过这两个变换把M变回来在mod n的意义下.
所以只需要M这个10进制的数的位数*e小于6即可.感觉很迷.这个限制也太严格了.
使用流程:
1.Bob给Alice发信息:发P(M)
2.Alice接受到然后用S转化即可.
'''
'''
伪素数:如果n是一个合数,并且a**(n-1)≡1(mod n)我们就称n是一个基为a的伪素数.
'''
def pseudoprime(n):
if modular_expotion(2,n-1,n)%n!=1:
return 'composite'
else:
return 'prime'
print (pseudoprime(857))#但是无论加多少个基都无法保证psudo变成真素数.因为有carmicheal数的存在.
#但是1024位的话,基于2的伪素数有1/10的40次幂还要少的概率.
#利用这个就可以随机生成一个大素数.根据素数定理:只需要测试ln n个数即可找到一个prime.
def find_a_big_prime(n):#n要保证很大,比如一亿多效果才好,下面为了测试需要就弄个小的数.
#因为n越大,伪素数的出现概率越低.
#这个函数输入n,返回一个比n大的伪素数.
import math
for i in range(2*int(math.log(n))):
n+=1
if pseudoprime(n)=='prime':
return n
print(find_a_big_prime(560))#这里面故意写一个范例,输出561,他是一个camichale数
#他是3的倍数,所以他不是一个素数.
'''
一些实用的关于整除的trick:
1.一个整数被3整除 等价于 他的各个数位的数字相加被3整数.(证明显然)
2.一个整数被9整除 等价于 他的各个数位的数字相加被9整数.(证明显然)
3.一个整数被5整除 等价于 他的个数位的数字=0,5.(证明显然)
4.设该数有两位以上,若个位数+十位数*2能被4整除,那么该数能被4整除.
5.被7整除等价于:若一个整数的个位数字截去,再从余下的数中,减去个位数的
2倍,如果差是7的倍数,则原数能被7整除.如果差太大或心算不易看出是否7的倍数
, 就需要继续上述「截尾、倍大、相减、验差」的过程,直到能清楚判断为止.(证明显然)
6.能被6整除的数的特征末尾是0、2、4、6、8且各位上数字的和能被3整除
7.显然1000能被8整除,所以1000的倍数能被8整除,所以看一个数能否被8整除,只需看它的后3位.
'''
print (322223333)
#下面通过构造witness函数把pseudoprime进行一次重写.得到更高的效率.
def witness(a,n):#函数目标是验证n是否是以a为基地的伪素数.
#返回True就表示是合数,False代表素数
if n %2==0:
return True
#先把n-1进行2的次方分解:
beifenjie=n-1
t=0
for i in range(beifenjie):
if beifenjie%2==0:
t+=1
beifenjie//=2
else:
break
u=beifenjie# 得到了n-1=2**t*u
x0=modular_expotion(a,u,n)
#下面这行的操作在算法导论第三版chs.pdf里面p568页写的,实在是太漂亮的结果了.
#大概写一下,具体应该翻书多看看怎么判断的合数.
#思路大概是这样,你经过t次平方得到的数.然后如果你已经发现xi-1是一个以n为模1的
#非平凡平方根,因为xi!=+-1(mod n),但是xi=1 (mod n),所以n必然是合数,因为
#已经证明过了只有n是合数时候才可以存在以n为模的1的非平凡平方根.(非平凡就是非1)!!!!!!
#我直接就惊了!用这个操作居然能排除很多calmicakal数.比如561.下面写代码.
x_old=x0
for i in range(t):
x_new=x_old**2%n
if x_new==1 and x_old!=1 and x_old!=n-1:
return True
x_old=x_new
if x_old!=1:
return True
return False
'''
下面就是简单的miller-rabin算法:
效率非常高,然后几乎保证了prime的判断是严格的!!!!!!!!
'''
#下面函数就是用s个基地来判断n是否是素数的终极程序
def miller_rabin(n,s):
#经过这个代码发现import 一个库最好少用
#比如我如果用numpy库会明显感觉卡,用random明显速度快.
#我估计他是把代码都贴过来,所以少import没用的库会加速很多.
import random
for j in range(s):
j=random.randint(1, n-1)
if witness(j,n):
return 'Composite'
return 'prime'
print (miller_rabin(2336,5))
'''
下面的是开放问题整数的因子分解.
也就是破译ras的算法问题.没有高速方法,只有启发式方法.
时间未知,也就是基本不是多项式时间的.
'''
#感觉很神秘,解释也比较奇怪的算法.
def pollard_rho(n):
import random
#注意rand这个函数里面的上下标都是闭括号.都能取到.
x=random.randint(0, n-1)
y=x
while 1:
x=(x**2-1)%n
d=Extended_Euclid(y-x,n)[0]
if d!=1 and d!=n:
print (d)
import random
'''
33章:高大上的计算几何. #顺便学习了运算符的重载.
1.线段的操作:
1.向量的叉积:正负,用来判断两个的相对顺时针和逆时针.
等于0表示2个向量共线.
'''
#下面写一个判断2个线段相交的函数.下面都是在R2中讨论,更高维的可以类似推广.
class point():
def __init__(self,x,y):
self.x=x
self.y=y
def __sub__(self,other):#other不用事先制定类型,可以直接当做一个point来用
#本身已经写好了调用的方法,所以还是用-号来调用即可.
return point(self.x-other.x,self.y-other.y)
#上面为了方便,我们引入运算符的重载.
def direction(p1,p2,p3):#这函数是求p1,p2这个向量和p1,p3这个向量的叉积.
#用行列式法则就能计算这种简单的2维情况.
return (p2-p1).x*(p3-p1).y-(p2-p1).y*(p3-p1).x #这个数如果大于0就表示从向量p1,p2转逆时针才能到p1,p3
#反之类似.
def segments_intersect(p1,p2,p3,p4):#函数的自变量可以取一个对象
#我们判断的线段是p1,p2为端点的 线段和p3,p4为端点的线段
#python 代码太长后面写即可 一行运行多个语句用,即可.
d1=direction(p1,p2,p3)
d2=direction(p1,p2,p4)
d3=direction(p3,p4,p1)
d4=direction(p3,p4,p2)
if direction(p1,p2,p3)*direction(p1,p2,p4)<0 and
direction(p3,p4,p1)*direction(p3,p4,p2)<0:
return True
if d1==0 and on_segment(p1,p2,p3):
return True
if d2==0 and on_segment(p1,p2,p4):
return True
if d3==0 and on_segment(p3,p4,p1):
return True
if d4==0 and on_segment(p3,p4,p2):
return True
return False
p1=point(3,4)
p2=point(3,5)
p3=point(2,5)
p4=point(4,-999)
def on_segment(a,b,c):
if min(a.x,b.x)<=c.x<=max(a.x,b.x) and min(a.y,b.y)<=c.y<=max(a.y,b.y):
return True
else:
return False
print (direction(p1,p2,p3)),print (324)
print (segments_intersect(p1,p2,p3,p4))
'''
计算几何记录下思路吧:
1.寻找凸包.1.graham方法:1.首先栈里面放入最低的点里面最左的点.作为原点
2.放入跟这个原点,夹角最小的2个点.append进去
3.继续找夹角更大一丁点的点,然后跟栈顶和栈倒数第二个栈顶比,
看是否弹出栈顶.
4.做到最后站里面的元素就是凸包的端点
速度O(Nlogn)
2.jarvis方法:1.相同
2.找夹角最小的和最大的.得到p1,p2两个点
3.再继续以p1找夹角最小的,和p2找夹角最大的.
4.循环即可.O(nh) h为凸包端点数
2.找最近的点对.时间O(nlogn)
1.分治法
2.合并时候还要多考虑中间地带2delta贷款里面的点.
'''