gap> Multiplicity(S3);
13
```
Update:
Looks like this is the code: [numericalsgps/gap/basics2.gi at master · gap-packages/numericalsgps · GitHub](https://github.com/gap-packages/numericalsgps/blob/master/gap/basics2.gi#L97-L120)
As I understand it they compute all elements below the Frobenius number for both semigroups and intersect them. I.e. basically brute force, so I guess not so much to find there after all.
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [October 24, 2024, 8:19am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/9 "2024-10-24T08:19:54Z")
For what it’s worth, I compared solutions of my brute force approach and the ILP provided by Eric. I randomly generated a, b, c, d \in \{ 1, \ldots, 10\,000 \}, obtained a solution from the ILP and my program and measured how long it took, compared the solutions and recorded if they differed.
After 500 experiments, I got the following results. ILP took 0.475s on average with standard deviation of 1.101s. My approach took 0.0004s on average with standard deviation of 0.0006s. I found one problem which resulted in different n: a=1985, b=1212, c=9941, d=7008. However, ILP gave higher (worse) n: 341420 = 172 \* 1985 + 0 \* 1212 = 28 \* 9941 + 9 \* 7008; my program gave 337994 = 166 \* 1985 + 7 \* 1212 = 34 \* 9941 + 0 \* 7008. I believe this could be due to the time limit I imposed on the ILP solver (20s).
Also, there can be multiple optimizers. For example:
```julia
# Format is:
# Problem
# ILP solution of the form n = i * a + j * b = k * c + l * d
# My solution of the same form
a=5 b=11 c=2381 d=982
982 = 192 * 5 + 2 * 11 = 0 * 2381 + 1 * 982
982 = 5 * 5 + 87 * 11 = 0 * 2381 + 1 * 982
a=1785 b=4418 c=17 d=20
1785 = 1 * 1785 + 0 * 4418 = 105 * 17 + 0 * 20
1785 = 1 * 1785 + 0 * 4418 = 5 * 17 + 85 * 20
a=8 b=35 c=401 d=3682
401 = 37 * 8 + 3 * 35 = 1 * 401 + 0 * 3682
401 = 2 * 8 + 11 * 35 = 1 * 401 + 0 * 3682
a=8945 b=6362 c=62 d=15
6362 = 0 * 8945 + 1 * 6362 = 91 * 62 + 48 * 15
6362 = 0 * 8945 + 1 * 6362 = 1 * 62 + 420 * 15
a=7122 b=3427 c=45 d=58
3427 = 0 * 7122 + 1 * 3427 = 71 * 45 + 4 * 58
3427 = 0 * 7122 + 1 * 3427 = 13 * 45 + 49 * 58
```
> **Code**
>
> ```julia
> using JuMP, HiGHS
> using Statistics
>
> function solve(x0, y0, a, b, c)
> @assert x0 > 0 || y0 > 0
>
> x0 = x0 * c
> y0 = y0 * c
>
> if x0 < 0
> t = ceil(Int, -x0 / b)
> x = x0 + b * t
> y = y0 - a * t
> elseif y0 < 0
> t = ceil(Int, -y0 / a)
> x = x0 - b * t
> y = y0 + a * t
> else
> x = x0
> y = y0
> end
>
> if x ≥ 0 && y ≥ 0
> return x, y
> else
> return nothing
> end
> end
>
> function solve_brute(a, b, c, d)
> _, i0, j0 = gcdx(a, b)
> _, k0, l0 = gcdx(c, d)
>
> n_min = max(min(a, b), min(c, d))
> n_max = a * c + b * d
>
> for n in n_min:n_max
> sol_ab = solve(i0, j0, a, b, n)
> isnothing(sol_ab) && continue
>
> sol_cd = solve(k0, l0, c, d, n)
> isnothing(sol_cd) && continue
>
> i, j = sol_ab
> k, l = sol_cd
>
> @assert a * i + b * j == n
> @assert c * k + d * l == n
>
> return (; n, i, j, k, l)
> end
>
> return nothing
> end
> function solve_ilp(a, b, c, d)
> model = Model(HiGHS.Optimizer)
> set_silent(model)
> set_time_limit_sec(model, 20)
> @variables(model, begin
> i >= 0, Int
> j >= 0, Int
> k >= 0, Int
> l >= 0, Int
> n >= 1, Int
> end)
> @constraint(model, n == i * a + j * b)
> @constraint(model, n == k * c + l * d)
> @objective(model, Min, n)
> optimize!(model)
> is_solved_and_feasible(model) || return nothing
>
> n = round(Int, value(n))
> i = round(Int, value(i))
> j = round(Int, value(j))
> k = round(Int, value(k))
> l = round(Int, value(l))
>
> @assert a * i + b * j == n
> @assert c * k + d * l == n
>
> return (; n, i, j, k ,l)
> end
>
> function generate_abcd(range)
> a, b, c, d = rand(range, 4)
> b += (a == b)
> d += (c == d)
> p = gcd(a, b)
> q = gcd(c, d)
> a, b = a ÷ p, b ÷ p
> c, d = c ÷ q, d ÷ q
> return a, b, c, d
> end
>
> function test(;m=10, high=10)
> times_ilp = Float64[]
> times_brute = Float64[]
>
> tricky_tasks = NTuple{12, Int64}[]
> wrong_tasks = NTuple{12, Int64}[]
>
> for i in 1:m
> a, b, c, d = generate_abcd(1:high)
>
> result_ilp = @timed solve_ilp(a, b, c, d)
> result_brute = @timed solve_brute(a, b, c, d)
>
> push!(times_ilp, result_ilp.time)
> push!(times_brute, result_brute.time)
>
> sol_ilp = result_ilp.value
> sol_brute = result_brute.value
>
> if sol_ilp != sol_brute
> if sol_ilp.n != sol_brute.n
> push!(wrong_tasks, (a, b, c, d,
> sol_ilp.i, sol_ilp.j, sol_ilp.k, sol_ilp.l,
> sol_brute.i, sol_brute.j, sol_brute.k, sol_brute.l))
> else
> push!(tricky_tasks, (a, b, c, d,
> sol_ilp.i, sol_ilp.j, sol_ilp.k, sol_ilp.l,
> sol_brute.i, sol_brute.j, sol_brute.k, sol_brute.l))
> end
> end
>
> if (i % 10) == 0
> println("$i tasks done")
> end
> end
>
> println("ILP times: ", round(mean(times_ilp), digits=4), " ± ", round(std(times_ilp), digits=4))
> println("BRT times: ", round(mean(times_brute), digits=4), " ± ", round(std(times_brute), digits=4))
> println()
> println("Wrong tasks:")
> for (a, b, c, d, ilp_i, ilp_j, ilp_k, ilp_l, brute_i, brute_j, brute_k, brute_l) in wrong_tasks
> ilp_n = ilp_i * a + ilp_j * b
> brute_n = brute_i * a + brute_j * b
> println("a=$a b=$b c=$c d=$d")
> println(" $ilp_n = $ilp_i * $a + $ilp_j * $b = $ilp_k * $c + $ilp_l * $d")
> println(" $brute_n = $brute_i * $a + $brute_j * $b = $brute_k * $c + $brute_l * $d")
> end
> println()
> println("Tricky tasks:")
> for (a, b, c, d, ilp_i, ilp_j, ilp_k, ilp_l, brute_i, brute_j, brute_k, brute_l) in tricky_tasks
> ilp_n = ilp_i * a + ilp_j * b
> brute_n = brute_i * a + brute_j * b
> println("a=$a b=$b c=$c d=$d")
> println(" $ilp_n = $ilp_i * $a + $ilp_j * $b = $ilp_k * $c + $ilp_l * $d")
> println(" $brute_n = $brute_i * $a + $brute_j * $b = $brute_k * $c + $brute_l * $d")
> end
> end
>
> ```
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [October 24, 2024, 10:06am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/10 "2024-10-24T10:06:26Z")
Here’s another brute force solution by iteratively generating the semigroups until the smallest common element is found.
```julia
using DataStructures
function solve_heap(a, b, c, d)
S1 = BinaryMinHeap([a, b])
S2 = BinaryMinHeap([c, d])
m1 = min(a, b)
m2 = min(c, d)
while m1 != m2
if m1 < m2
while first(S1) == m1
pop!(S1)
end
push!(S1, m1 + a)
push!(S1, m1 + b)
m1 = first(S1)
end
if m2 < m1
while first(S2) == m2
pop!(S2)
end
push!(S2, m2 + c)
push!(S2, m2 + d)
m2 = first(S2)
end
end
return m1
end
```
(`Set` can of course be used instead of heap, but is considerably slower.)
Edit: Fixed corner case.
---
### Author:  [@ambiso](https://discourse.julialang.org/u/ambiso)
#### Post date: [October 24, 2024, 10:28am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/11 "2024-10-24T10:28:47Z")
There’s probably some connection to lattice reduction.
Here’s a first attempt using sagemath:
```python
from sage.modules.free_module_integer import IntegerLattice
a = random_prime(2**14)
b = random_prime(2**14)
c = random_prime(2**14)
d = random_prime(2**14)
print(f"{a=} {b=} {c=} {d=}")
N = 1000 # some large enough multiplier
L = identity_matrix(4).augment(N * vector([a, b, -c, -d]))
print(L)
L = IntegerLattice(L)
v = L.shortest_vector() # may have negative coefficients
print(f"{v=}")
print(L)
B = L.basis()
er = 10
smallest_n = Infinity
# some terrible enumeration to (try to) eliminate negative coefficients
for x1 in range(-er, er+1):
for x2 in range(-er, er+1):
for x3 in range(-er, er+1):
for x4 in range(-er, er+1):
v = x1*B[0] + x2*B[1] + x3*B[2] + x4*B[3]
l = v[0]*a + v[1] * b
r = v[2]*c + v[3] * d
if l != r:
continue
n = l
if n > 0 and all(x >= 0 for x in v[:4]) and n < smallest_n:
smallest_n = n
print(n, v)
```
which prints:
```julia
a=1231 b=10357 c=5051 d=2087
[1 0 0 0 1231000]
[0 1 0 0 10357000]
[0 0 1 0 -5051000]
[0 0 0 1 -2087000]
v=(-1, 4, 3, 12, 0)
Free module of degree 5 and rank 4 over Integer Ring
User basis matrix:
[-1 4 3 12 0]
[-19 3 4 -6 0]
[5 15 39 -17 0]
[-8 -1 -4 0 -1000]
221367 (163, 2, 43, 2, 0)
123100 (100, 0, 19, 13, 0)
72712 (17, 5, 2, 30, 0)
```
the last line shows the best solution that was found.
edit: in the reduced basis only the last row vector has an entry in the last column, it should be possible to skip it during enumeration (as an optimization). A sufficiently large multiplier would maybe ensure this? I’m not sure.
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [October 24, 2024, 11:11am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/12 "2024-10-24T11:11:44Z")
> [@barucden](#):
>
> ```julia
> function optimize(a, b, c, d)
> n_min = max(min(a, b), min(c, d))
> n_max = a * c + b * d
>
> ```
ia+jb=lc+kd has 2 symmetric solutions i=c,j=d,l=a,k=b or i=d,j=c,l=b,k=a, so maybe `n_max = min(a*c+b*d, b*c+a*d)` is more appropriate?
But since i,j,k,l can be 0, there are more easy symmetric solutions:
i=0,l=0, n=bd \\ i=0,k=0, n=bc\\ j=0,l=0, n=ad\\ j=0,k=0, n=ac
among which is a value lower than `n_max` above.
You start the search at `n_min` and would hit one of these easy symmetric solutions no matter how much larger `n_max` is, so I don’t expect this to help runtime at all.
Tried a toy example that would, but `optimize(2, 3^23, 3, 2^45) == nothing`. The solution should be n=6=3\*2+0\*3^{23}=2\*3+0\*2^{45} if I’m reading this right
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [October 24, 2024, 12:12pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/13 "2024-10-24T12:12:29Z")
Try the version from the benchmark (`solve_brute`). I guess I made some changes, which I have forgotten from yesterday 😄
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [October 24, 2024, 12:30pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/14 "2024-10-24T12:30:51Z")
Ah I missed that details block, otherwise I would’ve just commented on the same lines from `solve_brute`. `solve_brute(2, 3^23, 3, 2^45) == nothing` as well.
I neglected to mention that the older `solve(2, 3^23, 6) == (3, 0)` and `solve(3, 2^45, 6) == (2, 0)`, which is why I was confused that `optimize` resulted in `nothing`. It seems like it should’ve reached the solution n=6, but I haven’t tried debugging it.
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [October 24, 2024, 12:33pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/15 "2024-10-24T12:33:42Z")
It fails due to integer overflow from `3^23 * 2^45`. But there’s no need for the upper bound, easier to just set it to `typemax(Int)`.
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [October 24, 2024, 12:40pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/16 "2024-10-24T12:40:33Z")
Switching to the symmetric `n_max = b*c+a*d`, now `optimize(2, 3^23, 3, 2^45) == (3, 0, 2, 0)`. That’s not going to hold up in unit testing, though. The least known upper bound could definitely be left in a comment.
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [October 24, 2024, 12:56pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/17 "2024-10-24T12:56:56Z")
If you want a trivial upper bound from cross products, you can go with `min(a, b) * min(c, d)`.
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [October 24, 2024, 1:43pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/18 "2024-10-24T13:43:23Z")
I’ll summarize my approach so that it is more convenient for people to improve it, if someone is interested.
We are given positive integers a, b, c, d such that \gcd(a, b) = 1 and \gcd(c, d) = 1. We are trying to find smallest positive n\>0 such that
n = i a + j b = k c + ld,
where i, j, k, l \in \mathbb{N} are non-negative variables.
My approach is based on finding natural solutions of the linear equations. Let’s first focus on the first one, n = ia + jb.
From the Bézout’s identity, we can find i\_0, j\_0 \in \mathbb{Z} such that
ai\_0 + b j\_0 = 1.
Multiplying by some n \in \mathbb{N}, we get
a (ni\_0) + b(n j\_0) = n.
Since either i\_0 or j\_0 is negative (unless a = b = 1), we don’t have a natural solution. Suppose i\_0 is negative. We can write
a(n i\_0 + b m) + b(n j\_0 - am) = n,
where m \in \mathbb{N} is some number. The goal now is to find the smallest n \in \mathbb{N} such that there is some m \in \mathbb{N} that makes both ni\_0 + bm and n j\_0 - am non-negative. Formally:
n i\_0 + bm \geq 0 \rightarrow m \geq -\frac{ni\_0}{b}, \\ n j\_0 - am \geq 0 \rightarrow m \leq \frac{nj\_0}{a},
which gives together
n s \leq m \leq n t,
where s = -\frac{i\_0}{b} and t = \frac{j\_0}{a}. Now, if j\_0 was negative, we would get s = -\frac{j\_0}{a} and t = \frac{i\_0}{b}.
So, our task transforms into finding the smallest n \in \mathbb{N} for which there is an m \in \mathbb{N} such that n s \leq m \leq n t. Subtract ns to get
0 \leq m - ns \leq n ( t - s ),
in which we set m = \lceil ns \rceil to satisfy the first inequality. It leads to finding the smallest n \in \mathbb{N} such that
\begin{aligned} \lceil ns \rceil - ns &\leq n ( t - s) \\ \lceil ns \rceil &\leq n t \end{aligned}
If we make the very same argument but for the equation n = kc + ld, we obtain
\lceil nu \rceil \leq nv,
where u = \begin{cases} -k\_0 / d & \text{if } k\_0 \< 0 \\ -l\_0 / c & \text{if } l\_0 \< 0 \end{cases} and v = \begin{cases} l\_0 / c & \text{if } k\_0 \< 0 \\ k\_0 / d & \text{if } l\_0 \< 0 \end{cases} with the Bézout’s coefficients k\_0, l\_0 such that c k\_0 + d l\_0 = 1.
So the algorithm is:
1. Initialize n := 1
2. Find smallest n\_1 \geq n for which \lceil n\_1 s \rceil \leq n\_1 t holds
3. Find smallest n\_2 \geq n for which \lceil n\_2 u \rceil \leq n\_2 v holds
4. Set n := \max\{ n\_1, n\_2 \}
5. If n\_1 \neq n\_2, go to 2.
6. Compute solution.
The solution can be computed as
\begin{aligned} i &= i\_0 n - \text{sign}(i\_0) b m\_1, \\ j &= j\_0 n + \text{sign}(i\_0) a m\_1, \\ k &= k\_0 n - \text{sign}(k\_0) d m\_2, \\ l &= l\_0 n + \text{sign}(k\_0) c m\_2, \end{aligned}
where m\_1 = \lceil ns \rceil and m\_2 = \lceil nu \rceil.
The problem with this approach is that for some inputs, we have to iterate over too many n before finding a match n\_1 = n\_2. We could improve that by initializing n with a value that is closer to the optimum. Another improvement could be in steps 2 and 3, where now I sequentially check if the condition holds (maybe we could do skip some values).
Code here:
> **Code**
>
> ```julia
> function next_lowest(u, v, n)
> while n <= typemax(Int)
> if n * v - ceil(n * u) ≥ 0
> return n
> end
> # Improvement: do larger steps
> n += 1
> end
> throw(ErrorException("not found"))
> end
>
> function solve(a, b, c, d)
> _, i0, j0 = gcdx(a, b)
> _, k0, l0 = gcdx(c, d)
>
> s = ifelse(i0 < 0, -i0 // b, -j0 // a)
> t = ifelse(i0 < 0, j0 // a, i0 // b)
>
> u = ifelse(k0 < 0, -k0 // d, -l0 // c)
> v = ifelse(k0 < 0, l0 // c, k0 // d)
>
> # Improvement: Initialize better
> n = 1
>
> while true
> n1 = next_lowest(s, t, n)
> n2 = next_lowest(u, v, n)
>
> n = max(n1, n2)
>
> if n1 == n2
> break
> end
> end
>
> println("Solved for n=$n")
>
> m1 = ceil(n * s)
> m2 = ceil(n * u)
>
> i = Int(i0 * n - sign(i0) * b * m1)
> j = Int(j0 * n + sign(i0) * a * m1)
> k = Int(k0 * n - sign(k0) * d * m2)
> l = Int(l0 * n + sign(k0) * c * m2)
>
> println("i=", i)
> println("j=", j)
> println("k=", k)
> println("l=", l)
>
> end
>
> ```
---
### Author:  [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)
#### Post date: [October 24, 2024, 2:32pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/19 "2024-10-24T14:32:52Z")
> [@StefanKarpinski](#):
>
> It’s easy to put an upper bound on this
is an “easy” upper bound `lcm(a+b,c+d)*(a+b)?`
or
`n<= min(a,b)*min(c,d)`
```julia
a
### Author:  [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)
#### Post date: [October 24, 2024, 2:37pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/20 "2024-10-24T14:37:13Z")
I think moving this to general usage would be a good idea. ILP is likely neither an exclusive nor optimal way to solve the problem.
[Next page](https://discourse.julialang.org/t/small-integer-programming-problem/121647.md?page=2)