Carl Love

## 26892 Reputation

11 years, 308 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Simple efficient recursive procedure...

A fundamental motivation of the discipline of Combinatorics is that counting combinatorial structures can often be done using far less computational effort than generating the structures. Counting partitions is like that.

It can be done quite efficiently with a very short recursive procedure:

```restart:
CountSums:= (k::posint, N::set(posint))->
((k,n)-> (thisproc(args):=
`if`(k*n > 0, thisproc(k, n-1) + thisproc(k-N[n], n), `if`(k=0, 1, 0))
))(k, nops(N))
:
CodeTools:-Usage(CountSums(1000, {2,3,6,9}));
memory used=421.16KiB, alloc change=0 bytes,
cpu time=15.00ms, real time=7.00ms, gc time=0ns

526848
CountSums(9, {2,3,6,9});
4
CountSums(8, {2,3,6,9});
3
#Verify correctness:
CodeTools:-Usage(combinat:-numbpart(100));
memory used=2.14MiB, alloc change=0 bytes,
cpu time=31.00ms, real time=33.00ms, gc time=0ns

190569292

CodeTools:-Usage(CountSums(100, {\$100}));
memory used=1.18MiB, alloc change=0 bytes,
cpu time=16.00ms, real time=13.00ms, gc time=0ns

190569292

```

So, for this example, my procedure's time beats the library procedure combinat:-numbpart.

@Pepini You wrote:

