## 20114 Reputation

15 years, 277 days

## Serpentine paths in matrices and genera...

Maple 2018

This post is closely related to the previous one  https://www.mapleprimes.com/posts/210930-Numbrix-Puzzle-By-The-Branch-And-Bound-Method  which presents the procedure  NumbrixPuzzle   that allows you to effectively solve these puzzles (the text of this procedure is also available in the worksheet below).
This post is about generating these puzzles. To do this, we need the procedure  SerpentinePaths  (see below) , which allows us to generate a large number of serpentine paths in a matrix of a specified size, starting with a specified matrix element. Note that for a square matrix of the order  n , the number of such paths starting from [1,1] - position is the sequence  https://oeis.org/search?q=1%2C2%2C8%2C52%2C824&language=english&go=Search .

The required parameter of  SerpentinePaths procedure is the list  S , which defines the dimensions of the matrix. The optional parameter is the list  P  - this is the position of the number 1 (by default P=[1,1] ).
As an example below, we generate 20 puzzles of size 6 by 6. In exactly the same way, we can generate the desired number of puzzles for matrices of other sizes.

 > restart;
 > SerpentinePaths:=proc(S::list, P::list:=[1,1]) local OneStep, A, m, F, B, T, a; OneStep:=proc(A::listlist) local s, L, B, T, k, l; s:=max[index](A); L:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]]; T:=table(); k:=0; for l in L do if l[1]>=1 and l[1]<=S[1] and l[2]>=1 and l[2]<=S[2] and A[op(l)]=0 then k:=k+1; B:=subsop(l=a+1,A); T[k]:=B fi; od; convert(T, list); end proc; A:=convert(Matrix(S[], {(P[])=1}), listlist); m:=S[1]*S[2]-1; B:=[\$ 1..m]; F:=LM->ListTools:-FlattenOnce(map(OneStep, `if`(nops(LM)<=30000,LM,LM[-30000..-1]))); T:=[A]; for a in B do T:=F(T); od; map(convert, T, Matrix); end proc:
 > NumbrixPuzzle:=proc(A::Matrix) local A1, L, N, S, MS, OneStepLeft, OneStepRight, F1, F2, m, L1, p, q, a, b, T, k, s1, s, H, n, L2, i, j, i1, j1, R; uses ListTools; S:=upperbound(A); N:=nops(op(A)[3]); MS:=`*`(S); A1:=convert(A, listlist); for i from 1 to S[1] do for j from 1 to S[2] do for i1 from i to S[1] do for j1 from 1 to S[2] do if A1[i,j]<>0 and A1[i1,j1]<>0 and abs(A1[i,j]-A1[i1,j1])e<>0, Flatten(A1))); L1:=[`if`(L[1]>1,seq(L[1]-k, k=0..L[1]-2),NULL)]; L2:=[seq(seq(`if`(L[i+1]-L[i]>1,L[i]+k,NULL),k=0..L[i+1]-L[i]-2), i=1..nops(L)-1), `if`(L[-1]=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 then k:=k+1; T[k]:=subsop(m=a-1,A1); fi; od; convert(T, list); end proc; OneStepRight:=proc(A1::listlist) local s, M, m, k, T, s1; uses ListTools; s:=Search(a, Matrix(A1));  s1:=Search(a+2, Matrix(A1));   M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]]; T:=table(); k:=0; for m in M do if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 and `if`(a+2 in L, `if`(is(abs(s1[1]-m[1])+abs(s1[2]-m[2])>1),false,true),true) then k:=k+1; T[k]:=subsop(m=a+1,A1); fi; od; convert(T, list);    end proc; F1:=LM->ListTools:-FlattenOnce(map(OneStepLeft, LM)); F2:=LM->ListTools:-FlattenOnce(map(OneStepRight, LM)); T:=[A1]; for a in L1 do T:=F1(T); od; for a in L2 do T:=F2(T); od; R:=map(t->convert(t,Matrix), T); if nops(R)=0 then return `no solutions` else R fi; end proc:

Simple examples

 > SerpentinePaths([3,3]);  # All the serpentine paths for the matrix  3x3, starting with [1,1]-position SerpentinePaths([3,3],[1,2]);  # No solutions if the start with [1,2]-position SerpentinePaths([4,4]):  # All the serpentine paths for the matrix  4x4, starting with [1,1]-position nops(%); nops(SerpentinePaths([4,4],[1,2]));  # The number of all the serpentine paths for the matrix  4x4, starting with [1,2]-position nops(SerpentinePaths([4,4],[2,2]));  # The number of all the serpentine paths for the matrix  4x4, starting with [2,2]-position
 (1)

Below we find 12,440 serpentine paths in the matrix  6x6 starting from various positions (the set  L )

 > k:=0:  n:=6: for i from 1 to n do for j from i to n do k:=k+1; S[k]:=SerpentinePaths([n,n],[i,j])[]; od: od: L1:={seq(S[i][], i=1..k)}: L2:=map(A->A^%T, L1): L:=L1 union L2: nops(L);
 (2)

Further, using the list  L, we generate 20 examples of Numbrix puzzles with the unique solutions

 > T:='T': N:=20: M:=[seq(L[i], i=combinat:-randcomb(nops(L),N))]: for i from 1 to N do for k from floor(n^2/4) do T[i]:=Matrix(n,{seq(op(M[i])[3][j], j=combinat:-randcomb(n^2,k))}); if nops(NumbrixPuzzle(T[i]))=1 then break; fi; od:  od: T:=convert(T,list): T1:=[seq([seq(T[i+j],i=1..5)],j=[0,5,10,15])]: DocumentTools:-Tabulate(Matrix(4,5, (i,j)->T1[i,j]), fillcolor = "LightYellow", width=95):

The solutions of these puzzles

 > DocumentTools:-Tabulate(Matrix(4,5, (i,j)->NumbrixPuzzle(T1[i,j])[]), fillcolor = "LightYellow", width=95):
 >

For some reason, these 20 examples and their solutions did not load here.

Edit. I separately inserted these generated 20 puzzles as a picture:

## Numbrix Puzzle by the branch and bound m...

Maple 2018

