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