• I thought about use parametrization of torus:
plot3d([(R + r*cos(t)*cos(s), (R + r*cos(t)*sin(s), r*sin(t)], t = 0 .. 2*Pi, s = 0 .. 2*Pi)

You are missing some parentheses in that parametrization: The R also needs to be multiplied by the cos(s) or sin(s). That omission may be why your plot wasn't satisfactory as a surface of revolution. Here's a correction:

r:=  t-> 1 + sin(8 + sin(16 + sin(32*t)))/4:  #This works for any polar curve in form r = f(t)
R:= 2:
plot3d(
[(R + (r*cos)(t))*~(cos,sin)(s), (r*sin)(t)], t= -Pi..Pi, s= -Pi..Pi,
scaling= constrained,
#usually a better option than using 'view'
grid= [101,49]   #Increase points used for t to show wrinkling
);

This is a figure homeomorphic to an ordinary torus, not the "flat torus" that @sand15 gave links to a PDF about and you showed a picture of. And it has only one level and direction of what those authors call "corrugation" (you might call it "wrinkling"). I think that 2 directions would be the minimum for a reasonable plot, but 1 direction is the most that can be accommodated by a surface of revolution.

## Matrix, erf...

Here's a Matrix with all the information you need. You could display it more elegantly with DocumentTools, I suppose.

 > restart:
 > interface(rtablesize= [37,11]):
 > <     ^%T,     <          |         Matrix((36,10), (i,j)-> evalf[4]((1+erf((10*i+j-11)/100/sqrt(2)))/2))     > >;

 >

## From 1st principles, no packages...

Perhaps you will understand better if things are made more explicit rather than using the geometry package. MaplePrimes won't let me display this worksheet inline at the moment, but you can download it...

...or just copy and paste this code into a Maple session:

```Dist:= (P,Q)-> sqrt(add((P-Q)^~2)):  #Euclidean distance between points P and Q
Mid:= (P,Q)->  (P+Q)/2:              #Midpoint between P and Q

P:= <x,y>:                       #An arbitrary symbolic point on the ellipse
F:= [<0,-5>, <5, 0>]:  MA:= 12:  #Foci and major-axis length

(* An ellipse's defining property: The sum of the distances from any of its
points to its foci is constant. I'll call that constant C for the moment.
We'll use that other piece of given information, the length of the major axis,
to figure out C in a moment. *)

Ell_eqn:= C = add(map(Dist, F, P));

# Perhaps you'll accept without further explanation that the line through
# the foci (which is collinear with the major axis) is
MA_line:= y = x - 5;

(* The next step to getting C is finding numerically any point on the ellipse.
Two such points are the points on that line that are at distance MA/2 from the
midpoint of the foci. We solve for them: *)

%solve({MA_line, Dist(P, Mid(F[])) = MA/2}, explicit);
sols:= value(%);

# As suspected, there are 2 solutions. We only need 1, and I arbitrarily choose
# the first. To find C, we just use those x and y in the ellipse's equation:

eval(Ell_eqn, sols[1]);
eval(Ell_eqn, %):  #Put that in for C

# Subtract one side of the equation from the other to make an expression
# implicitly equated to 0:

(lhs-rhs)(%);
evala(Norm(%));      #Convert that to a polynomial
(primpart/sign)(%);  #Cancel common integer factors of the coefficients
#Final:
119*x^2 - 50*x*y + 119*y^2 - 720*x + 720*y - 1584

```

## Simplify, or UseHardwareFloats:= false...

Your expression involves extreme underflows and overflows that can't be handled by the default hardware floats. So, the easy correction is

UseHardwareFloats:= false:

An alternative is to simplify the expression to

f2:= exp(5002.646494 + 2499*ln(t) - 5000.*sqrt(t));

Then hardware floats can be used by plot.

## `if` as index...

My preferred syntax for situations like this is decisions[`if`(1 > T, 1, 2)].

## GAMMA...

I told you before that you need to change Gamma to GAMMA. You said that you did, but I still see Gamma.

## Try A[b,c]...

Try entering it as A[b,c]. This makes it an indexed variable. The output will be subscripted like you want. An indexed variable is not quite the same thing as the atomic variable that you were trying to create, but it'll work if you do not assign values to Ab, or c.

The characters that you see are not apostrophes. Note that they lean to the left. They're called "back quotes" or "accents graves" or, in Maple, "name quotes". They let you make any string of characters, no matter how weird, into a variable.

However, you shouldn't be seeing those back quotes in your output. If I enter

`A__b,c`;

then my output appears without quotes as Ab,c

## plots:-setoptions3d...

You can set defaut options for all 3-D plots with the command

plots:-setoptions3d('axesfont'= ["TimesNewRoman", 24], 'labelfont'= ["TimesNewRoman", 24]);

These options can be overridden, if desired, on any individual plot command simply by including a different version of the option in the individual command.

You can use plots:-setoptions to do the same thing for 2-D plots.

These can be put into an initialization file which is automatically run whenever a new worksheet is loaded or a restart is executed.

## rcurry...

One way:

(Int = rcurry(int, 'AllSolutions'))(1/x, x= a..2);

## Done, by multithreaded use of the Iterat...

I have verified conclusively that the bondage number of your example graph is 7 by showing that a certain 7-subset of edges is a bondage set (that part is extremely easy), and using multithreaded processing with the Iterator package to verify exhaustively that none of the 439,381,916 edge subsets with less than 7 edges is a bondage set. With careful multithreading (using 8 threads on my computer), I reduced the average real time to verify whether an edge subset is a bondage set to about 17.3 microseconds (equivalent to about 3.4 million subsets checked per minute).

I'll put the main code and the worksheet separately because the main code is in a Code Edit Region within the worksheet.

```MinBondageSet:= module()
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Jun-22`;
local
#causes the creation of a separate copy for each thread.
;
export
#All of this module's exports are essentially "static", but they're not declared
#"static" because this module is not an "object"
VL, n, V, A0, #VL = vertex labels, n= # of vertices, V={\$n}, A0 = list of vertex neighborhoods
edges, NE,
DS, #minimum dominating sets of input graph, as a listlist

#returns the set of minimum dominating sets of G as a listlist:
MinDominatingSets:= module()
local
N, s,
Scans:= subs~(
_s=~ [0, s, s[-1]],
proc()
local N0:= eval(N), Ns, v;
N:= evaln(N); #garbage collection
(
for s, Ns in eval(N0) do
(
for v from _s+1 to n do
if (N[s,v]:= Ns union A0[v]) = V then [s,v] fi
od
)
od
)
end proc
),
ModuleApply:= proc(\$)
local Scan, R;
N:= table([()= {}]);
for Scan in Scans do if (R:= Scan()) <> () then return [R] fi od;
[do R:= Scan() until R <> ()]
end proc
;
end module,

GetEdges:= ()-> local k; seq['reduce'= `union`](`{}`~(k, select(`>`, A0[k], k)), k= 1..n-1),

GetGraphData:= proc(G::Graph) #Assign exports specific to this graph:
VL:= op(3,G);
n:= nops(VL);
V:= {\$n};
A0:= [seq](op(4,G) union~ `{}`~(V));
edges:= GetEdges();
NE:= nops(edges);
DS:= MinDominatingSets()
end proc,

#This proc essentially makes a local copy of the modified graph,
#but does so in an efficient way:
RemoveEdges:= proc(E::{set,list}(set), \$)
local A1:= rtable(A0), e;
for e in E do A1[e[1]] minus= {e[2]}; A1[e[2]] minus= {e[1]} od;
A:= seq(A1); #This is why A is thread_local
return
end proc,

#returns true iff no dominating set of the original graph is still dominating after
#edge removal:
CheckDS:= proc() option threadsafe; local d; andseq(`union`(A[d]) <> V, d= DS) end proc,

#much slower than CheckDS, but potenially useful for some "greedy" algorithms: selects
#all of the original dominating sets that are still dominating after edge removal:
ReduceDS:= ()-> select(d-> `union`(A[d]) = V, DS),

ModuleApply:= module()
local
k,        #size of edge subsets being checked
Min, Max,    #min and max values of k to use
C,        #combinations Iterator for k-subsets
st,        #time at start of an iteration
found?,     #Flag set in one thread to make the others stop: Has a bondage set been found?

TimeReport:= proc() #per-iteration userinfo timing report and estimate:
local T:= time['real']() - st;
if k = Min or T = 0 then return fi;
userinfo(
1, 'MinBondageSet', 'NoName',
sprintf(
"Time to check %d-subsets of edges: %9.1f s; "
"expected time to check %d-subsets: %9.1f s.",
k, T, k+1, T*(NE-k)/(k+1)
)
)
end proc
;
export
CheckSplit:= proc(rnk, num) #main multithreaded procedure:
local (h,g):= ModuleIterator(Object(C, 'rank'= rnk)), J:= g(), B;
to num while h() do
RemoveEdges((B:= edges[[seq](J+~1)]));
until found? #can be set true by any thread
end proc,

ModuleApply:= proc(
G::Graph,
{ #keyword options:
#Return after fully processing the input graph but
#before examining any edge subset:
setup_only::truefalse:= false,

#Skip the setup phase because the setup of G is already stored
#in the module:
no_setup::truefalse:= false,

#sizes of edge subsets to check for being bondage sets:
set_sizes::range({posint, identical()}):= ..
},
\$
)
local ne, R;
if not no_setup then GetGraphData(G) fi;
if NE = 0 then return {} fi;
if setup_only then return fi;

Min:= lhs(set_sizes); if Min=() then Min:= 1 fi;
Max:= rhs(set_sizes); if Max=() then Max:= NE fi;
ne:= binomial(NE, Min);
found?:= false;
for k from Min to Max do
st:= time['real']();
C:= It:-RevolvingDoorCombination(NE, k);
if R <> () then return R fi;
ne*= (NE-k)/(k+1); #update binomial(NE, k) -> binomial(NE, k+1)
TimeReport()
od
end proc
;
end module,

#utility that lets the user convert vertex-label indices (default output) back
#to the original vertex labels,
LabelVertices:= e-> evalindets(e, integer[1..n], j-> VL[j])
;
end module
:```

Minimum bondage sets of graphs[1]:

A multi-threaded use of the Iterator package

Author: Carl Love <carl.j.love@gmail.com> 2023-Jun-22

[1] We only consider here so-called simple graphs: no directed edges, no self-loops, no multiple edges. In other words, a simple graph's edge set is nothing more than a subset of the set of 2-subsets of its vertices. In the following worksheet, "graph" will always mean "simple graph" without further comment. Likewise, we only consider graphs with a finite number of vertices.

Definitions:

Of course, it's expected that the reader has some familiarity with graph theory; so, the first two definitions here are primarily to clarify my notation for the reader.

1. Graph, Vertices, Edges: A graph  is an ordered pair of sets G = (V, E), where V is any nonempty set (called the vertices), and E is any set of 2-subsets of V. E is called the edges.

Henceforth, let G = (V, E) be a graph.

2. Neighborhood: Given v 2V, the (closed} neighborhood  of v is the union of all edges containing v. The neighborhood of v contains v itself. If it's necessary to make that distinction (it isn't in this worksheet), it's called a closed neighborhood. For A3V, the neighborhood of A is the union of the neighborhoods of the elements of A.

3. Dominating set: A dominating set of G is a subset of V whose neighborhood is V.

4. Minimum dominating set: A minimum dominating set of G is a dominating set of the minimum size of all its dominating sets. (There is a closely related concept called a minimAL dominating set, which is not used in this worksheet.)

5. Domination number: The  (lower) domination number of G is the size of any of its minimum dominating sets. (A similar concept applied to minimAL dominating sets is called the upper domination number. We have no need of this here, so I'll just say domination number for the lower domination number.) The domination number of G is denoted by γ(G).

6. Bondage set: A bondage set (I prefer the term unbinding set) of G is a B 3 E such that the domination number of (V, E \ B) is greater than the domination number of G.

7. Minimum bondage set: A minimum bondage set of G is a bondage set of the minimum size of all its bondage sets.

8. Bondage number: The bondage  number of G is the size of any of its minimum bondage sets.

The following proposition is fundamental to the algorithm below, making it possible to function within a reasonable amount of time and memory:

Proposition: Let S be the set of all minimum dominating sets of G. A B 3E is a bondage set of G if and only if no member of S is a dominating set of (V, E \ B).

Proof: (0) If B is a bondage set, then γ((V, E \ B)) > γ(G). Thus every dominating set of (V, E \ B) has greater than γ(G) elements. Thus no dominating set of (V, E \ B) is in S.

(*) Removing an edge from G cannot decrease its domination number. Since S is all dominating sets of G of size γ(G), γ((V, E \ B)) >  γ(G). Thus B is a bondage set.

The bondage set produced as a result of the proposition is not necessarily minimum.

The reason that this can vastly improve the efficiency of finding a bondage set is that there are likely to be relatively few minimum dominating sets of G, and they'll always be the same regardless of the B currently being considered as a bondage set. In the main example below, there are 300 minimum bondage sets (of size 4) for a graph with |V| = 24. Checking whether any of these is a dominating set of (V, E \ B) is trivial. Indeed, we can stop as soon as any of the 300 is found to be dominating, thus verifying that B is not a bondage set.

Furthermore, that checking can be done by shared-memory parallel processing via Maple's Threads:-Task model using very little memory by using the Iterator package.

 > restart :
MinBondageSet:= module()
 > # # Example: A 7-regular graph with 24 vertices # g:= (GT:= GraphTheory):-Graph(     (`union`@op@(`{}`~@op@`[]`)~)(         [\$24],         [             {2,3,4,5,7,9,11}, {3,5,6,7,12,14},  {4,7,8,9,16},  {5,9,10,11,17},             {6,11,12,19},     {7,12,13,14,21},  {8,14,15},     {9,14,15,16,22},             {10,16,17},       {11,17,18,19,23}, {12,18,19},    {13,19,20},             {14,19,20,21,24}, {15,21},          {16,21,22,24}, {17,22,23},             {18,22,23},       {19,20,23,24},    {20},          {21,23,24},             {22,24},          {23,24},          {24},          {}         ]     ) );

 > MBS:= MinBondageSet: infolevel[MBS]:= 1:
 > CodeTools:-Usage(MBS(g, setup_only));

memory used=13.97MiB, alloc change=37.00MiB, cpu time=109.00ms, real time=77.00ms, gc time=62.50ms

 > # Just out of curiosity, display and count the dominating sets: # MBS:-DS;

 > nops(%);

 > # # Remove a specific subset of 7 edges (all edges adjacent to vertex 1), and prove that the # domination number is <= 7: # MBS:-RemoveEdges(`{}`~(1, {2,3,4,5,7,9,11})); MBS:-CheckDS();

 > # # That means that none of g's original 300 minimum dominating sets is still dominating; # thus the domination number must have increased; thus the bondage number is <= 7. # # Prove that the bondage number is > 5 by exhaustive multi-threaded search: # CodeTools:-Usage(MBS(g, no_setup, set_sizes= ..5));

Time to check 2-subsets of edges:       0.6 s; expected time to check 3-subsets:      15.9 s.

Time to check 3-subsets of edges:       1.6 s; expected time to check 4-subsets:      32.8 s.

Time to check 4-subsets of edges:      29.2 s; expected time to check 5-subsets:     467.9 s.

Time to check 5-subsets of edges:     537.7 s; expected time to check 6-subsets:    7079.2 s.

memory used=182.78GiB, alloc change=267.31MiB, cpu time=70.33m, real time=9.50m, gc time=5.44m

 > # time per edge subset checked: (0.6 + 1.6 + 29.2 + 537.7)/add(binomial(84,k), k= 1..5);

 > 60/%;

 > # So, it's using about 16 microseconds per edge subset checked. # I'm using an Intel Core i7-7700HQ with 4 cores @ 2.8 GHz. There are 2 hyperthreads # per core, so a total of 8 threads in the parallel process.

I also ran MinBondageSet on the 6-subsets of edges, finding no bondage sets. This takes about 2 hours on my computer. Memory usage is trivial. This confirms that the bondage number is 7.

 >

## 1st suggestion: Remove inner circles...

Suggestions:

The circles for r = 1/3 and r = 2/3 serve no didactic purpose, and thus I think that they're detrimental to the pedagogy.

You need rays for theta = Pi/4, 3*Pi/4, etc.

The vertical lines should not extend beyond the circle. If they go outside the circle, it's easy to mistake them as stray marks. Their color should be more visible.

## No free variables...

Your series has no free variables. The concept of radius of convergence doesn't make sense without a free variable such that the convergence of the series depends on that free variable being within a certain "radius" of some point. The concept of power series requires a free variable which is raised to "powers".

## [...] can't be used in place of (...)...

In your first expression, TPP, you have square brackets used algebraically; in other words, they're used as a substitute for parentheses. While that is indeed standard in mathematical typesetting, and I agree that it improves readibility, it's not alllowed in Maple (nor in any other computer language that I'm aware of). The reason is that with the severe limitation of the number of special characters available on a standard keyboard, they need to all have their own meanings.

If I replace those square brackets with paremtheses, the error goes away.

## Embedded expression for-do loop...

In relatively recent versions of Maple, and only using 1-D input, a for-do loop (or, indeed, any do loop) can be used as an embedded expression whose output is an expression sequence:

```A:= select(isprime, {\$3..101});
A := {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
59, 61, 67, 71, 73, 79, 83, 89, 97, 101}

B:= {for a in A do if a mod 3 = 1 then a fi od};
B := {7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97}

```

This is not just a syntactic nicety; it's often the most efficient way to solve set-building problems. However, in this case, acer's select works efficiently also, and it's what I would use.

 First 14 15 16 17 18 19 20 Last Page 16 of 386
﻿