Carl Love

## 24805 Reputation

10 years, 113 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

## Prime testing is much faster than factor...

It is surprising---but nonetheless true---that determining whether or not a large number is prime can be done in much less time than factoring a number of similar size that has at least two large prime factors. That fact is fundamental to the usefulness of the RSA encryption method. So the fact in the first sentence is fairly well known because that encryption method is of immense practical importance.

Let's see just how fast it is. First, let's estimate the number of 116-digit primes:

evalf(Int(1/ln(x), x= 1e115..1e116));
3.37895657556787*10^113

That's far more than the number of atoms in the Universe or even the volume of the observable Universe in cubic nanometers. So, if we were to choose one of these numbers at random, the chance that it could've been pre-stored is infinitesimal, right?  Surely it'd take more than one atom to store a 116-digit number.

```R:= rand(10^115..10^116-1):
randprime:= proc()
local n, p:= R();
for n while not isprime(p) do p:= R() od;
n, p
end proc
:
(n,p):= CodeTools:-Usage(randprime());
memory used=0.71MiB, alloc change=0 bytes,
cpu time=0ns, real time=5.00ms, gc time=0ns

n, p := 210, 24729845748629285111963221649430106890373391188328886700762826051452963665879923038495438270815820598184513197776291
```

So, it took 5 milliseconds to generate and test 210 random 116-digit numbers, the last of which was the first prime among them. So generating and testing 116-digit numbers takes about 25 microseconds apiece.

## Need to solve the 2 ODEs together...

It's possible to solve your ode1 by itself because it doesn't contain theta(x). But it's not possible to solve ode2 by itself because it contains 2 unknown functions, f(x) and theta(x). Nor is there any need to solve it by itself; just solve them together:

b5:= subs(a5, {ode1, ode2});
c5:= dsolve({b5[], bcs, bcs1}, numeric);

And do likewise for b6 and c6.

## Using matrices...

By using matrix multiplication and matrix stacking, the code can be shortened to this:

```restart:
Digits := 15:
T5Scheme:= proc(XY0::Matrix, a::Matrix, L::nonnegint)
local i, k, N:= upperbound(XY0)[1], XY:= rtable(XY0);
for k to L do XY:= <seq(a.XY[i-2..i+2], i= 3..3^(k-1)*(N-6)+4)> od;
XY
end proc
:
## Initial polygon
q2:= -sqrt(2): q3:= sqrt(3):
XY:= <
-2,  1, 2,  0, q2, -2,  1, 2,  0, q2, -2,  1, 2,  0, q2  | #x
0, q3, 0, -2, q2,  0, q3, 0, -2, q2,  0, q3, 0, -2, q2    #y
>/2:
st2:= <-2, 1, 2, 0, q2, -2  |  0, q3, 0, -2, q2, 0>/2:
L:= [13, -22, -81, -107, 179, 729, 1233]/1296:
a:= Array(-6..8, [L[], ListTools:-Reverse(L)[], 0]):
J:= [seq](6..-6, -3):
A:= <a[J], a[J+~1], a[J+~2]>:
st1:= T5Scheme(XY, A, 3):
plot([st1, st2], linestyle= [1, 3], color= [black, red]);
```

## Use the numerators...

In almost all Maple commands, using an expression that's not an equation as an equation is equivalent to setting that expression equal to 0. Solving an expression set equal to 0 is almost equivalent to just solving its numerator set to 0, the only difference being that you might want to validate the solutions by checking whether they make the denominator 0, possibly under some restrictions. So, this works for your case:

pdsolve(numer~([eq1, eq2]));

Regarding your literal question "What does this error mean?": I don't know. It's very often much easier for me to correct or prevent an error in Maple than to understand it. It's a good question, but answering it might take several hours of investigating, digging through the code.

## Procedure for any of the 2^8 graph produ...

Using the same 1--8 numbering scheme as @lcz posted (from  Barik, S., Bapat, R.B., Pati, S.: "On the Laplacian spectra of product graphs" Appl. Anal. Discrete Math. 9, 39–58 (2015)), here is a procedure for any of the 2^8 = 256 possible products of two graphs:

