Add computaional method

master
mbinary 2018-10-02 21:24:06 +08:00
parent 0b1f116be9
commit 0c08851597
54 changed files with 1221 additions and 50 deletions

2
.gitignore vendored
View File

@ -1,7 +1,5 @@
.*
*.o
*.exe
tree_link.py
__pycache__/**
!.gitignore
todo

View File

@ -1,7 +1,7 @@
# Algorithm and Data Structures
>Notes and codes for learning algorithm and data structures :smiley:
Some pictures and idead are from `<<Introduction to Algotithm>>
Some pictures and ideas are from `<<Introduction to Algotithm>>
I use python 3.6+ and c++ to implements them.
Since I used f-Strings in python, you may use python 3.6+ to run the following python scripts.
@ -14,36 +14,67 @@ So,if you wannt to view the notes which contain latex math formulas and are in m
# Index
* [.](.)
* [notes](./notes)
* [alg-general.md](./notes/alg-general.md)
* [hashTable.md](./notes/hashTable.md)
* [red-black-tree.md](./notes/red-black-tree.md)
* [sort.md](./notes/sort.md)
* [tree.md](./notes/tree.md)
* [b-tree.md](./notes/b-tree.md)
* [graph](./notes/graph.md)
* [fibonacci-heap](./notes/fib-heap.md)
* [algorithm](./algorithm)
* [8Astar.py](./algorithm/8Astar.py)
* [cantor.cc](./algorithm/cantor.cc)
* [manacher.py](./algorithm/manacher.py)
* [markov.py](./algorithm/markov.py)
* [sort](./algorithm/sort)
* [sunday.py](./algorithm/sunday.py)
* [computationalMethod](./computationalMethod)
* [interplotion.py](./computationalMethod/interplotion.py)
* [iteration.py](./computationalMethod/iteration.py)
* [least_square.py](./computationalMethod/least_square.py)
* [linear_equation.py](./computationalMethod/linear_equation.py)
* [numerical_differential.py](./computationalMethod/numerical_differential.py)
* [numerical_integration.py](./computationalMethod/numerical_integration.py)
* [README.md](./computationalMethod/README.md)
* [solve-linear-by-iteration.py](./computationalMethod/solve-linear-by-iteration.py)
* [tongyu_equation.py](./computationalMethod/tongyu_equation.py)
* [vector_norm.py](./computationalMethod/vector_norm.py)
* [dataStructure](./dataStructure)
* [redBlackTree.py](./dataStructure/redBlackTree.py)
* [bTree.py](./dataStructure/bTree.py)
* [hashTable.py](./dataStructure/hashTable.py)
* [splayTree.py](./dataStructure/splayTree.py)
* [allOone](./dataStructure/allOone)
* [binaryHeap.py](./dataStructure/binaryHeap.py)
* [binaryTree.py](./dataStructure/binaryTree.py)
* [bTree.py](./dataStructure/bTree.py)
* [graph](./dataStructure/graph)
* [hashTable.py](./dataStructure/hashTable.py)
* [huffman](./dataStructure/huffman)
* [leftHeap.py](./dataStructure/leftHeap.py)
* [loserTree.py](./dataStructure/loserTree.py)
* [map.cc](./dataStructure/map.cc)
* [polynomial.cpp](./dataStructure/polynomial.cpp)
* [polynomial.py](./dataStructure/polynomial.py)
* [redBlackTree.py](./dataStructure/redBlackTree.py)
* [splayTree.py](./dataStructure/splayTree.py)
* [trie.py](./dataStructure/trie.py)
* [winnerTree.py](./dataStructure/winnerTree.py)
* [docs](./docs)
* [algorithm-general.md](./docs/algorithm-general.md)
* [b-tree.md](./docs/b-tree.md)
* [fib-heap.md](./docs/fib-heap.md)
* [graph.md](./docs/graph.md)
* [hashTable.md](./docs/hashTable.md)
* [README.md](./docs/README.md)
* [red-black-tree.md](./docs/red-black-tree.md)
* [sort.md](./docs/sort.md)
* [tree.md](./docs/tree.md)
* [dynamicProgramming](./dynamicProgramming)
* [lcs.hs](./dynamicProgramming/lcs.hs)
* [lcs.py](./dynamicProgramming/lcs.py)
* [matrix-multiply.py](./dynamicProgramming/matrix-multiply.py)
* [splitStripe.hs](./dynamicProgramming/splitStripe.hs)
* [splitStripe.py](./dynamicProgramming/splitStripe.py)
* [testVec2d.hs](./dynamicProgramming/testVec2d.hs)
* [Vec2d.hs](./dynamicProgramming/Vec2d.hs)
* [math](./math)
* [cantor.cc](./math/cantor.cc)
* [isPrime.py](./math/isPrime.py)
* [num_weight.py](./math/num_weight.py)
* [README.md](./README.md)
* [search](./search)
* [8Astar.py](./search/8Astar.py)
* [sort](./sort)
* [binaryTree.py](./sort/binaryTree.py)
* [heapSort.py](./sort/heapSort.py)
* [quickSort.py](./sort/quickSort.py)
* [radixSort.py](./sort/radixSort.py)
* [select.py](./sort/select.py)
* [shellSort.py](./sort/shellSort.py)
* [string](./string)
* [manacher.py](./string/manacher.py)
* [markov.py](./string/markov.py)
* [sunday.py](./string/sunday.py)

View File

@ -1,4 +0,0 @@
def adjustOrd(sizes):
''' adjust the chain-multiply of matrix, sizes=[row1,row2,..,rown,coln]'''
n = len(sizes)
if n<3: return

View File

@ -0,0 +1,35 @@
# 计算方法
>一些计算方法的算法,比如插值,拟合 近似计算,解方程组等
有些功能numpy, sympy可能已经就有,但是为了学习各种计算方法,我就重新写了一遍,主要使用的是numpy的数组,矩阵,sympy的符号运算
# 需要
* python3
* python modules
- sympy
- numpy
# 目录说明
## interplotion
插值, 有Lagrange插值, Newton插值
## iteration
迭代, 二分迭代, 不动点迭代,差商迭代, 弦截法迭代
## least_square
最小二乘拟合, 解矛盾方程组
## linear_equation
解线性方程组,用到
* doolittle分解
* crout分解
* ldlt分解
* 列主元消元法
## vector-norm
计算向量,矩阵的各种范数
## tongyu_equation
解同余方程
## LICENCE
[MIT](LICENCE.txt)
## 联系我
* mail: <img style="display:inline" src="http://ounix1xcw.bkt.clouddn.com/gmail.png"></img>
* QQ : 414313516

View File

@ -0,0 +1,84 @@
''' mbinary
#########################################################################
# File : interplotion.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
#########################################################################
# File : interplotion.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.github.io
# Github: https://github.com/mbinary
# Created Time: 2018-05-18 09:29
# Description: 插值计算,有牛顿插值,拉格朗日插值,以及通过插值得到的多项式估计新的函数值
#########################################################################
import sympy
from collections import namedtuple
from functools import reduce
from operator import mul
X = sympy.Symbol ('x')
point = namedtuple('point',['x','y'])
class interplotion:
def __init__(self,points):
self.points = [point(x,y) for x,y in points]
self.xs= [i for i,j in points]
self.poly,self.rem = self.newton(self.points,0,len(self.points)-1)
def newton(self,li,a,b):
'''li:[(x,f(x))...]'''
qs = [li[0].y]
def quoDiff(begin,end):
if begin == end:return li[begin].y
q = (quoDiff(begin+1,end)-quoDiff(begin,end-1))/(li[end].x-li[begin].x)
if begin == a:qs.append(q)
return q
quoDiff(a,b)
poly ,base = 0, 1
for i,q in enumerate(qs):
poly += q*base
base*= X-li[i].x
return poly, base*qs[-1]
def lagrange(self,points=None):
xs = None
if points is None:
xs = self.xs
points = self.points
else: xs =[x for x,y in points]
product = reduce(mul,[X-x for x in xs],1)
poly = 0
for x,y in points:
tmp = product/(X-x)
coef = y/(tmp.subs(X,x))
poly+= coef *tmp
return poly
def predict(self,val,poly = None):
if poly is None:poly = self.poly
return poly.subs(X,val) # note the func subs
if __name__ == '__main__':
f = interplotion([(81,9),(100,10),(121,11)])
p = f.lagrange()
print(p.subs(X,105))
print(p)
intor = interplotion([(0,11),(0.02,9),(0.04,7),(0.06,10)])
p = intor.lagrange()
print(p)
res = intor.predict(0.08)
print(res)

View File

@ -0,0 +1,108 @@
''' mbinary
#########################################################################
# File : iteration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
import sympy
import numpy as np
from math import sqrt
def newton(y:sympy.core,x0:float,epsilon:float=0.00001,maxtime:int=50) ->(list,list):
'''
newton 's iteration method for finding a zeropoint of a func
y is the func, x0 is the init x val: int float epsilon is the accurrency
'''
if epsilon <0:epsilon = -epsilon
ct =0
t = y.free_symbols
varsymbol = 'x' if len(t)==0 else t.pop()
x0= float(x0)
y_diff = y.diff()
li = [x0]
vals = []
while 1:
val = y.subs(varsymbol,x0)
vals.append(val)
x = x0- val/y_diff.subs(varsymbol,x0)
li.append(x)
ct+=1
if ct>maxtime:
print("after iteration for {} times, I still havn't reach the accurrency.\
Maybe this function havsn't zeropoint\n".format(ct))
return li ,val
if abs(x-x0)<epsilon:return li,vals
x0 = x
def secant(y:sympy.core,x0:float,x1:float,epsilon:float =0.00001,maxtime:int=50) ->(list,list):
'''
弦截法, 使用newton 差商计算,每次只需计算一次f(x)
secant method for finding a zeropoint of a func
y is the func , x0 is the init x val, epsilon is the accurrency
'''
if epsilon <0:epsilon = -epsilon
ct =0
x0,x1 = float(x0),float(x1)
li = [x0,x1]
t = y.free_symbols
varsymbol = 'x' if len(t)==0 else t.pop()
last = y.subs(varsymbol,x0)
vals = [last]
while 1:
cur = y.subs(varsymbol,x1)
vals.append(cur)
x = x1-cur*(x1-x0)/(cur-last)
x0 ,x1= x1,x
last = cur
li.append(x)
ct+=1
if ct>maxtime:
print("after iteration for {} times, I still havn't reach the accurrency.\
Maybe this function havsn't zeropoint\n".format(ct))
return li,vals
if abs(x0-x1)<epsilon:return li,vals
x0 = x
def solveNonlinearEquations(funcs:[sympy.core],init_dic:dict,epsilon:float=0.001,maxtime:int=50)->dict:
'''solve nonlinear equations:'''
li = list(init_dic.keys())
delta = {i:0 for i in li}
ct = 0
while 1:
ys = np.array([f.subs(init_dic) for f in funcs],dtype = 'float')
mat = np.matrix([[i.diff(x).subs(init_dic) for x in li] for i in funcs ],dtype = 'float')
delt = np.linalg.solve(mat,-ys)
for i,j in enumerate(delt):
init_dic[li[i]] +=j
delta[li[i]] = j
if ct>maxtime:
print("after iteration for {} times, I still havn't reach the accurrency.\
Maybe this function havsn't zeropoint\n".format(ct))
return init_dic
if sqrt(sum(i**2 for i in delta.values()))<epsilon:return init_dic
if __name__ =='__main__':
x,y,z = sympy.symbols('x y z')
res,res2= newton(x**5-9,2,0.01)
print(res,res2)
res,res2 = secant (x**3-3*x-2,1,3,1e-3)
print(res,res2)
funcs=[x**2+y**2-1,x**3-y]
init = {x:0.8,y:0.6}
res_dic = solveNonlinearEquations(funcs,init,0.001)
print(res_dic)

View File

@ -0,0 +1,104 @@
''' mbinary
#########################################################################
# File : least_square.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
'''*************************************************************************
> File Name: least-square.py
> Author: mbinary
> Mail: zhuheqin1@gmail.com
> Created Time: Sat 07 Apr 2018 09:55:25 PM DST
> Blog: https://mbinary.github.io
> Description: 最小二乘法解线性方程组, 解矛盾方程组
************************************************************************'''
import re
import numpy as np
def solveConflitEqualtion(A,y):
'''solve equation like this: Av = y,
A:m*n v:n*1 y:m*1
return vector v
'''
A = np.matrix(A)
y = np.array(y)
ata = A.T*A
print('AtA')
print(ata)
return np.linalg.solve(ata,A.T*y) # note that is numpy.linalg.solve
def solveLinear(point,index):
y = [[i[1]] for i in point]
x = [[i[0]] for i in point]
A = []
for i in x:
A.append([i[0]**j for j in index])
res = solveConflitEqualtion(A,y)
print('the solution is : \n',res)
print('namely: ')
items = ['{:.4f}x^{}'.format(res[i,0],j) for i, j in enumerate(index)]
print('phi(x) = ',' + '.join(items))
def handleInput(s=None,y=None):
# numPt = re.compile (r'\d*\.{0,1}\d+')
if not s: s = input('input matrix A:m*n //m>=n\n')
s = s.replace(' ','')
li = re.findall(r'(\[(\d+)(,(\d+))+\])',s)
li = [parseLst(i[0]) for i in li]
if not y:y = input('input a vector y:n*1\n')
y = parseLst(y)
print('Equation: Av = y:')
print('y is as follows: ')
print(y)
print('A is as follows: ')
for i in li:
for j in i:
print('{}'.format(j).rjust(5),end='')
print('')
print('result v is as follows: ')
res = solveConflitEqualtion(li,y)
print(res)
def parseLst(s):
s = s.strip('[]')
li = s.split(',')
li = [float(j) for j in li]
return li
if __name__ == '__main__':
'''
li = '[[23,2],[2,5],[2,6]]'
y = '[1,3]'
while True:
handleInput(li,y)
s = input('input y to continue, n for exit')
if s!='y':break
'''
point = [(-3,14.3),(-2,8.3),(-1,4.7),(2,-8.3),(4,-22.7)]
lst = [0,3]
solveLinear(point,lst)
point= [(-3,14.3),(-2,8.3),(-1,4.7),(2,8.3),(4,22.7)]
lst = [0,2]
solveLinear(point,lst)
A = [[1,2],[2,1],[1,1]]
y = [[5],[6],[4]]
res = solveConflitEqualtion(A,y)
print(res)
A = [[1,-2],[1,5],[2,1],[1,1]]
y = [[1],[13.1],[7.9],[5.1]]
print(solveConflitEqualtion(A,y))

View File

@ -0,0 +1,118 @@
''' mbinary
#########################################################################
# File : linear_equation.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
#coding: utf-8
'''************************************************************************
> File Name: doolittle.py
> Author: mbinary
> Mail: zhuheqin1@gmail.com
> Blog: https://mbinary.github.io
> Created Time: 2018-04-20 08:32
************************************************************************'''
import numpy as np
def getLU(A):
'''doolittle : A = LU,
L is in down-triangle form,
U is in up-triangle form
'''
m,n = A.shape
if m!=n:raise Exception("this matrix is not inversable")
L = np.zeros([m,m])
U = np.zeros([m,m])
L = np.matrix(L)
U = np. matrix(U)
U[0] = A[0]
L[:,0] = A[:,0] / A[0,0]
for i in range(1,m):
for j in range(i,m):
U[i,j]= A[i,j] - sum(L[i,k]*U[k,j] for k in range(i))
L[j,i] = (A[j,i] - sum(L[j,k]*U[k,i] for k in range(i)))/U[i,i]
print(L)
print(U)
return L,U
def gauss_prior_elimination(A):
'''using guass elimination,get up_trianglge form of A'''
m,n = A.shape
if m!=n:raise Exception("[Error]: matrix is not inversable")
B = np.matrix(A,dtype=float) # necessary,otherwise when the dtype of A is int, then it will be wrong
for i in range(m-1):
col = abs(B[i:,i]) # note using abs value, return a matrix in (m-i)x1 form
mx = col.max()
if mx==0: raise Exception("[Error]: matrix is not inversable")
pos = i+col.argmax()
if pos != i : B[[pos,i],:] = B[[i,pos],:] # note how to swap cols/rows
B[i,:] = 1/mx*B[i,:]
for j in range(i+1,m):
#print(B)
B[j,:] -= B[j,i] * B[i,:]
print(B)
return B
def solveDown(A,b):
'''A is a matrix in down-triangle form'''
sol = np.zeros(b.shape)
for i in range(b.shape[0]):
sol[i,0] = (b[i,0]-sum(A[i,j]*sol[j,0] for j in range(i)))/A[i,i]
return sol
def solveUp(A,b):
'''A is a matrix in up-triangle form'''
sol = np.zeros(b.shape)
n = b.shape[0]
for i in range(n-1,-1,-1):
sol[i,0] = (b[i,0]-sum(A[i,j]*sol[j,0] for j in range(n-1,i,-1)))/A[i,i]
return sol
def doolittle(A,b):
L,U = getLU(A)
y = solveDown(L,b)
x = solveUp(U,y)
print(y)
print(x)
return x
def ldlt(A,b):
L,U = getLU(A)
D = np.diag(np.diag(U))
print(D,"D")
z = np.linalg.solve(L,b)
print(z,"z")
y = np.linalg.solve(D,z)
print(y,"y")
x = np.linalg.solve(L.T,y)
print(x,"x")
return x
if __name__ == '__main__':
A = np.matrix([[10,5,0,0],
[2,2,1,0],
[0,10,0,5],
[0,0,2,1]])
b = np.matrix([[5],[3],[27],[6]])
gauss_prior_elimination(A)
'''ldlt
A = np.matrix([[-6,3,2],
[3,5,1],
[2,1,6]])
b = np.matrix([[-4],[11],[-8]])
ldlt(A,b)
'''
'''
A = np.matrix([[2,1,1],
[1,3,2],
[1,2,2]])
b = np.matrix([[4],[6],[5]])
doolittle(A,b)
'''

View File

@ -0,0 +1,13 @@
''' mbinary
#########################################################################
# File : numerical_differential.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''

View File

@ -0,0 +1,70 @@
''' mbinary
#########################################################################
# File : numerical_integration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
#########################################################################
# File : numerical integration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.github.io
# Github: https://github.com/mbinary
# Created Time: 2018-05-11 08:58
# Description:
# numerical intergration: using Newton-Cotes integration, and Simpson
# 数值积分, 使用 牛顿-科特斯积分, 辛普森
#########################################################################
import numpy as np
def trapezoidal(a,b,h,fs):
'''梯形积分公式'''
xs = [i for i in np.arange(a,b+h,h)]
print(xs)
ret = h*(sum(fs)-fs[0]/2 - fs[-1]/2)
print(ret)
return ret
def simpson(a,b,h,fs):
'''辛普森积分公式'''
xs = [i for i in np.arange(a,b+h,h)]
print(xs)
ret = h/3*(4* sum(fs[1::2])+ 2*sum(fs[2:-1:2]) + fs[0]+fs[-1])
print(ret)
return ret
def romberg(a,b,f,epcilon):
'''romberg(龙贝格) 数值积分'''
h = b-a
lst1=[h*(f(a)+f(b))/2]
print(lst1)
delta = epcilon
k=1
while delta >= epcilon:
h/=2
k+=1
lst2=[]
lst2.append((lst1[0]+h*2*sum(f(a+(2*i-1)*h) for i in range(1,2**(k-2)+1)))/2)
for j in range(0,k-1):
lst2.append(lst2[j]+(lst2[j]-lst1[j])/(4**(j+1)-1))
delta = abs(lst2[-1]-lst1[-1])
lst1=lst2
print(lst1)
if __name__=='__main__':
a,b,h = 0.6,1.8,0.2
fs=[5.7,4.6,3.5,3.7,4.9,5.2,5.5]
trapezoidal(a,b,h,fs)
simpson(a,b,h,fs)
romberg(1,2,lambda x:sin(x**4),1e-4)

View File

@ -0,0 +1,130 @@
''' mbinary
#########################################################################
# File : solve-linear-by-iteration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
'''
#########################################################################
# File : solve-linear-by-iteration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.github.io
# Github: https://github.com/mbinary
# Created Time: 2018-05-04 07:42
# Description:
#########################################################################
'''
import numpy as np
from operator import le,lt
def jacob(A,b,x,accuracy=None,times=6):
''' Ax=b, arg x is the init val, times is the time of iterating'''
A,b,x = np.matrix(A),np.matrix(b),np.matrix(x)
n,m = A.shape
if n!=m:raise Exception("Not square matrix: {A}".format(A=A))
if b.shape !=( n,1) : raise Exception('Error: {b} must be {n} x1 in dimension'.format(b = b,n=n))
D = np.diag(np.diag(A))
DI = np.zeros([n,n])
for i in range(n):DI[i,i]= 1/D[i,i]
R = np.eye(n) - DI * A
g = DI * b
print('R =\n{}'.format(R))
print('g =\n{}'.format(g))
last = -x
if accuracy != None:
ct=0
while 1:
ct+=1
tmp = x-last
last = x
mx = max ( abs(i) for i in tmp)
if mx<accuracy:return x
x = R*x+g
print('x{ct} =\n{x}'.format(ct = ct,x=x))
else:
for i in range(times):
x = R*x+g
print('x{ct} = \n{x}'.format(ct=i+1,x=x))
print('isLimitd: {}'.format(isLimited(A)))
return x
def gauss_seidel(A,b,x,accuracy=None,times=6):
''' Ax=b, arg x is the init val, times is the time of iterating'''
A,b,x = np.matrix(A),np.matrix(b),np.matrix(x)
n,m = A.shape
if n!=m:raise Exception("Not square matrix: {A}".format(A=A))
if b.shape !=( n,1) : raise Exception('Error: {b} must be {n} x1 in dimension'.format(b = b,n=n))
D =np. matrix(np.diag(np.diag(A)))
L = np.tril(A) - D # L = np.triu(D.T) - D
U = np.triu(A) - D
DLI = (D+L).I
S = - (DLI) * U
f = (DLI)*b
print('S =\n{}'.format(S))
print('f =\n{}'.format(f))
last = -x
if accuracy != None:
ct=0
while 1:
ct+=1
tmp = x-last
last = x
mx = max ( abs(i) for i in tmp)
if mx<accuracy:return x
x = S*x+f
print('x{ct} =\n{x}'.format(ct=ct,x=x))
else:
for i in range(times):
x = S*x+f
print('x{ct} = \n{x}'.format(ct=i+1,x=x))
print('isLimitd: {}'.format(isLimited(A)))
return x
def isLimited(A,strict=False):
'''通过检查A是否是[严格]对角优来判断迭代是否收敛, 即对角线上的值是否都大于对应行(或者列)的值'''
diag = np.diag(A)
op = lt if strict else le
if op(A.max(axis=0),diag).all(): return True
if op(A.max(axis=1), diag).all(): return True
return False
testcase=[]
def test():
for func,A,b,x,*args in testcase:
acc =None
times = 6
if args !=[] :
if isinstance(args[0],int):times = args[0]
else : acc = args[0]
return func(A,b,x,acc,times)
if __name__ =='__main__':
A = [[2,-1,-1],
[1,5,-1],
[1,1,10]
]
b = [[-5],[8],[11]]
x = [[1],[1],[1]]
#testcase.append([gauss_seidel,A,b,x])
A = [[2,-1,1],[3,3,9],[3,3,5]]
b = [[-1],[0],[4]]
x = [[0],[0],[0]]
#testcase.append([jacob,A,b,x])
A = [[5,-1,-1],
[3,6,2],
[1,-1,2]
]
b= [[16],[11],[-2]]
x = [[1],[1],[-1]]
testcase.append([gauss_seidel,A,b,x,0.001])
test()

View File

@ -0,0 +1,163 @@
''' mbinary
#########################################################################
# File : tongyu_equation.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
# created by mbinary @2018-3-4
# description: solve tongyu equation
# notice that i use -- to repr tongyu symbol
from isPrime import isPrime,primes
from operator import mul, and_
from functools import reduce,partial
import re
def primeFactoize(x):
'''质因数分解 , ret {p:r}'''
if isPrime(x):return {x:1}
mp={}
for i in primes:
if x==1:break
ct=0
while x%i==0:
ct+=1
x//=i
if ct!=0:mp[i]=ct
return mp
def xgcd(a,b):
'''ax+by=gcd(a,b) ,用辗转相除法得到gcd,x,y'''
def _xgcd(a,b):
if b==0:return a,1,0
gcd,x,y=_xgcd(b,a%b)
return gcd,y,x-y*a//b
if a<b:
g,x,y = _xgcd(b,a)
return g,y,x
return _xgcd(a,b)
def gcd(a,b):
return a if b==0 else gcd(b,a%b)
def lcm(a,b):
return a*b//gcd(a,b)
def euler(x):
mp = primeFactoize(x)
fac = [1-1/i for i in mp]
return round(reduce(mul,fac,x))
def ind(m,g):
''' mod m ,primary root g -> {n:indg n}'''
return {j:i for i in range(m-1) \
for j in range(m) if (g**i-j)%m==0}
def gs(m,num=100):
'''return list of m's primary roots below num'''
phi = euler(m)
mp = primeFactoize(phi)
checkLst = [phi//i for i in mp]
return [i for i in range(2,num) \
if reduce(and_,[(i**n-1)%m !=0 for n in checkLst])]
def minG(m):
phi = euler(m)
mp = primeFactoize(phi)
checkLst = [phi//i for i in mp]
i=2
while 1:
if reduce(and_,[(i**n-1)%m !=0 for n in checkLst]):return i
i+=1
class solve:
def __init__(self,equ=None):
self.linearPat= re.compile(r'\s*(\d+)\s*--\s*(\d+)[\s\(]*mod\s*(\d+)')
self.sol = []
#self.m = m
#self.ind_mp = ind(m,minG(m))
def noSol(self):
print('equation {equ} has no solution'.format(equ=self.equ))
def error(self):
print("Error! The divisor m must be postive integer")
def solveLinear(self,a,b,m):
'''ax--b(mod m): solve linear equation with one unknown
return ([x1,x2,...],m)
'''
a,b,m = self.check(a,b,m)
g,x,y=xgcd(a,m)
if a*b%g!=0:
self.noSol()
return None
sol=x*b//g
m0 = m//g
return ([(sol+i*m0)%m for i in range(g)],m)
def check(self,a,b,m):
if m<=0:
self.error()
return None
if a<0:a,b=-a,-b ## important
if b<0:b+= -b//m * m
return a,b,m
#def solvePoly(self,m,mp):
''' mod m, mp:{index:coef} is a dict of the polynomials' coefficient and index'''
''' g = minG(m)
ind_mp = ind(m,g)
li = []
for i in mp:
solve
'''
def solveHigh(self,a,n,b,m):
''' ax^n -- b (mod m) ind_mp is a dict of m's {n: indg n}'''
ind_mp = ind(m,minG(m))
tmp = ind_mp[b] - ind_mp[a]
if tmp < 0:tmp+=m
print(n,tmp)
sol = self.solveLinear(n,tmp,euler(m))
re_mp = {j:i for i ,j in ind_mp.items()}
print(sol)
return [re_mp[i] for i in sol[0]],m
def solveGroup(tups):
'''tups is a list of tongyu equation groups, like
[(a1,b1,m1),(a2,b2,m2)...]
and, m1,m2... are all primes
'''
mp = {}
for a,b,m in tups:
if m in mp:
if mp[m][0]*b!=mp[m][1]*a:
self.noSol()
return
else:mp[m] = (a,b)
product= reduce(lambda i,j:i*j, mp.keys(), 1)
sol = 0
for i in mp:
x = self.solveLinear(product//i*mp[i][0],1,i)
sol+= x*product//i*mp[i][1]
sol%=m
return ([sol],m)
def __call__(self):
s=input('输入同余方程,用--代表同于号形如3--5(mod 7)代表3x模7同余于5')
li= self.linearPat.findall(s)
li = [(int(a),int(b),int(m)) for a,b,m in li]
print(self.solveLinear(li[0]))
if __name__ == '__main__':
solver = solve()
res = solver.solveLinear(3,6,9)
print(res)
print(solver.solveHigh(1,8,3,11))

View File

@ -0,0 +1,134 @@
''' mbinary
#########################################################################
# File : vector_norm.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-10-02 21:14
# Description:
#########################################################################
'''
from random import randint,random
import numpy as np
from operator import neg,and_
from functools import reduce
class obj():
def __init__(self,data):
self.data=np.array(data)
def __add__(self,x):
data = x.data if self.__class__ == x.__class__ else x
return self.__class__(self.data + data)
def __radd__(self,x):
data = x.data if self.__class__ == x.__class__ else x
return self.__class__(data +self.data)
def __iadd__(self,x):
data = x.data if self.__class__ == x.__class__ else x
self.data += data
def __mul__(self,x):
data = x.data if self.__class__ == x.__class__ else x
return self.__class__(self.data * data)
def __imul__(self,x):
data = x.data if self.__class__ == x.__class__ else x
self.data *= data
def __rmul__(self,x):
data = x.data if self.__class__ == x.__class__ else x
return self.__class__(data * self.data)
def __neg__(self):
return neg(self)
def __abs__(self):
return abs(self.data)
'''
@property
def data(self):
return self.data
@data.setter
def data(self,s):
self.data = s
'''
def norm(self,n=0):
'''the default is +oo norm'''
absolute = abs(self.data)
if n < 1 :return max(absolute)
return (sum(absolute**n))**(1/n)
def hasNorm(self):
'''check norm's three necessary conditions:
1. not neg
2. homogenious (qici)
3. triangle inequlity
there is much probably wrong
'''
bl = reduce(and_,[self.norm(i)>=0 for i in range(3)])
if bl:
n = randint(2,100)
bl = reduce(and_,[n*(self.norm(i))==(n*self).norm(i) for i in range(3)])
if bl:
another = self*randint(2,10)-randint(1,100)
return reduce(and_,[(another+self).norm(i)<=another.norm(i)+self.norm(i) for i in range(3)])
return False
class vector(obj):
def __init__(self,arr):
''' arr: iterable'''
self.data =np.array(arr)
def innerProduct(self,x):
return sum(self.data*x)
def outerProduct(self,x):
pass
class matrix(obj):
def __init__(self,s):
'''s is a list of lists'''
self.data=np.mat(s)
self.T = None
self. I = None
'''
@property
def T(self):
if self.T==None:self.T = self.data.T
return self.T
@T.setter
def T(self,s):
self.T = s
@property
def I(self):
if self.I == None: self.I = self.data.I
return self.I
@I.setter
def I(self,s):
self.I = s
'''
def E(self,n=None):
if n is None: n = self.data.shape[0]
return np.eye(n)
def norm(self,n=0):
absolute = abs(self.data)
if n < 1:
# max of one row sum
return max([sum(i) for i in absolute])
if n==1:return self.norm1()
elif n==2:return self.norm2()
def norm1(self):
''' max of sum of cols'''
absolute = abs(self.data)
return max(absolute.sum(axis=0))
def norm2(self):
''' max of sum of rows'''
absolute = abs(self.data)
return max(absolute.sum(axis=1))
def norm_f(self):
return sum((self.data**2).sum(axis=1))**0.5
if __name__ =='__main__':
v1 = vector([1,-2,3,4])
v2 = vector([0,2,0,5])
m1 = matrix([v1,v2,v2,v1])
print([v1.norm(i) for i in range(3)])

View File

@ -1,3 +1,15 @@
''' mbinary
#########################################################################
# File : bTree.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-08-29 12:49
# Description:
#########################################################################
'''
class node:
def __init__(self,keys=None,isLeaf = True,children=None):
if keys is None:keys=[]

View File

@ -1,3 +1,15 @@
''' mbinary
#########################################################################
# File : redBlackTree.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-07-14 16:15
# Description:
#########################################################################
'''
'''
#########################################################################
# File : redBlackTree.py

View File

@ -1,3 +1,29 @@
---
title: 『数据结构』Fibonacci-heap
date: 2018-09-06 19:09
categories: 数据结构与算法
tags: [数据结构,斐波那契堆]
keywords: 数据结构,斐波那契堆
mathjax: true
description:
---
<!-- TOC -->
- [1. 结构](#1-结构)
- [2. 势函数](#2-势函数)
- [3. 最大度数](#3-最大度数)
- [4. 操作](#4-操作)
- [4.1. 创建一个斐波那契堆](#41-创建一个斐波那契堆)
- [4.2. 插入一个结点](#42-插入一个结点)
- [4.3. 寻找最小结点](#43-寻找最小结点)
- [4.4. 合并两个斐波那契堆](#44-合并两个斐波那契堆)
- [4.5. 抽取最小值](#45-抽取最小值)
- [4.6. 关键字减值](#46-关键字减值)
- [4.7. 删除结点](#47-删除结点)
- [5. 最大度数的证明](#5-最大度数的证明)
<!-- /TOC -->
![](https://upload-images.jianshu.io/upload_images/7130568-22531846a72b0d83.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
@ -12,16 +38,15 @@
<a id="markdown-2-势函数" name="2-势函数"></a>
# 2. 势函数
下面用势函数来分析摊还代价, 如果你不明白, 可以看[摊还分析](https://www.jianshu.com/p/052fbe9d92a4)
$\Phi(H) = t(H) + 2m(h)$
t 是根链表中树的数目,m(H) 表示被标记的结点数
最初没有结点
<a id="markdown-3-最大度数" name="3-最大度数"></a>
# 3. 最大度数
$D(n)\leqslant \lfloor lgn \rfloor$
![](https://upload-images.jianshu.io/upload_images/7130568-c9e0cd3be4e98c4b.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
结点的最大度数(即孩子数)$D(n)\leqslant \lfloor lgn \rfloor$, 证明放在最后
<a id="markdown-4-操作" name="4-操作"></a>
# 4. 操作
<a id="markdown-41-创建一个斐波那契堆" name="41-创建一个斐波那契堆"></a>
@ -133,19 +158,13 @@ def cascading-cut(H,y):
<a id="markdown-47-删除结点" name="47-删除结点"></a>
## 4.7. 删除结点
```python
decrease(H,nd, MIN)
<!-- TOC -->
extract-min(H)
```
- [1. 结构](#1-结构)
- [2. 势函数](#2-势函数)
- [3. 最大度数](#3-最大度数)
- [4. 操作](#4-操作)
- [4.1. 创建一个斐波那契堆](#41-创建一个斐波那契堆)
- [4.2. 插入一个结点](#42-插入一个结点)
- [4.3. 寻找最小结点](#43-寻找最小结点)
- [4.4. 合并两个斐波那契堆](#44-合并两个斐波那契堆)
- [4.5. 抽取最小值](#45-抽取最小值)
- [4.6. 关键字减值](#46-关键字减值)
- [4.7. 删除结点](#47-删除结点)
<!-- /TOC -->
<a id="markdown-5-最大度数的证明" name="5-最大度数的证明"></a>
# 5. 最大度数的证明
这也是`斐波那契`这个名字的由来,
$D(n)\leqslant \lfloor lgn \rfloor$
![](https://upload-images.jianshu.io/upload_images/7130568-c9e0cd3be4e98c4b.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

View File

@ -1,3 +1,13 @@
---
title: 图算法
date: 2018-09-06 19:10
categories: 数据结构与算法
tags: [图,算法]
keywords: 图,算法
mathjax: true
description:
---
<!-- TOC -->
- [1. 图](#1-图)
@ -161,7 +171,7 @@ for edge as u,v in sorted(G.E):
return edges
```
如果并查集的实现采用了 按秩合并与路径压缩技巧, 则 find 与 union 的时间接近常数
所以时间复杂度在于排序边, 即 $O(ElgE)$, 而 $E<V^2$, $lgE = O(lgV)$, $O(ElgV)$
所以时间复杂度在于排序边, 即 $O(ElgE)$, 而 $ E\< V^2 $, 所以 $lgE = O(lgV)$, 时间复杂度为 $O(ElgV)$
<a id="markdown-32-prim-算法" name="32-prim-算法"></a>
## 3.2. Prim 算法
用了 BFS, 类似 Dijkstra 算法
@ -190,7 +200,7 @@ while not que.isempty():
<a id="markdown-4-单源最短路" name="4-单源最短路"></a>
# 4. 单源最短路
求一个结点到其他结点的最短路径, 可以用 Bellman-ford算法, 或者 Dijkstra算法.
求一个结点到其他结点的最短路径, 可以用 Bellman-Ford算法, 或者 Dijkstra算法.
定义两个结点u,v间的最短路
$$
\delta(u,v) = \begin{cases}
@ -205,7 +215,7 @@ $$
<a id="markdown-41-负权重的边" name="41-负权重的边"></a>
## 4.1. 负权重的边
Dijkstra 算法不能处理, 只能用 Bellman-Ford 算法,
Dijkstra 算法不能处理负权边, 只能用 Bellman-Ford 算法,
而且如果有负值圈, 则没有最短路, bellman-ford算法也可以检测出来
<a id="markdown-42-初始化" name="42-初始化"></a>
## 4.2. 初始化
@ -288,7 +298,7 @@ $$
$$
\delta(i,j) = l_{ij}^{(|V|-1)} = l_{ij}^{(|V|)} =l_{ij}^{(|V| + 1)}= ...
$$
所以客户处自底向上计算, 如下
所以可以自底向上计算, 如下
输入权值矩阵 $W(w_{ij})), L^{(m-1)}$,输出$ L^{(m)}$, 其中 $L^{(1)} = W$,
```python
n = L.rows
@ -461,4 +471,4 @@ def ford-fulkerson(G,s,t):
<a id="markdown-7-参考资料" name="7-参考资料"></a>
# 7. 参考资料
[^1]: 算法导论
[^2]: 图论, 王树禾
[^2]: 图论, 王树禾

View File

Before

Width:  |  Height:  |  Size: 117 KiB

After

Width:  |  Height:  |  Size: 117 KiB

View File

Before

Width:  |  Height:  |  Size: 38 KiB

After

Width:  |  Height:  |  Size: 38 KiB

View File

Before

Width:  |  Height:  |  Size: 104 KiB

After

Width:  |  Height:  |  Size: 104 KiB

View File

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 22 KiB

View File

Before

Width:  |  Height:  |  Size: 29 KiB

After

Width:  |  Height:  |  Size: 29 KiB

View File

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 16 KiB

View File

Before

Width:  |  Height:  |  Size: 23 KiB

After

Width:  |  Height:  |  Size: 23 KiB

View File

@ -1,3 +1,15 @@
''' mbinary
#########################################################################
# File : lcs.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-08-25 12:00
# Description:
#########################################################################
'''
def lcs(a,b):
'''time: O(mn); space: O(mn)'''
m,n= len(a),len(b)

View File

@ -0,0 +1,16 @@
''' mbinary
#########################################################################
# File : matrix-multiply.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-08-24 21:24
# Description:
#########################################################################
'''
def adjustOrd(sizes):
''' adjust the chain-multiply of matrix, sizes=[row1,row2,..,rown,coln]'''
n = len(sizes)
if n<3: return

View File

@ -1,3 +1,15 @@
''' mbinary
#########################################################################
# File : splitStripe.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-08-24 17:07
# Description:
#########################################################################
'''
'''
There is stripe which length is n,
priceMap contains a map for different length of stripe and its price

44
math/isPrime.py Normal file
View File

@ -0,0 +1,44 @@
''' mbinary
#########################################################################
# File : isPrime.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-05-19 21:34
# Description:
#########################################################################
'''
# created by mbinary @2018-3-4
# description: judge a num if it's a prime. It will be more efficient when judging many times
primes = [2,3,5,7,11,13]
def isPrime(x):
global primes
if x>primes[-1]:
return genPrime(x)
return twoDivideFind(x,primes)
def genPrime(x):
global primes
while x>primes[-1]:
left = primes[-1]
right = (left+1)**2
lst = []
for i in range(left,right):
for j in primes:
if i%j==0:break
else:lst.append(i)
primes+=lst
else:return twoDivideFind(x,lst)
def twoDivideFind(x,primes):
a,b = 0, len(primes)
while a<=b:
mid = (a+b)//2
if primes[mid]<x:a=mid+1
elif primes[mid]>x: b= mid-1
else:return True
return False

50
math/num_weight.py Normal file
View File

@ -0,0 +1,50 @@
''' mbinary
#########################################################################
# File : num_weight.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.coding.me
# Github: https://github.com/mbinary
# Created Time: 2018-05-19 21:36
# Description:
#########################################################################
'''
def covert(s,basefrom=10,baseto=2):
return d2n(n2d(s,basefrom),baseto)
def n2d(s,base=16):
''' num of base_n(n<36) to decimal'''
dic = {chr(i+ord('0')):i for i in range(10)}
s=s.upper()
if base>10:
dic.update({chr(i+ord('A')):i+10 for i in range(26)})
#if base in [16,8,2] :
# p=max(map(s.find,'OBX'))
# s=s[p+1:] #remove prefix of hex or bin or oct
rst=0
for i in s:
rst=dic[i]+rst*base
return rst
def d2n(n,base=16):
''' num of base_n(n<36) to decimal'''
dic = {i:chr(i+ord('0')) for i in range(10)}
if base>10:
dic.update({i+10:chr(i+ord('A')) for i in range(26)})
rst=[]
while n!=0:
i=int(n/base)
rst.append(dic[n-i*base])
n=i
return ''.join(rst[::-1])
'''
>>> n2d(str(d2n(4001)))
4001
>>> d2n(n2d(str(4001)),2)
'100000000000001'
>>> covert('4001',16,2)
'100000000000001'
'''