继续实现算法导论里面的题目和算法

1.动态规划:切钢条问题

#算法导论第四章:动态规划题目:切钢条
'''问题:给定一个长度为n的钢条,和一个价格表,问切割方案s.t.受益最大
切割表:
长度:1 2 3 4 5  6  7   8  9 10
价格:1 5 8 9 10 17 17 20 24 30

'''
#作为动态规划的本质还是要先写出地柜的版本.然后改地推即可.
#首先做出数据:
a=[]
a.append([1,2,3,4,5,6,7,8,9,10])
a.append([1,5,8,9,10,17,17,20,24,30])
print (a)
#先写一个水的版本,把方案利润最后的答案打印出来即可.
def fangan(n):
 tmp=0#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=(a[1][i])
                break
 for i in range(1,n//2+1):
      b=fangan(i)+fangan(n-i)
      if b>tmp:
       tmp=b
 return tmp
print (fangan(5))

#下面写把切割方案也写出来.这个书中没有体现,感觉有点难的!
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege
 return tmp
print (fangan(10))
#然而利用函数的多返回值还是成功写出来了.注意这里的返回值用list实现,用tuple修改不了所以不行.


#然而如果老板要求获得所有最好的方案呢?比如fangan(7)就应该有2个解7=1+6=2+2+3



#继续实现上面的目标:
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege

 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b==tmp[0]:
          tmp.append(sorted(qiege))
 #去重
 k=[]
 for i in range(2,len(tmp)):
         
         if tmp[i] in tmp[:i]:
             
              k.append(i)#记录需要弹出的index


 #不知道为什么这个去重有bug!??但是基本不难了
 return tmp
print (fangan(7))
View Code

 更新版:改了上面的一个bug

#算法导论第四章:动态规划题目:切钢条
'''问题:给定一个长度为n的钢条,和一个价格表,问切割方案s.t.受益最大
切割表:
长度:1 2 3 4 5  6  7   8  9 10
价格:1 5 8 9 10 17 17 20 24 30

'''
#作为动态规划的本质还是要先写出地柜的版本.然后改地推即可.
#首先做出数据:
a=[]
a.append([1,2,3,4,5,6,7,8,9,10])
a.append([1,5,8,9,10,17,17,20,24,30])
print (a)
#先写一个水的版本,把方案利润最后的答案打印出来即可.
def fangan(n):
 tmp=0#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=(a[1][i])
                break
 for i in range(1,n//2+1):
      b=fangan(i)+fangan(n-i)
      if b>tmp:
       tmp=b
 return tmp
print (fangan(5))

#下面写把切割方案也写出来.这个书中没有体现,感觉有点难的!
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege
 return tmp
print (fangan(10))
#然而利用函数的多返回值还是成功写出来了.注意这里的返回值用list实现,用tuple修改不了所以不行.


#然而如果老板要求获得所有最好的方案呢?比如fangan(7)就应该有2个解7=1+6=2+2+3



#继续实现上面的目标:
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege

 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b==tmp[0]:
          tmp.append(sorted(qiege))
 #去重
 k=[]
 for i in range(2,len(tmp)):
         
         if tmp[i] in tmp[:i]:
             
              k.append(i)#记录需要弹出的index

## for i in k:
##     tmp.pop(i)
 #不知道为什么这个去重有bug!??但是基本不难了知道错哪里了,pop一次index就变了!!!
 #所以还是要用for来重新建立就行了.
 c=[]
 for i in range(len(tmp)):
     if i not in k:
         c.append(tmp[i])
 return c
print (fangan(7))
View Code

 继续修改成动态规划版本,并且用一个储存技巧,轻松实现递归到动态规划的转化

#算法导论第四章:动态规划题目:切钢条
'''问题:给定一个长度为n的钢条,和一个价格表,问切割方案s.t.受益最大
切割表:
长度:1 2 3 4 5  6  7   8  9 10
价格:1 5 8 9 10 17 17 20 24 30

'''
#作为动态规划的本质还是要先写出地柜的版本.然后改地推即可.
#首先做出数据:
a=[]
a.append([1,2,3,4,5,6,7,8,9,10])
a.append([1,5,8,9,10,17,17,20,24,30])
print (a)
#同样对他进行计数
#先写一个水的版本,把方案利润最后的答案打印出来即可.
count=0
def fangan(n):
 global count
 count+=1
 tmp=0#用来储存所有的方案的利益
 
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=(a[1][i])
                break
 for i in range(1,n//2+1):
      b=fangan(i)+fangan(n-i)
      if b>tmp:
       tmp=b
 return tmp
print (fangan(7))
print ('递归版本的计数是:')
print (count)
#下面写把切割方案也写出来.这个书中没有体现,感觉有点难的!
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege
 return tmp
print (fangan(10))
#然而利用函数的多返回值还是成功写出来了.注意这里的返回值用list实现,用tuple修改不了所以不行.


#然而如果老板要求获得所有最好的方案呢?比如fangan(7)就应该有2个解7=1+6=2+2+3



#继续实现上面的目标:
def fangan(n):
 tmp=[0,[0]]#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=[a[1][i],[n]]
                
                break
 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b>tmp[0]:
       tmp[0]=b
       tmp[1]=qiege

 for i in range(1,n//2+1):
      b=fangan(i)[0]+fangan(n-i)[0]
      qiege=fangan(i)[1]+fangan(n-i)[1]
      if b==tmp[0]:
          tmp.append(sorted(qiege))
 #去重
 k=[]
 for i in range(2,len(tmp)):
         
         if tmp[i] in tmp[:i]:
             
              k.append(i)#记录需要弹出的index

## for i in k:
##     tmp.pop(i)
 #不知道为什么这个去重有bug!??但是基本不难了知道错哪里了,pop一次index就变了!!!
 #所以还是要用for来重新建立就行了.
 c=[]
 for i in range(len(tmp)):
     if i not in k:
         c.append(tmp[i])
 return c
print (fangan(7))

#想法是直接在递归过程中,建立表,然后每一次递归调用都先搜索表中是否存在,
#从而实现动态规划.下面来实现.



print ('下面是动态规划版本')

#首先还是做这个只返回最大利益的这个简单例子.b作为表格来储存每一个递归结果.
#简单的写出来,然后为了讨论效率加一个count来计数
b=[0]*len(a[0])
count=0
def fangan(n):
 global count
 count+=1
 global b
 if b[n-1]!=0:
     return b[n-1]
 tmp=0#用来储存所有的方案的利益
 if n in a[0]:
      for i in range(len(a[0])):
           if a[0][i]==n:
                tmp=(a[1][i])
                break
 for i in range(1,n//2+1):
      c=fangan(i)+fangan(n-i)
      if c>tmp:
       tmp=c
 
 b[n-1]=tmp
 return tmp
print (fangan(7))
print('动态规划版本的计数是')
print (count)
View Code

 2.贪心问题:会议安排问题

'''贪心算法:算法导论第十六章'''

#活动选择问题:
'''

题目是:下面的表表示会议i在si点开始在fi点结束.那么如何安排会议可以让
会议开的数目最大

i  1    2   3    4   5    6   7    8   9    10  11

si 1    3   0   5   3   5   6   8   8   2   12

fi 4    5   6   7   9   9   10  11  12  14  16




'''

a=[]

a.append([1,2,3,4,5,6,7,8,9,10,11])
a.append([1,3,0,5,3,5,6,8,8,2,12])


a.append([4,5,6,7,9,9,10,11,12,14,16])

print (a)
#先写递归版本:理清思路再考虑效率问题.基本是一个On^3的

def shuchu(n):#函数代表如果会议只让最晚开到n点,包含n点的话,那么
              #最大的会议数量是几?这么写显然不行,因为如果一个会议已经用过了
               #但是你还是记录他的一个上届来做n显然会重复,
             #这时候的思想就是细分,加强这个条件,来让子问题不重叠.
            #类似上升子序列的思路:把shuchu(n)这个函数的n定义改成:
        #严格的开到n点,这个情况下最大的会议数量.例如比如不考虑这个题目的话
       #shuchu(5.5) return一个0.
    c=[0]
    
    for i in range(1,n):
        d=0
        for j in range(len(a[0])):
            if a[1][j]>=i and a[2][j]==n:
                c.append(shuchu(i)+1)
            else:
                c.append(shuchu(i))
    return max(c)     
            
print (shuchu(6))
 







#改成动态规划:上面实在是太卡了
b=[-1]*2*len(a[0])
def shuchu(n):#函数代表如果会议只让最晚开到n点,包含n点的话,那么
              #最大的会议数量是几?这么写显然不行,因为如果一个会议已经用过了
               #但是你还是记录他的一个上届来做n显然会重复,
             #这时候的思想就是细分,加强这个条件,来让子问题不重叠.
            #类似上升子序列的思路:把shuchu(n)这个函数的n定义改成:
        #严格的开到n点,这个情况下最大的会议数量.例如比如不考虑这个题目的话
       #shuchu(5.5) return一个0.
    global b
    if b[n]!=-1:
        return b[n]
    c=[0]
    
    for i in range(1,n):
        d=0
        for j in range(len(a[0])):
            if a[1][j]>=i and a[2][j]==n:
                c.append(shuchu(i)+1)
            else:
                c.append(shuchu(i))
    b[n]=max(c)
    return b[n]   
            
print (shuchu(16))

#贪心算法:直接O(n) 每一次选结束时间最早的并且跟前面无冲突的即可,然后扫一遍数据就完事了...
View Code

 3.用indegree方法来求解有向图的拓扑排序问题

'''
图的表示:算法导论22章
图又2中表示方式:
1.是链表:一个节点跟他相连的节点都用一个链表穿起来.
         然后n个节点就用n个链表来储存这个图的信息
2.是矩阵:用节点数n,n*n的矩阵来储存,两个点如果有边就存1.
优势和劣势:
1.链表:适用于稀疏图,但是很难判断一个边是否存在
2.矩阵:反之.

'''
'''
拓扑排序

一.定义

对一个有向无环图(Directed Acyclic Graph简称DAG)G进
行拓扑排序(Topological Sort),是将G中所有顶点排成一个线
性序列,使得对图中任意一对顶点u和v,若<u,v>∈E(G),则u在
线性序列中出现在v之前。通常将这样的线性序
列称为满足拓扑次序(Topolgical Order)的序列,简称拓扑序列。'''




'''
用链表来拓扑排序
入度是图论算法中重要的概念之一。它通常指有向
图中某点作为图中边的终点的次数之和。'''


'''
入度表

对于DAG的拓扑排序,显而易见的办法!!!!!!!!
这是一个好方法:

找出图中0入度的顶点;
依次在图中删除这些顶点,删除后再找出0入度的顶点;
然后再删除……再找出……
直至删除所有顶点,即完成拓扑排序
为了保存0入度的顶点,我们采用数据结构栈(亦可用队列);算法的可视化可参看这里。

图用邻接表(adjacency list)表示,用数组inDegreeArray[]记录结点的入度变化情况。C实现:'''


class Graph():
    def __init__(self,n):
        self.vertex_num=n
        self.edge=[]
        for i in range(n):#2-dim-empty list cannot use *99
         self.edge.append([])
        self.indegree=[0]*n
#vertex_index=0,1,2,3,4....
    def addedge(self,n,m):
        self.edge[n].append(m)
        self.indegree[m]+=1
    def topological_sort(self):
        a=[]
        b=len(self.indegree)

        while (b):
         for i in range(len(self.indegree)):
            if 0 not in self.indegree:
                return False
            if self.indegree[i]==0:
                print (i)
                a.append(i)
                self.indegree[i]-=1
                
                for j in self.edge[i]:
                   
                    self.indegree[j]-=1
                b-=1
                break
         
        return a










            
g=Graph(6)
g.addedge(5,2)
g.addedge(5,0)
g.addedge(4,0)
g.addedge(4,1)
g.addedge(5,4)
g.addedge(0,5)
g.addedge(2,3)
g.addedge(3,1)
print (g.edge)
print (g.indegree)
print (g.topological_sort())
View Code

 4.数论,素数,线性规划,密码的学习

'''
下面学习算法导论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、246、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
View Code

 5.计算几何:

'''
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))


























        
View Code

 

'''
计算几何记录下思路吧:
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贷款里面的点.
'''
View Code

 

posted on 2018-03-15 22:19  张博的博客  阅读(449)  评论(0编辑  收藏  举报

导航