-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
3e5493c
commit cf50371
Showing
26 changed files
with
1,113 additions
and
585 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,118 @@ | ||
local matrix_mt = {} | ||
matrix_mt.__index = matrix_mt | ||
function matrix(t) | ||
return setmetatable(t,matrix_mt) | ||
end | ||
local function is_matrix(a) | ||
return getmetatable(a)==matrix_mt | ||
end | ||
function matrix_mt:size() | ||
return #self,#self[1] | ||
end | ||
function matrix_mt:new(m,n) -- Define new Matrix | ||
local mt = {} -- Make new array for Matrix | ||
for i=1,m do mt[i] = {} end -- Null matrix elements | ||
return matrix(mt) -- Set Matrix metatable | ||
end | ||
function matrix_mt.__mul(a,b) | ||
if not is_matrix(a) then a,b = b,a end | ||
local m,n = a:size() | ||
local res = a:new(m,n) | ||
if not is_matrix(b) then | ||
for i=1,m do for j=1,n do res[i][j] = b*a[i][j] end end | ||
else | ||
local mb,nb = b:size() | ||
assert(n==mb) | ||
for i=1,m do | ||
for j=1,nb do | ||
local sum = 0 | ||
for k=1,n do sum = sum + a[i][k]*b[k][j] end | ||
res[i][j] = sum | ||
end | ||
end | ||
end | ||
return res | ||
end | ||
function matrix_mt.__tostring(mx) -- Convert to string for printing | ||
local m,n = mx:size() | ||
local s = '{{' | ||
for i=1,m do | ||
for j=1,n-1 do | ||
s = s..(tostring(mx[i][j]) or '..')..', ' end | ||
s = s..(tostring(mx[i][n]) or '..')..'}' | ||
if i==m then s = s..'}' else s = s..',\n{' end | ||
end | ||
return s -- Printable string | ||
end | ||
--[[ | ||
require"num.complex" | ||
require"num.filter_design" | ||
function MM(k) | ||
return matrix{{Zpow(1),k*Zpow(-1)},{k*Zpow(1),Zpow(-1)}}*(1/(k+1)) | ||
end | ||
function KLJunctionGain(kaes,rG,rL) | ||
kaes = TA(kaes):reverse() | ||
local res = matrix{{1},{rL}} | ||
for i,k in ipairs(kaes) do | ||
res = MM(k)*res | ||
end | ||
res = matrix{{Zpow(1),-rG*Zpow(-1)}}*res | ||
return 1/res[1][1] | ||
end | ||
function KLJunctionGain2(kaes,rG,rL) | ||
local res = matrix{{Zpow(1),-rG*Zpow(-1)}}--matrix{{1},{rL}} | ||
for i,k in ipairs(kaes) do | ||
res = res*MM(k) | ||
end | ||
res = res*matrix{{1},{rL}} | ||
return 1/res[1][1] | ||
end | ||
require"sc.utils" | ||
require"sc.utilsstream" | ||
Tract = require"num.Tract"(18,nil,nil,false) | ||
function Areas2KK(AA) | ||
local Kaes = {} | ||
for i=2,#AA do | ||
local num = (AA[i-1]-AA[i]) | ||
local den = (AA[i]+AA[i-1]) | ||
if AA[i]==math.huge then | ||
Kaes[i-1] = 1 | ||
else | ||
if num == 0 or den == 0 then | ||
Kaes[i-1] = 0 | ||
else | ||
Kaes[i-1] = num/den | ||
end | ||
end | ||
end | ||
return Kaes | ||
end | ||
areas = Tract.areas.E | ||
kaes = Areas2KK(areas) | ||
print(KLJunctionGain({1,2,3},0.95,0.97)) | ||
print(KLJunctionGain2({1,2,3},0.95,0.97)) | ||
fil = KLJunctionGain(kaes,0.95,0.97) | ||
print(MM(3)) | ||
pp = FilterPoly({1,1},{1})--*Zpow(-1) | ||
print(pp) | ||
pp2 = RatPoly({1,1},{1})--*Zpow(-1) | ||
print(pp2) | ||
aa = matrix{{1,1},{0,1}} | ||
vec = matrix{{2,3}} | ||
vec2 = matrix{{2},{3}} | ||
print(aa:size()) | ||
print(aa*vec2) | ||
--]] | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,78 @@ | ||
--[[/* | ||
** find rational approximation to given real number | ||
** David Eppstein / UC Irvine / 8 Aug 1993 | ||
** | ||
** With corrections from Arno Formella, May 2008 | ||
** | ||
** usage: a.out r d | ||
** r is real number to approx | ||
** d is the maximum denominator allowed | ||
** | ||
** based on the theory of continued fractions | ||
** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...))) | ||
** then best approximation is found by truncating this series | ||
** (with some adjustments in the last term). | ||
** | ||
** Note the fraction can be recovered as the first column of the matrix | ||
** ( a1 1 ) ( a2 1 ) ( a3 1 ) ... | ||
** ( 1 0 ) ( 1 0 ) ( 1 0 ) | ||
** Instead of keeping the sequence of continued fraction terms, | ||
** we just keep the last partial product of these matrices. | ||
*/ | ||
--]] | ||
function GCD(a,b) | ||
local r | ||
if a < b then b,a = a,b end | ||
while true do | ||
r = a%b | ||
if r == 0 then return b end | ||
--print(a,b,r) | ||
a,b = b,r | ||
end | ||
end | ||
--always gives reduced fractions, dont need GCD | ||
function torational(x,maxden) | ||
local startx = x | ||
-- initialize matrix */ | ||
local m00, m11 = 1,1; | ||
local m01, m10 = 0,0; | ||
|
||
-- loop finding terms until denom gets too big */ | ||
local t | ||
local ai = math.floor(x) | ||
while (m10 * ai + m11 <= maxden) do | ||
print(ai) | ||
t = m00 * ai + m01; | ||
m01 = m00; | ||
m00 = t; | ||
t = m10 * ai + m11; | ||
m11 = m10; | ||
m10 = t; | ||
if(x == ai) then break end -- AF: division by zero | ||
x = 1/(x - ai); | ||
if(x >0x7FFFFFFF) then break end -- AF: representation failure | ||
ai = math.floor(x) | ||
end | ||
|
||
-- now remaining x is between 0 and 1/ai */ | ||
-- approx as either 0 or 1/m where m is max that will fit in maxden */ | ||
-- first try zero */ | ||
-- print( string.format("%d/%d, error = %e\n", m[0][0], m[1][0],startx - ( m[0][0] / m[1][0]))); | ||
return m00,m10,startx - ( m00 / m10) | ||
|
||
-- now try other possibility */ | ||
--ai = (maxden - m[1][1]) / m[1][0]; | ||
--m[0][0] = m[0][0] * ai + m[0][1]; | ||
--m[1][0] = m[1][0] * ai + m[1][1]; | ||
--print(string.format("%d/%d, error = %e\n", m[0][0], m[1][0], startx - ( m[0][0] / m[1][0]))); | ||
end | ||
|
||
--print(torational(2*5*7/(11*3*13),10000)) | ||
--print(torational((3/2)^12,1e3)) | ||
--print(GCD(1e17*13,1e17*11)) | ||
--[[ | ||
local N,M = torational(math.pi,1000) | ||
local gcd = MCD(N,M) | ||
print(N,M) | ||
print(N/gcd,M/gcd) | ||
--]] |
Oops, something went wrong.