59 lines
1.8 KiB
Text
59 lines
1.8 KiB
Text
requires(64,true)
|
||
|
||
function congruence(atom a, c, m)
|
||
// Solves ax = c mod m
|
||
atom aa = mod(a,m),
|
||
cc = mod(c+a*m,m)
|
||
if aa == 1 then
|
||
return cc
|
||
end if
|
||
atom y = congruence(m, -cc, aa);
|
||
return (m*y+cc)/aa
|
||
end function
|
||
|
||
function frobenius(sequence unsorted_a)
|
||
sequence a = sort(deep_copy(unsorted_a))
|
||
assert(a[1]>=1);
|
||
if gcd(a)>1 then return "Undefined" end if
|
||
atom d12 = gcd(a[1],a[2]),
|
||
d13 = gcd(a[1]/d12,a[3]),
|
||
d23 = gcd(a[2]/d12,a[3]/d13),
|
||
{a1,a2,a3} = sort({a[1]/d12/d13,
|
||
a[2]/d12/d23,
|
||
a[3]/d13/d23})
|
||
atom rod = -1
|
||
if a1!=1 then
|
||
// Rødseth’s Algorithm
|
||
atom w1 = a1,
|
||
s0 = congruence(a2,a3,a1);
|
||
sequence s = {a1},
|
||
q = {}
|
||
while s0 != 0 do
|
||
s &= s0
|
||
atom s1 = iff(s0==1 ? 0 : s0-mod(w1,s0)),
|
||
q1 = (w1+s1)/s0;
|
||
q &= q1
|
||
w1 = s0
|
||
s0 = s1
|
||
end while
|
||
sequence p = {0, 1}
|
||
atom r = (s[2]*a2-p[2]*a3)/a1;
|
||
integer i = 1;
|
||
while r > 0 do
|
||
atom p_next = q[i]*p[i+1]-p[i];
|
||
p &= p_next
|
||
r = (s[i+2]*a2-p_next*a3)/a1;
|
||
i += 1;
|
||
end while
|
||
rod = -a1+a2*(s[i]-1)+a3*(p[i+1]-1)
|
||
- min(a2*s[i+1],a3*p[i])
|
||
end if
|
||
return sprintf("%d",rod*d12*d13*d23+a[1]*(d23-1)+a[2]*(d13-1)+a[3]*(d12-1))
|
||
end function
|
||
|
||
constant tests = {{6,9,20}, {6,7,20}, {1,3,20}, {10,5,18}, {5,17,44}, {2,4,6}, {3,6,15},
|
||
{12,14,17},{12,13,34},{5,9,21},{10,18,21},{71,98,99},{4,30,16},{12,12,13},
|
||
{7074047,8214596,9098139},{582795988,1753241221,6814151015},{6,15,1}}
|
||
for t in tests do
|
||
printf(1,"Maximum non-Mcnugget number using %v is: %s\n",{t,frobenius(t)})
|
||
end for
|