```RelationGraph:= proc(f, V::list({integer, symbol, string}))
local n:= {\$1..nops(V)}, i, F:= i-> j-> f(V[i],V[j]);
GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
end proc
:
GraphProducts:= proc(
G::Graph, H::Graph,
Rel::{
set(integer[1..8]),
identical(
"Cartesian", "box", "strong", "normal", "AND", "modular", "homomorphic",
"lexicographic", "co-normal", "disjunctive", "OR",
"tensor", "Kronecker", "categorical", "direct"
)
}:= "Cartesian",
{
vertex_format::{string, symbol}:= "%A:%A",
labeltype::identical(string, symbol):= string
}
)
local
V:= GraphTheory:-Vertices~([G,H]),
E:= op~(4, [G,H]),
i, j, P:= [seq](seq([i,j], j= 1..nops(V[2])), i= 1..nops(V[1])),
Vgh:= (
curry(`if`(labeltype=string, sprintf, nprintf), vertex_format)
@ ((i,j)-> (V[1][i],V[2][j])) @ op
)~(P),
J:= op~(table(Vgh=~ P)),
Edge:= (u,v,k)-> J[v][k] in E[k][J[u][k]],
Eq:= (u,v,k)-> J[v][k] = J[u][k],
Rs:= table([
1= ((u,v)->     Edge(u,v,1) and       Eq(u,v,2)),
2= ((u,v)-> not Edge(u,v,1) and       Eq(u,v,2)),
3= ((u,v)->     Edge(u,v,1) and     Edge(u,v,2)),
4= ((u,v)-> not Edge(u,v,1) and     Edge(u,v,2)),
5= ((u,v)->     Edge(u,v,1) and not Edge(u,v,2)),
6= ((u,v)-> not Edge(u,v,1) and not Edge(u,v,2)),
7= ((u,v)->       Eq(u,v,1) and     Edge(u,v,2)),
8= ((u,v)->       Eq(u,v,1) and not Edge(u,v,2))
]),
Rels:= table([
("Cartesian", "box")=~ ``({1,7}),
("strong", "normal", "AND")=~ ``({1,3,7}),
"lexicographic"= ``({1,3,5,7}),
"modular"= ``({3,6}),
"homomorphic"= ``({5,7,8}),
("co-normal", "disjunctive", "OR")=~ ``({1,3,4,5,7}),
("tensor", "Kronecker", "categorical", "direct")=~ ``({3})
])
;
RelationGraph(
subs(_R= op(Rels[Rel]), (u,v)-> u<>v and ormap(k-> Rs[k](u,v), _R)),
Vgh
)
end proc
:
#Usage:
G:= GraphTheory:-CycleGraph(5): H:= GraphTheory:-CycleGraph(4):
GH:= GraphProducts(G, H, "co-normal", labeltype= symbol);
GH := Graph 3: an undirected unweighted graph with 20 vertices and 140 edge(s
```

The third argument can be any subset of {1,2,3,4,5,6,7,8}; seven of the subsets are given string names (with some synonyms) for convenience. The labeltype of the product's vetrtices defaults to string.

## ||...

What you want can be done with the || operator, which is very similar to cat, but without the evaluation issue that's causing your problem.

x||(a,b,c):= y =~ x^~(1,2,3);
rhs~([x||(a,c)])[];

## Can't be done...

Of course it can't be done. For it to be possible at all, at least the positions of the 0s must be the same.

## Plotting the numbers...

Here is a procedure that plots the numbers themselves on the triangle's sides with the text orientation parallel to sides.

```restart
:
RandPythTriple:= proc({upper::posint:= 2^62})
local
ub:= isqrt(iquo(upper, 2)),
n:= rand(1..iquo(ub,2))(),
m:= rand(1+n..ub)()
;
(min,max)(m^2-n^2, 2*m*n), m^2+n^2
end proc
:
PlotTriple:= proc(a::posint, b::posint, c::posint)
uses P= plots, PT= plottools;
local prim:= igcd(a,b,c)=1, A:= a/c, B:= b/c;
P:-display(
PT:-polygon([[0,0], [0,A], [B,0]], 'color'= `if`(prim, 'red', 'blue')),
P:-textplot(
[
[0, A/2, a, 'rotation'= Pi/2, 'align'= 'left'],
[B/2, 0, b, 'align'= 'below'],
[B/2, 1.1*A/2, c, 'rotation'= -arctan(a/b)]
],
'font'= ['times', 'bold', 16]
),
'title'= cat(`if`(prim, "Primitive ", ""), "Pythagorean Triple"),
'titlefont'= ['times', 'bold', 20],
'view'= [0..B, 0..A], 'axes'= 'none', 'scaling'= 'constrained', _rest
)
end proc
:
PlotTriple(RandPythTriple(upper= 99));
```