In this post, the Numbrix Puzzle is solved by the branch and bound method (see the details of this puzzle in  https://www.mapleprimes.com/posts/210643-Solving-A-Numbrix-Puzzle-With-Logic). The main difference from the solution using the  Logic  package is that here we get not one but all possible solutions. In the case of a unique solution, the  NumbrixPuzzle procedure is faster than the  Numbrix  one (for convenience, I inserted the code for Numbrix procedure into the worksheet below). In the case of many solutions, the  Numbrix  procedure is usually faster (see all the examples below).

 > restart;
 > NumbrixPuzzle:=proc(A::Matrix) local A1, L, N, S, MS, OneStepLeft, OneStepRight, F1, F2, m, L1, p, q, a, b, T, k, s1, s, H, n, L2, i, j, i1, j1, R; uses ListTools; S:=upperbound(A); N:=nops(op(A)[3]); MS:=`*`(S); A1:=convert(A, listlist); for i from 1 to S[1] do for j from 1 to S[2] do for i1 from i to S[1] do for j1 from 1 to S[2] do if A1[i,j]<>0 and A1[i1,j1]<>0 and abs(A1[i,j]-A1[i1,j1])e<>0, Flatten(A1))); L1:=[`if`(L[1]>1,seq(L[1]-k, k=0..L[1]-2),NULL)]; L2:=[seq(seq(`if`(L[i+1]-L[i]>1,L[i]+k,NULL),k=0..L[i+1]-L[i]-2), i=1..nops(L)-1), `if`(L[-1]=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 then k:=k+1; T[k]:=subsop(m=a-1,A1); fi; od; convert(T, list); end proc;   OneStepRight:=proc(A1::listlist) local s, M, m, k, T, s1; uses ListTools; s:=Search(a, Matrix(A1));  s1:=Search(a+2, Matrix(A1));   M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]]; T:=table(); k:=0; for m in M do if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 and `if`(a+2 in L, `if`(is(abs(s1[1]-m[1])+abs(s1[2]-m[2])>1),false,true),true) then k:=k+1; T[k]:=subsop(m=a+1,A1); fi; od; convert(T, list);    end proc; F1:=LM->ListTools:-FlattenOnce(map(OneStepLeft, LM)); F2:=LM->ListTools:-FlattenOnce(map(OneStepRight, LM)); T:=[A1]; for a in L1 do T:=F1(T); od; for a in L2 do T:=F2(T); od; R:=map(t->convert(t,Matrix), T); if nops(R)=0 then return `no solutions` else R[] fi; end proc:
 > Numbrix := proc( M :: ~Matrix, { inline :: truefalse := false } ) local S, adjacent, eq, i, initial, j, k, kk, m, n, one, single, sol, unique, val, var, x;     (m,n) := upperbound(M);     initial := &and(seq(seq(ifelse(M[i,j] = 0                                    , NULL                                    , x[i,j,M[i,j]]                                   )                             , i = 1..m)                         , j = 1..n));     adjacent := &and(seq(seq(seq(x[i,j,k] &implies &or(NULL                                                        , ifelse(i>1, x[i-1, j, k+1], NULL)                                                        , ifelse(i1, x[i, j-1, k+1], NULL)                                                        , ifelse(j

Two simple examples

 > A:=<0,0,5; 0,0,0; 0,0,9>; # The unique solution NumbrixPuzzle(A); A:=<0,0,5; 0,0,0; 0,8,0>; # 4 solutions NumbrixPuzzle(A);
 (1)

Comparison with Numbrix procedure. The example is taken from

 > A:=<0, 0, 0, 0, 0, 0, 0, 0, 0;  0, 0, 46, 45, 0, 55, 74, 0, 0;  0, 38, 0, 0, 43, 0, 0, 78, 0;  0, 35, 0, 0, 0, 0, 0, 71, 0;  0, 0, 33, 0, 0, 0, 59, 0, 0;  0, 17, 0, 0, 0, 0, 0, 67, 0;  0, 18, 0, 0, 11, 0, 0, 64, 0;  0, 0, 24, 21, 0, 1, 2, 0, 0;  0, 0, 0, 0, 0, 0, 0, 0, 0>; CodeTools:-Usage(NumbrixPuzzle(A)); CodeTools:-Usage(Numbrix(A));
 memory used=7.85MiB, alloc change=-3.01MiB, cpu time=172.00ms, real time=212.00ms, gc time=93.75ms
 memory used=1.21GiB, alloc change=307.02MiB, cpu time=37.00s, real time=31.88s, gc time=9.30s
 (2)

In the example below, which has 104 solutions, the  Numbrix  procedure is faster.

 > C:=Matrix(5,{(1,1)=1,(5,5)=25}); CodeTools:-Usage(NumbrixPuzzle(C)): nops([%]); CodeTools:-Usage(Numbrix(C)):
 memory used=0.94GiB, alloc change=-22.96MiB, cpu time=12.72s, real time=11.42s, gc time=2.28s
 memory used=34.74MiB, alloc change=0 bytes, cpu time=781.00ms, real time=783.00ms, gc time=0ns
 >

## Three bucket problem and its generalizat...

Here is a classic puzzle:
You are camping, and have an 8-liter bucket which is full of fresh water. You need to share this water fairly into exactly two portions (4 + 4 liters). But you only have two empty buckets: a 5-liter and a 3-liter. Divide the 8 liters in half in as short a time as possible.

This is not an easy task and requires at least 7 transfusions to solve it.

The procedure  Pouring  solves a similar problem for the general case. Given n buckets of known volume and the amount of water in each bucket is known. Buckets can be partially filled, be full or be empty (of course the case when all is empty or all is full is excluded). With each individual transfusion, the bucket from which it is poured must be completely free, or the bucket into which it is poured must be completely filled. It is forbidden to pour water anywhere other than the indicated buckets.

Formal parameters of the procedure: BucketVolumes  is a list of bucket volumes,  WaterVolumes  is a list of water volumes in these buckets. The procedure returns all possible states that can occur during transfusions in the form of a tree (the initial state  WaterVolumes  is its root).

```restart;
Pouring:=proc(BucketVolumes::list(And(positive,{integer,float,fraction})),WaterVolumes::list(And(nonnegative,{integer,float,fraction})), Output:=graph)
local S, W, n, N, OneStep, j, v, H, G;
uses ListTools, GraphTheory;

n:=nops(BucketVolumes);
if nops(WaterVolumes)<>n then error "The lists should be the same length" fi;
if n<2 then error "Must have at least 2 buckets" fi;
if not `or`(op(WaterVolumes>~0)) then error "There must be at least one non-empty bucket" fi;
if BucketVolumes=WaterVolumes then error "At least one bucket should not be full" fi;
if not `and`(seq(WaterVolumes[i]<=BucketVolumes[i], i=1..n)) then error "The amount of water in each bucket cannot exceed its volume" fi;
W:=[[WaterVolumes]];

OneStep:=proc(W)
local w, k, i, v, V, k1, v0;
global L;
L:=convert(Flatten(W,1), set);
k1:=0;
for w in W do
k:=0; v:='v';
for i from 1 to n do
for j from 1 to n do
if i<>j and w[-1][i]<>0 and w[-1][j]<BucketVolumes[j] then k:=k+1; v[k]:=subsop(i=`if`(w[-1][i]<=BucketVolumes[j]-w[-1][j],0,w[-1][i]-(BucketVolumes[j]-w[-1][j])),j=`if`(w[-1][i]<=BucketVolumes[j]-w[-1][j],w[-1][j]+w[-1][i],BucketVolumes[j]),w[-1]); fi;
od;
od;
v:=convert(v,list);
if `and`(seq(v0 in L, v0=v)) then k1:=k1+1; V[k1]:=w else
for v0 in v do
if not (v0 in L) then k1:=k1+1; V[k1]:=[op(w),v0] fi;
od;
fi;
L:=L union convert(v,set);
od;
convert(V,list);
end proc:

S[0]:={};
for N from 1 do
H[N]:=(OneStep@@N)(W);
S[N]:=L;
if S[N-1]=S[N] then break fi;
od;
if Output=set then return L else
if Output=trails then interface(rtablesize=infinity);
return <H[N-1]> else
G:=Graph(seq(Trail(map(t->t[2..-2],convert~(h,string))),h=H[N-1]));
DrawGraph(G, style=tree, root=convert(WaterVolumes,string)[2..-2], stylesheet = [vertexcolor = "Yellow", vertexfont=[TIMES,20]], size=[800,500])  fi; fi;

end proc:
```

Examples of use:

Here is the solution to the original puzzle above. We see that at least 7 transfusions are
required to get equal volumes (4 + 4) in two buckets

Pouring([8,5,3], [8,0,0]);

With an increase in the number of buckets, the number of solutions is extremely
increased. Here is the solution to the problem: is it possible to equalize the amount of water (7+7+7+7) in the following example?

Pouring([14,10,9,9],[14,10,4,0]);
S:=Pouring([14,10,9,9],[14,10,4,0], set);
is([7,7,7,7] in S);
nops(S);

## Bugs in maximize and minimize commands...

Maple 2018

Yesterday, I accidentally discovered a nasty bug in a fairly simple example:

restart;
Expr:=a*sin(x)+b*cos(x);
maximize(Expr, x=0..2*Pi);
minimize(Expr, x=0..2*Pi);

I am sure the correct answers are  sqrt(a^2+b^2)  and  -sqrt(a^2+b^2)  for any real values  a  and  b .  It is easy to prove in many ways. The simplest method does not require any calculations and can be done in the mind. We will consider  Expr  as the scalar product (or the dot product) of two vectors  <a, b>  and  <sin(x), cos(x)>, one of which is a unit vector. Then it is obvious that the maximum of this scalar product is reached if the vectors are codirectional and equals to the length of the first vector, that is, sqrt(a^2+b^2).

Bugs in these commands were noted by users and earlier (see search by keywords bug, maximize, minimize) but unfortunately are still not fixed.

## Both area and perimeter in half...

In this post an interesting geometric problem is solved: for an arbitrary convex polygon, find a straight line that divides both the area and the perimeter in half. The following results on this problem are well known:
1. For any convex polygon there is such a straight line.
2. For any convex polygon into which a circle can be inscribed, in particular for any triangle, the desired line must pass through the center of the inscribed circle.
3. For a triangle, the number of solutions can be 1, 2, or 3.
4. If the polygon has symmetry with respect to a point, then any straight line passing through this point is a solution.

The procedure called  InHalf  (the code below) symbolically solves this problem. The formal parameter of the procedure is the list of coordinates of the vertices of a convex polygon (vertices must be passed opposite or clockwise). The procedure returns all solutions in the form of a list of pairs of points, lying on the perimeter of the polygon, that are the ends of segments that implement the desired dividing.

 > restart;
 > InHalf:=proc(V::listlist) local L, n, a, b, M, N, i, j, P, Q, L1, L2, Area, Area1, Area2, Perimeter, Perimeter1, Perimeter2, sol, m, k, Sol; uses LinearAlgebra, ListTools; L:=map(convert,[V[],V[1]],rational); n:=nops(L)-1; a:=<(V[2]-V[1])[1],(V[2]-V[1])[2],0>; b:=<(V[n]-V[1])[1],(V[n]-V[1])[2],0>; if is(CrossProduct(a,b)[3]<0) then L:=Reverse(L) fi; M:=[seq([L[i],L[i+1]], i=1..n)]: N:=0; for i from 1 to n-1 do for j from i+1 to n do P:=map(t->t*(1-s),M[i,1])+map(t->t*s,M[i,2]); Q:=map(s->s*(1-t),M[j,1])+map(s->s*t,M[j,2]); L1:=[P,L[i+1..j][],Q,P]; L2:=[Q,L[j+1..-1][],L[1..i][],P,Q]; Area:=L->(1/2)*add(L[k, 1]*L[k+1, 2]-L[k, 2]*L[k+1, 1], k = 1 .. nops(L)-1); Area1:=Area(L1); Area2:=Area(L2); Perimeter:=L->add(sqrt((L[k,1]-L[k+1,1])^2+(L[k,2]-L[k+1,2])^2), k=1..nops(L)-2); Perimeter1:=Perimeter(L1); Perimeter2:=Perimeter(L2); sol:=[solve({Area1=Area2,Perimeter1=Perimeter2,s>=0,s<1,t>=0,t<1}, {s,t}, explicit)] assuming real; if sol<>[] then m:=nops(sol); for k from 1 to m do N:=N+1; if nops(sol[k])=2 then Sol[N]:=simplify(eval([P,Q],sol[k])) else Sol[N]:=simplify(eval([P,Q],s=t)) fi; od; fi; od; od; Sol:=convert(Sol, list); `if`(indets(Sol)={},Sol,op([Sol,t>=0 and t<1])); end proc:

Examples of use

 > # For the Pythagorean triangle with sides 3, 4, 5, we have a unique solution L:=[[4,3],[4,0],[0,0]]: P:=InHalf(L); plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), scaling=constrained);
 > # For an isosceles right triangle, there are 3 solutions. We see that all the cuts pass through the center of the inscribed circle L:=[[0,0],[4,0],[4,4]]: InHalf(L); P:=InHalf(L); r:=(4+4-4*sqrt(2))/2: a:=4-r: b:=r: plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), plot([r*cos(t)+a,r*sin(t)+b, t=0..2*Pi], color=blue), scaling=constrained);
 > # There are 3 solutions for the quadrilateral below L:=[[0,0],[4.5,0],[4,3],[0,2]]: P:=InHalf(L); plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), scaling=constrained);
 > # There are infinitely many solutions for a polygon with a center of symmetry. Any cut through the center solves the problem. The picture shows 2 solutions. L:=[[1,0],[1+2*sqrt(3),2],[2*sqrt(3),sqrt(3)+2],[0,sqrt(3)]]: P:=InHalf(L); plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(eval(P[1],t=1/3),  color=red), scaling=constrained);
 >