with(LinearAlgebra); with(simplex); with(combinat); Elig := proc (l) local i; for i to nops(l) do if evalf(l(i)) < 0 then return false end if end do; true end proc; Ver := proc (l) local i, s; s := {}; for i to nops(l) do if 0 < evalf(l[i]) then s := `union`(s, {i}) end if end do; s end proc; mop := proc (x, n) modp(x-1, n)+1 end proc; pm := proc (l, S) local d, i; d := add(l[i]^2, i in S); sqrt(d) end proc; IsCell := proc (s, L) local i; for i to nops(L) do if `subset`(s, L[i]) then return true end if end do; false end proc; Par := proc (l1, l2, L) Partition(l1, l2, [Ver(l1)], [Ver(l2)], L) end proc; Partition := proc (l1, l2, la, lb, L) local A, a, B, b, c, d, F, i, j, i1, j1, k, m, P, Q, s, t; k := nops(la); if k <> nops(lb) or k = 0 then error "invalid input" end if; F := `union`(seq(L[i], i = 1 .. nops(L))); a := `union`(seq(la[i], i = 1 .. k)); b := `union`(seq(lb[i], i = 1 .. k)); if `intersect`(a, b) <> {} then return Partition(l1, l2, map(proc (x) options operator, arrow; `minus`(x, b) end proc, la), map(proc (x) options operator, arrow; `minus`(x, a) end proc, lb), L) end if; for m to k do A := la[m]; B := lb[m]; if `not`(`subset`(A, F)) or `not`(`subset`(B, F)) then error "invalid input" end if; if pm(l1, A) = 0 or pm(l2, B) = 0 then next end if; P := powerset(A); Q := powerset(B); d := 1; c := 0; for i in P do i1 := `minus`(A, i); for j in Q do j1 := `minus`(B, j); if `not`(IsCell(`union`(seq(la[s], s = m+1 .. k), seq(lb[s], s = 1 .. m-1), i1, j), L)) then next end if; if evalf(pm(l1, i)^2/pm(l1, A)^2+pm(l2, j1)^2/pm(l2, B)^2-d) < 0 then d := pm(l1, i)^2/pm(l1, A)^2+pm(l2, j1)^2/pm(l2, B)^2; c := 1; t := [i, i1, j, j1] end if end do end do; if c = 1 then return Partition(l1, l2, [op(la[1 .. m-1]), t[1], t[2], op(la[m+1 .. k])], [op(lb[1 .. m-1]), t[3], t[4], op(lb[m+1 .. k])], L) end if end do; [la, lb] end proc; empdeg := proc (l) local d, i; d := 0; for i to nops(l) do if l[i] = {} then d := d+1 end if end do; d end proc; Geo := proc (l1, l2, L) local a1, A, b1, B, d1, d2, i, j, k, l, m, n, P, t, v, V; A := Ver(l1); B := Ver(l2); n := nops(l1); if IsCell(`union`(A, B), L) then return [l1, l2] end if; P := Partition(l1, l2, [A], [B], L); k := nops(P[1]); V := Array([seq(0, i = 1 .. k+2)]); V[1] := l1; V[k+2] := l2; for t from 2 to k+1 do m := V[t-1]; v := Array([seq(0, i = 1 .. n)]); a1 := P[1][t-1]; b1 := P[2][t-1]; d1 := pm(m, a1); d2 := pm(l2, b1); for i to n do if i in `intersect`(A, B) or i in `union`(seq(P[2][j], j = 1 .. t-2)) then v[i] := simplify((d2*m[i]+d1*l2[i])/(d1+d2)) end if; for j from t to k do if i in P[1][j] then v[i] := simplify(m[i]*(d2*pm(m, P[1][j])-d1*pm(l2, P[2][j]))/(pm(m, P[1][j])*(d1+d2))) end if end do end do; l := convert(v, list); V[t] := l end do; l := convert(V, list); l := l[empdeg(P[1])+1 .. -empdeg(P[2])-1]; l end proc;