Xenharmonic Wiki

Mathematics of MOS: Difference between revisions

← Older editNewer edit →
VisualWikitext
Wikispaces>FREEZE
No edit summary
Spt3125 (talk | contribs)
autopatrolled, Sysops
523 edits
→Algorithms: reformatted code (damaged by MediaWiki import)
Line 72: Line 72:


=Algorithms=
=Algorithms=
Below is some Maple code for various mathematical routines having to do with MOS. If you have access to Maple, you can of course copy and run these programs. Even if you do not, since Maple code makes better pseudocode than most languages or computer algebra packages afford, it can be used as pseudocode. For that purpose, it will be helpful to know that "modp(x, n)" means reducing x mod the integer n to 0, 1, ..., n-1 not only when x is an integer, but also when it is a rational number with denominator prime to n. In that case, p/q mod n = r means p = qr mod n.
Below is some Maple code for various mathematical routines having to do with MOS. If you have access to Maple, you can of course copy and run these programs. Even if you do not, since Maple code makes better pseudocode than most languages or computer algebra packages afford, it can be used as pseudocode. For that purpose, it will be helpful to know that <code>modp(x, n)</code> means reducing x mod the integer n to 0, 1, ..., n-1 not only when x is an integer, but also when it is a rational number with denominator prime to n. In that case, p/q mod n = r means p = qr mod n.


log2 := proc(x)
Note that some procedures below depend on others.  Special procedure names (i.e. not built into Maple) are shown in <code>'''bold type'''</code>.


<ol><li>logarithm base 2</li></ol>evalf(ln(x)/ln(2)) end:
'''log2''' := proc(x)
# logarithm base 2
evalf(ln(x)/ln(2)) end:


nextfarey := proc(q, n)
'''nextfarey''' := proc(q, n)
# next in row n of Farey sequence from q, 0 <= q <= 1
local a, b, r, s;
if q >= (n-1)/n then RETURN(1) fi;
a := numer(q);
b := denom(q);
s := n - modp(n + 1/a, b);
r := modp(1/b, s);
r/s end:


<ol><li>next in row n of Farey sequence from q, 0 &lt;= q &lt;= 1</li></ol>local a, b, r, s;
'''prevfarey''' := proc(q, n)
# previous in row n of Farey sequence from q, 0 <= q <= 1
local a, b, r, s;
if q=0 then RETURN(0) fi;
if n=0 then RETURN(0) fi;
a := numer(q);
b := denom(q);
s := n - modp(n - 1/a, b);
r := modp(-1/b, s);
r/s end:


if q &gt;= (n-1)/n then RETURN(1) fi;
'''fareypair''' := proc(q)
# Farey pair with q as its mediant
local n;
n := denom(q);
['''prevfarey'''(q, n), '''nextfarey'''(q, n)] end:


a := numer(q);
'''mediant''' := proc(u, v)
# mediant of two rational numbers u and v
(numer(u) + numer(v))/(denom(u) + denom(v)) end:


b := denom(q);
'''convergents''' := proc(z)
# convergent list for z
local q;
convert(z,confrac,'q');
q end:


s := n - modp(n + 1/a, b);
'''exlist''' := proc(l)
# expansion of a convergent list to semiconvergents
local i, j, s, d;
if nops(l)<3 then RETURN(l) fi;
d[1] := l[1];
d[2] := l[2];
s := 3;
for i from 3 to nops(l)-1 do
for j from 1 to (numer(l[i])-numer(l[i-2]))/numer(l[i-1]) do
d[s] :=
(j*numer(l[i-1])+numer(l[i-2]))/(j*denom(l[i-1])+denom(l[i-2]));
s := s+1 od od;
convert(convert(d, array), list) end:


r := modp(1/b, s);
'''semiconvergents''' := proc(z)
# semiconvergent list for z
'''exlist'''('''convergents'''(z)) end:


r/s end:
'''penult''' := proc(q)
# penultimate convergent to q
local i, u;
u := '''convergents'''(q);
if nops(u)=1 then RETURN(u[1]) fi;
for i from 1 to nops(u) do
if u[i]=q then RETURN(u[i-1]) fi od;
end:


prevfarey := proc(q, n)
'''Ls''' := proc(q)
# large-small steps from mediant q
local u;
u := '''fareypair'''(q);
[denom(u[2]), denom(u[1])] end:


<ol><li>previous in row n of Farey sequence from q, 0 &lt;= q &lt;= 1</li></ol>local a, b, r, s;
'''medi''' := proc(u)
# mediant from Large-small steps
local q, r;
if u[2]=1 then RETURN(1/(u[1]+1)) fi;
r := igcd(u[1], u[2]);
if r>1 then RETURN(medi([u[1]/r, u[2]/r])) fi;
q := '''penult'''(u[1]/u[2]);
if q > u[1]/u[2] then RETURN((numer(q)+denom(q))/(u[1]+u[2])) fi;
(u[1]+u[2]-numer(q)-denom(q))/(u[1]+u[2]) end:


if q=0 then RETURN(0) fi;
'''Lsgen''' := proc(g, n)
# given generator g and scale size n determines large-small steps
local q, u, w;
q := round(n*g)/n;
w := n/denom(q);
u := '''fareypair'''(q);
if g<u[1] or g>u[2] or g=q then RETURN('false') fi;
if g<q then RETURN([w*denom(u[1]), w*denom(u[2])]) fi;
[w*denom(u[2]), w*denom(u[1])] end:


if n=0 then RETURN(0) fi;
'''revlist''' := proc(l)
# reverse of list
local i, v, e;
e := nops(l);
for i from 1 to e do
v[i] := l[e-i+1] od;
convert(convert(v,array),list) end:


a := numer(q);
'''invcon''' := proc(l)
# inverse continued fraction
local d, i, h, k;
h[-2] := 0;
h[-1] := 1;
k[-2] := 1;
k[-1] := 0;
for i from 0 to nops(l)-1 do
h[i] := l[i+1]*h[i-1] + h[i-2];
k[i] := l[i+1]*k[i-1] + k[i-2];
d[i+1] := h[i]/k[i] od;
convert(convert(d, array), list) end:


b := denom(q);
'''quest''' := proc(x)
# Minkowski ? function
local i, j, d, l, s, t;
l := convert(x, confrac);
d := nops(l);
s := l[1];
for i from 2 to d do
t := 1;
for j from 2 to i do
t := t - l[j] od;
s := s + (-1)^i * 2^t od;
if type(x, float) then s := evalf(s) fi;
s end:


s := n - modp(n - 1/a, b);
'''Box''' := proc(x)
# inverse ? function
local d, e, i, n, w, y;
if type(x, integer) then RETURN(x) fi;
y := x-floor(x);
if y = 1/8 then RETURN(floor(x)+1/4) fi;
w := round('''log2'''(10)*Digits)-5;
n := round(2^w * y);
i :=0;
while n>0 do
i := i+1;
if modp(n,2)=0 then
d[i] := padic[ordp](n, 2);
n := n/2^d[i];
else
d[i] := padic[ordp](n+1, 2);
n := (n-2^d[i]+1)/2^d[i] fi od;
e := convert(convert(d, array), list);
e := subsop(1=NULL,e);
w := ceil(-'''log2'''(y));
e := [op(e), w];
e := [op(e), floor(x)];
e := '''revlist'''(e);
n := '''invcon'''(e);
w := n[nops(n)];
if type(x, rational) and modp(denom(x), 2)=0 then RETURN(w) fi;
evalf(w) end:


r := modp(-1/b, s);
r/s end:
fareypair := proc(q)
<ol><li>Farey pair with q as its mediant</li></ol>local n;
n := denom(q);
[prevfarey(q, n), nextfarey(q, n)] end:
mediant := proc(u, v)
<ol><li>mediant of two rational numbers u and v</li></ol>(numer(u) + numer(v))/(denom(u) + denom(v)) end:
convergents := proc(z)
<ol><li>convergent list for z</li></ol>local q;
convert(z,confrac,'q');
q end:
exlist := proc(l)
<ol><li>expansion of a convergent list to semiconvergents</li></ol>local i, j, s, d;
if nops(l)&lt;3 then RETURN(l) fi;
d[1] := l[1];
d[2] := l[2];
s := 3;
for i from 3 to nops(l)-1 do
for j from 1 to (numer(l[i])-numer(l[i-2]))/numer(l[i-1]) do
d[s] :=
(j*numer(l[i-1])+numer(l[i-2]))/(j*denom(l[i-1])+denom(l[i-2]));
s := s+1 od od;
convert(convert(d, array), list) end:
semiconvergents := proc(z)
<ol><li>semiconvergent list for z</li></ol>exlist(convergents(z)) end:
penult := proc(q)
<ol><li>penultimate convergent to q</li></ol>local i, u;
u := convergents(q);
if nops(u)=1 then RETURN(u[1]) fi;
for i from 1 to nops(u) do
if u[i]=q then RETURN(u[i-1]) fi od;
end:
Ls := proc(q)
<ol><li>large-small steps from mediant q</li></ol>local u;
u := fareypair(q);
[denom(u[2]), denom(u[1])] end:
medi := proc(u)
<ol><li>mediant from Large-small steps</li></ol>local q, r;
if u[2]=1 then RETURN(1/(u[1]+1)) fi;
r := igcd(u[1], u[2]);
if r&gt;1 then RETURN(medi([u[1]/r, u[2]/r])) fi;
q := penult(u[1]/u[2]);
if q &gt; u[1]/u[2] then RETURN((numer(q)+denom(q))/(u[1]+u[2])) fi;
(u[1]+u[2]-numer(q)-denom(q))/(u[1]+u[2]) end:
Lsgen := proc(g, n)
<ol><li>given generator g and scale size n determines large-small steps</li></ol>local q, u, w;
q := round(n*g)/n;
w := n/denom(q);
u := fareypair(q);
if g&lt;u[1] or g&gt;u[2] or g=q then RETURN('false') fi;
if g&lt;q then RETURN([w*denom(u[1]), w*denom(u[2])]) fi;
[w*denom(u[2]), w*denom(u[1])] end:
revlist := proc(l)
<ol><li>reverse of list</li></ol>local i, v, e;
e := nops(l);
for i from 1 to e do
v[i] := l[e-i+1] od;
convert(convert(v,array),list) end:
invcon := proc(l)
<ol><li>inverse continued fraction</li></ol>local d, i, h, k;
h[-2] := 0;
h[-1] := 1;
k[-2] := 1;
k[-1] := 0;
for i from 0 to nops(l)-1 do
h[i] := l[i+1]*h[i-1] + h[i-2];
k[i] := l[i+1]*k[i-1] + k[i-2];
d[i+1] := h[i]/k[i] od;
convert(convert(d, array), list) end:
quest := proc(x)
<ol><li>Minkowski ? function</li></ol>local i, j, d, l, s, t;
l := convert(x, confrac);
d := nops(l);
s := l[1];
for i from 2 to d do
t := 1;
for j from 2 to i do
t := t - l[j] od;
s := s + (-1)^i * 2^t od;
if type(x, float) then s := evalf(s) fi;
s end:
Box := proc(x)
<ol><li>inverse ? function</li></ol>local d, e, i, n, w, y;
if type(x, integer) then RETURN(x) fi;
y := x-floor(x);
if y = 1/8 then RETURN(floor(x)+1/4) fi;
w := round(log2(10)*Digits)-5;
n := round(2^w * y);
i :=0;
while n&gt;0 do
i := i+1;
if modp(n,2)=0 then
d[i] := padic[ordp](n, 2);
n := n/2^d[i];
else
d[i] := padic[ordp](n+1, 2);
n := (n-2^d[i]+1)/2^d[i] fi od;
e := convert(convert(d, array), list);
e := subsop(1=NULL,e);
w := ceil(-log2(y));
e := [op(e), w];
e := [op(e), floor(x)];
e := revlist(e);
n := invcon(e);
w := n[nops(n)];
if type(x, rational) and modp(denom(x), 2)=0 then RETURN(w) fi;
evalf(w) end:
[[Category:math]]
[[Category:math]]
[[Category:mos]]
[[Category:mos]]
[[Category:theory]]
[[Category:theory]]
Retrieved from "https://en.xen.wiki/w/Mathematics_of_MOS"