Mathematics of MOS: Difference between revisions

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

# logarithm base 2

evalf(ln(x)/ln(2)) end:

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

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

'''fareypair''' := proc(q)

# Farey pair with q as its mediant

local n;

n := denom(q);

['''prevfarey'''(q, n), '''nextfarey'''(q, n)] end:

'''mediant''' := proc(u, v)

# mediant of two rational numbers u and v

(numer(u) + numer(v))/(denom(u) + denom(v)) end:

'''convergents''' := proc(z)

# convergent list for z

local q;

convert(z,confrac,'q');

q end:

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

'''semiconvergents''' := proc(z)

# semiconvergent list for z

'''exlist'''('''convergents'''(z)) 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:

'''Ls''' := proc(q)

# large-small steps from mediant q

local u;

u := '''fareypair'''(q);

[denom(u[2]), denom(u[1])] end:

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

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

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

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

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