## simplify(..., zero)...

Like this:

simplify(ans, zero);

## Double-underscore subscripts...

Like this:

a__01 #Two underscores

## [seq]...

To return all the values of w, replace your for loop (all 3 lines) with this:

w:= [seq](gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1), i= 1..m);

That makes w a list. What you had makes w a table. There's nothing wrong with using a table, but they're slightly more difficult to return.

## General procedure for all such questions...

As promised, here's a general procedure for all such questions involving an equivalence relation between the rows of two matrices.

 > restart :
 > A:= <     1, 3, 3, 2, 3, 2, 2, 3, 3, 2, 0;     1, 3, 3, 2, 3, 1, 1, 2, 3, 3, 2;     1, 3, 2, 2, 2, 3, 1, 2, 3, 3, 2;     1, 3, 5, 3, 2, 2, 3, 3, 2, 0, 0;     1, 3, 5, 3, 1, 1, 2, 3, 3, 2, 0;     1, 3, 4, 2, 3, 1, 2, 3, 3, 2, 0;     1, 3, 6, 4, 2, 3, 3, 2, 0, 0, 0;     1, 3, 6, 3, 1, 2, 3, 3, 2, 0, 0;     1, 3, 5, 5, 5, 3, 2, 0, 0, 0, 0;     1, 3, 3, 3, 5, 4, 3, 2, 0, 0, 0;     1, 3, 4, 4, 4, 3, 3, 2, 0, 0, 0;     1, 3, 3, 4, 6, 5, 2, 0, 0, 0, 0;     1, 3, 5, 5, 6, 4, 0, 0, 0, 0, 0;     1, 3, 5, 6, 4, 3, 2, 0, 0, 0, 0;     1, 3, 3, 5, 5, 2, 3, 2, 0, 0, 0;     1, 3, 5, 5, 3, 2, 3, 2, 0, 0, 0;     1, 3, 3, 4, 3, 3, 2, 3, 2, 0, 0;     1, 3, 4, 4, 2, 3, 2, 3, 2, 0, 0;     1, 3, 5, 2, 1, 2, 3, 2, 3, 2, 0;     1, 3, 5, 3, 2, 3, 2, 3, 2, 0, 0;     1, 3, 4, 3, 1, 2, 3, 2, 3, 2, 0;     1, 3, 3, 3, 2, 2, 3, 2, 3, 2, 0;     1, 3, 3, 3, 1, 1, 2, 3, 2, 3, 2;     1, 3, 2, 2, 3, 1, 2, 3, 2, 3, 2 >:

For testing, I construct matrix B by permuting (at random) a subset of the rows of A and then doing a random cyclic permutation on the interior columns.

 > (r,c):= op(1,A);

 > C:= combinat:
 > cc:= rand(2..c-1)();

 > B:= A[C:-randperm(r)[..rand(1..r)()], [1, cc..-2, 2..cc-1, c]];

 > CycPermRep:= proc(r::posint, M::Matrix) description " This procedure returns (as a list) the lexicographically least cyclic permutation of M[r, 2..-2]. The time complexity is O(c) where c is the number of columns of M. "; local A:= [seq](M[r, 2..-2]), p;     {seq}(A[[p..-1, 1..p-1]], p= ListTools:-SearchAll({A[]}[1], A))[1] end proc :
 > MatrixRowMatch:= proc(A::Matrix, B::Matrix, ClRep) description " Author: Carl Love 2022-Aug-28 . For matrices A and B, this procedure produces a bipartite graph showing a relation between the rows of A and those of B. The relation is determined by the 3rd argument, procedure or operator ClRep. . The matrices can be any sizes; if one has fewer columns than the other, it'll be padded on the right with 0s. Then the matrices are stacked into a single matrix AB. ClRep's arguments will be a row number i of AB, AB itself, and any extra arguments passed to this procedure. It should return a canonical identifier (of any immutable type) of the relation's equivalence class that contains row AB[i]. (The return value should be immutable, such as a list; thus, it shouldn't be a vector.) . The equivalence classes are used to form a bipartite graph by removing connections from A to itself and B to itself. The graph's vertices are [cat(\"A\", 1..ra), cat(\"B\", 1..rb)] where ra and rb are the numbers of rows of A and B. . The time complexity of this algorithm is O(r*C) where r is the number of rows of AB and C is the time complexity of C1Rep. "; local AB, Cl, ra, ca, rb, cb, r, RepRem, rest:= _rest;     (ra, ca):= op(1, A); (rb, cb):= op(1, B); r:= ra+rb;     AB:= <, >;     RepRem:= proc(i::posint) option remember; ClRep(i, AB, rest) end proc;     Cl:= ListTools:-Classify(RepRem, [\$1..r]);     GraphTheory:-Graph(         [cat("A", 1..ra), cat("B", 1..rb)],         Array(1..r, i-> remove(`if`(i > ra, `>`, `<=`), Cl[RepRem(i)], ra))     ) end proc :
 > BG:= MatrixRowMatch(A, B, CycPermRep);

 > GT:= GraphTheory:
 > E:= GT:-Edges(BG);

 > for e in E do     print(A[parse(e[1][2..])], B[parse(e[2][2..])]) od;

 > GT:-DrawGraph(BG, layout= bipartite);

 >

## Procedure for it...

Here is a procedure for it (based on my limited understanding from a very quick lookup of the definition of the co-normal product of graphs). Let me know if I got it right. The definition that I used is that conormal(G,H) contains the edge {u1:v1, u2:v2} iff (G contains {u1, u2} or H contains {v1, v2}).

```restart
:
RelationGraph:= (f, V::list({integer, symbol, string}))->
local n:= {\$nops(V)}, i, F:= i-> j-> f(V[i],V[j]);
GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
:
ConormalProduct:= proc(
G::Graph, H::Graph, {vertex_format::{string, symbol}:= "%a:%a"}
)
local Vg, Vh, ng, nh, Ng, Nh, P, i, j, Vgh, K;
(Vg,Vh):= op(GraphTheory:-Vertices~([G,H]));
(ng,nh):= op(`..`~(1, nops~([Vg,Vh])));
(Ng,Nh):= op(op~(4, [G,H]));
P:= [seq](seq([i,j], j= nh), i= ng);
Vgh:= (curry(sprintf, vertex_format)@((i,j)-> (Vg[i],Vh[j]))@op)~(P);
K:= op~(table(Vgh=~ P));
RelationGraph(
proc(a, b, \$)
local u, v, x, y;
(u,v):= K[a]; (x,y):= K[b];
x in Ng[u] or y in Nh[v]
end proc,
Vgh
)
end proc
:
```

## I've had many problems with coerce...

I've had many problems using coerce over the years, although it seems much better in Maple 2022 than previously. That's why for a previous Question from you, I recommended using coercion procedures whose names begin rather than using the coerce parameter modifier.

I have no experience using Maple's debugger (or any other debugger for that matter), but as far as I can tell, the following works with the debugger and correctly does the coercion. Let me know if it works for you.

```restart;
~list:= (e, T::type)->
if e::{set, thistype}(T) then [`if`](e::set, e[], e)
else error "%1 not cercible to type list(%2)", A, T
fi
:
A:= module()
option object;
export
ode, ic,

ModuleCopy::static:= proc(_self, proto::A, ode::`=`, {ic::~list(`=`):= []}, \$)
print("ode=",ode);
print("ic=",ic)
end proc
;
end module;

foo:= proc(ode::`=`, {ic::~list(`=`):= []}, \$)
local o;
#DEBUG();
o:= Object(A, ode, ':-ic'= ic)
end proc:

foo(diff(y(x),x)=1); ```