Applications, Examples and Libraries

Share your work here

Seeking for fast approximate formulas to compute (a huge number of) quantiles of a Gaussian random variable (here the standard one, but its extension to any Gaussian RV is straightforward), I found a few of them in the Abramowitz and Stegun book, page 933, relations 26.2.22 and 26.2.23.
Each approximation model is expressed as a rational fraction, the second one being the more accurate.
Each model depends on (respectively 4 and 6) parameters that are estimated (I guess it was done this way) through a least-square-like method.

See here for an online access http://people.math.sfu.ca/~cbm/aands/page_933.htm.

These approximation, and specially the most accurate one (formula 26.2.23) seem to be still widely used today(1) (see for instance https://www.johndcook.com/blog/normal_cdf_inverse/ ).

As an amusement I decided to compute the best fit by using the Statistics:-NonLinearFit procedure and a sample of (probability, quantile) points where probability ranges in [0.5, 1-1/1000] (the range used in formulas 26.2.22 and 26.2.23 is (0, 0.5] but this is not a point).
Surprisingly Statistics:-NonLinearFit returned, for the two formulas, parameter estimations substantially different from the one given in the Abramowitz & Stegun's book. A reason could be that the points I used when I did the fits weren't the one they used (unfortunately they give no informations about this).

More interesting, whatever the formula I refitted,  NonLinearFit produced an approximation whose the absolute error was smaller by about two orders of magnitude to the onesprovided by Abramowitz and Stegun.
For instance they wrote that the most accurate formula (26.2.23) had an absolute approximation error less than 4.5*10-4 as I obtained a value around 10-6!

(1) To get an idea of the persistence of the use of the formula 26.2.23, just type the value 2.515517 of its parameter c[0] in any search engine.


In the plots below the gray rectangle refers to the region where the approximate ICDF is used for extrapolation.
 

restart:

with(Statistics):

cdf := unapply(evalf(CDF(Normal(0, 1), x)), x):
X   := [seq(0..5, 0.1)]:
A   := cdf~(X):
T  := alpha -> sqrt(-2*log(1-alpha)):
q  := Quantile~(Normal(0, 1), A):
Aq := convert([A,q], Matrix)^+:

r := 1:

J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):


model  := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);


# these lines are for estimating the performances
B  := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------
Y  := [seq(0..6, 0.01)]:
B  := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.5454311687345044)+HFloat(0.8058592540791468)*(-2*ln(1-alpha))^(1/2))/(1+HFloat(1.4689746699940707)*(-2*ln(1-alpha))^(1/2)-HFloat(0.34455942407858625)*ln(1-alpha)) end proc

 

memory used=170.31MiB, alloc change=76.01MiB, cpu time=3.06s, real time=3.05s, gc time=54.87ms

memory used=171.59MiB, alloc change=256.00MiB, cpu time=3.12s, real time=3.03s, gc time=154.77ms

memory used=8.24MiB, alloc change=0 bytes, cpu time=95.00ms, real time=95.00ms, gc time=0ns

 

 

r := 2:

 
J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):


model  := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);


# these lines are for estimating the performances
B  := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------


Y  := [seq(0..6, 0.01)]:
B  := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.9637294443959394)+HFloat(4.527738737327481)*(-2*ln(1-alpha))^(1/2)-HFloat(0.9571637188191973)*ln(1-alpha))/(1+HFloat(3.472400103322335)*(-2*ln(1-alpha))^(1/2)-HFloat(3.426536241250657)*ln(1-alpha)+HFloat(0.08875278252087411)*(-2*ln(1-alpha))^(3/2)) end proc

 

memory used=170.09MiB, alloc change=32.00MiB, cpu time=3.29s, real time=3.11s, gc time=268.60ms

memory used=170.85MiB, alloc change=0 bytes, cpu time=3.23s, real time=3.10s, gc time=201.52ms
memory used=10.76MiB, alloc change=0 bytes, cpu time=127.00ms, real time=127.00ms, gc time=0ns

 

 

# Optimized "r=2" computation

z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
z_fit := unapply(convert~(%, horner), z);

p := proc(alpha)
  local z:
  z := sqrt(-2*log(1-alpha)):
  z_fit(z):
end proc:

R3 := CodeTools:-Usage(p~(B)):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (z) options operator, arrow; (-2.963729444+(-3.527738737+(2.993818244+(1.713268121+0.8875278252e-1*z)*z)*z)*z)/(1.+(3.472400103+(1.713268121+0.8875278252e-1*z)*z)*z) end proc

 

memory used=1.67MiB, alloc change=0 bytes, cpu time=14.00ms, real time=15.00ms, gc time=0ns

 

 


AS stands for Abramowith & Stegun

J_AS := unapply(normal(eval(J(t), [a__0=2.515517, a__1=0.802853, a__2=0.010328, b__1=1.432788, b__2=0.189269, b__3=0.001308])), t):
J_AS(t);


# for comparison:

print():
z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
map(sort, %, z);

plot([z_fit(z), J_AS(z)], z=0.5..1, color=[blue, red], legend=[mmcdara, Abramowitz_Stegun], gridlines=true);

print():
R2_AS := CodeTools:-Usage(J_AS~(T~(B))):
print():


plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2_AS-~R1)), legend=Abramowitz_Stegun, gridlines=true, size=[700, 400]),
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

(0.1308000000e-2*t^4+.1892690000*t^3+1.422460000*t^2+.1971470000*t-2.515517000)/(0.1308000000e-2*t^3+.1892690000*t^2+1.432788000*t+1.)

 

 

(0.8875278252e-1*z^4+1.713268121*z^3+2.993818244*z^2-3.527738737*z-2.963729444)/(0.8875278252e-1*z^3+1.713268121*z^2+3.472400103*z+1.)

 

 

 

memory used=2.92MiB, alloc change=0 bytes, cpu time=25.00ms, real time=25.00ms, gc time=0ns

 

 

 


 

Download InverseNormalCDF.mw

 

 

Although the graph of a parametrized surface can be viewed and manipulated on the computer screen as a surface in 3D, it is not quite suitable for printing on a 3D printer since such a surface has zero thickness, and thus it does not correspond to physical object.

To produce a 3D printout of a surface, it needs to be endowed with some "thickness".  To do that, we move every point from the surface in the direction of that point's nomral vector by the amount ±T/2, where T is the desired thickness.  The locus of the points thus obtained forms a thin shell of thickness T around the original surface, thus making it into a proper solid. The result then may be saved into a file in the STL format and be sent to a 3D printner for reproduction.

The worksheet attached to this post provides a facility for translating a parametrized surface into an STL file.  It also provides a command for viewing the thickened object on the screen.  The details are documented within that worksheet.

Here are a few samples.  Each sample is shown twice—one as it appears within Maple, and another as viewed by loading the STL file into MeshLab which is a free mesh viewing/manipulation software.

 

Here is the worksheet that produced these:  thicken.mw

 

 

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])<abs(i-i1)+abs(j-j1) then return `no solutions` fi;
od; od; od; od;
L:=sort(select(e->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]<MS,seq(L[-1]+k,k=0..MS-L[-1]-1),NULL)];
OneStepLeft:=proc(A1::listlist)
local s, M, m, k, T;
uses ListTools;
s:=Search(a, 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 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

[Matrix(3, 3, {(1, 1) = 1, (1, 2) = 6, (1, 3) = 7, (2, 1) = 2, (2, 2) = 5, (2, 3) = 8, (3, 1) = 3, (3, 2) = 4, (3, 3) = 9}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 8, (1, 3) = 7, (2, 1) = 2, (2, 2) = 9, (2, 3) = 6, (3, 1) = 3, (3, 2) = 4, (3, 3) = 5}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 8, (1, 3) = 9, (2, 1) = 2, (2, 2) = 7, (2, 3) = 6, (3, 1) = 3, (3, 2) = 4, (3, 3) = 5}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 4, (1, 3) = 5, (2, 1) = 2, (2, 2) = 3, (2, 3) = 6, (3, 1) = 9, (3, 2) = 8, (3, 3) = 7}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 9, (2, 1) = 4, (2, 2) = 3, (2, 3) = 8, (3, 1) = 5, (3, 2) = 6, (3, 3) = 7}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 8, (2, 2) = 7, (2, 3) = 4, (3, 1) = 9, (3, 2) = 6, (3, 3) = 5}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 8, (2, 2) = 9, (2, 3) = 4, (3, 1) = 7, (3, 2) = 6, (3, 3) = 5}), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 6, (2, 2) = 5, (2, 3) = 4, (3, 1) = 7, (3, 2) = 8, (3, 3) = 9})]

 

[]

 

52

 

25

 

36

(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);

12440

(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:

 

Download SerpPathsinMatrix.mw

 

I experienced a significant obstacle while trying to generate independent random samples with Statistics:-Sample on different nodes of a Grid multi-processing environment. After many hours of trial-and-error, I discovered an astonishing workaround, and I achieved excellent time and memory performance. Since this seems like a generally useful computation, I thought that it was worthy of a Post.

This Post is also worth reading to learn how to use Grid when you need to initialize a substantial environment on each node before using Grid:-Map or Grid:-Seq.

All remaining details are in the following worksheet.
 

How to use Statistics:-Sample in the `Grid` environment

Author: Carl Love <carl.j.love@gmail.com> 1 August 2019

 

I experienced a significant obstacle while trying to generate indenpendent random samples with Statistics:-Sample on the nodes of a multi-processor Grid (on a single computer). After several hours of trial-and-error, I discovered that two things are necessary to do this:

1. 

The random number generator needs to be seeded differently in each node. (The reason for this is easy to understand.)

2. 

The random variables generated by Statistics:-RandomVariable need to have different names in each node. This one is mind-boggling to me. Afterall, each node has its own kernel and so its own memory It's as if the names of random variables are stored in a disk file which all kernels access. And also the generator has been seeded differently in each node.

 

Once these things were done, the time and memory performance of the computation were excellent.

restart
:

Digits:= 15
:

#Specify the size of the computation:
(n1,n2,n3):= (100, 100, 1000):
# n1 = size of each random sample;
# n2 = number of samples in a batch;
# n3 = number of batches.

#
#Procedure to initialize needed globals on each node:
Init:= proc(n::posint)
local node:= Grid:-MyNode();
   #This is wrapped in parse so that it'll act globally. Otherwise, an environment
   #variable would be reset when this procedure ends.
   parse("Digits:= 15;", 'statement');

   randomize(randomize()+node); #Initialize independent RNG for this node.
   #If repeatability of results is desired, remove the inner randomize().

   (:-X,:-Y):= Array(1..n, 'datatype'= 'hfloat') $ 2;

   #Perhaps due to some oversight in the design of Statistics, it seems necessary that
   #r.v.s in different nodes **need different names** in order to be independent:
   N||node:= Statistics:-RandomVariable('Normal'(0,1));
   :-TRS:= (X::rtable)-> Statistics:-Sample(N||node, X);
   #To verify that different names are needed, change N||node to N in both lines.
   #Doing so, each node will generate identical samples!

   #Perform some computation. For the pedagogical purpose of this worksheet, all that
   #matters is that it's some numeric computation on some Arrays of random Samples.
   :-GG:= (X::Array, Y::Array)->
      evalhf(
         proc(X::Array, Y::Array, n::posint)
         local s, k, S:= 0, p:= 2*Pi;
            for k to n do
               s:= sin(p*X[k]);  
               S:= S + X[k]^2*cos(p*Y[k])/sqrt(2-sin(s)) + Y[k]^2*s
            od
         end proc
         (X, Y, n)
      )      
   ;
   #Perform a batch of the above computations, and somehow numerically consolidate the
   #results. Once again, pedagogically it doesn't matter how they're consolidated.  
   :-TRX1:= (n::posint)-> add(GG(TRS(X), TRS(Y)), 1..n);
   
   #It doesn't matter much what's returned. Returning `node` lets us verify that we're
   #actually running this on a grid.
   return node
end proc
:

The procedure Init above uses the :- syntax to set variables globally for each node. The variables set are X, Y, N||node, TRS, GG, and TRX1. Names constructed by concatenation, such as N||node, are always global, so :- isn't needed for those.

#
#Time the initialization:
st:= time[real]():
   #Send Init to each node, but don't run it yet:
   Grid:-Set(Init)
   ;
   #Run Init on each node:
   Nodes:= Grid:-Run(Init, [n1], 'wait');
time__init_Grid:= time[real]() - st;

Array(%id = 18446745861500764518)

1.109

The only purpose of array Nodes is that it lets us count the nodes, and it lets us verify that Grid:-MyNode() returned a different value on each node.

num_nodes:= numelems(Nodes);

8

#Time the actual execution:
st:= time[real]():
   R1:= [Grid:-Seq['tasksize'= iquo(n3, num_nodes)](TRX1(k), k= [n2 $ n3])]:
time__run_Grid:= time[real]() - st

4.440

#Just for comparison, run it sequentially:
st:= time[real]():
   Init(n1):
time__init_noGrid:= time[real]() - st;

st:= time[real]():
   R2:= [seq(TRX1(k), k= [n2 $ n3])]:
time__run_noGrid:= time[real]() - st;

0.16e-1

24.483

R1 and R2 will be different because different random numbers were used, but they should have similar histograms.

plots:-display(
   Statistics:-Histogram~(
      <R1 | R2>, #side-by-side plots
      'title'=~ <<"With Grid\n"> | <"Without Grid\n">>,
      'gridlines'= false
   )
);

(Plot output deleted because MaplePrimes cannot handle side-by-side plots!)

They look similar enough to me!

 

Let's try to quantify the benefit of using Grid:

speedup_factor:= time__run_noGrid / time__run_Grid;

5.36319824753560

Express that as a fraction of the theoretical maximum speedup:

efficiency:= speedup_factor / num_nodes;

.670399780941950

I think that that's really good!

 

The memory usage of this code is insignificant, which can be verified from an external memory monitor such as Winodws Task Manager. It's just a little bit more than that needed to start a kernel on each node. It's also possible to measure the memory usage programmatically. Doing so for a Grid:-Seq computation is a little bit beyond the scope of this worksheet.

 


 

Download GridRandSample.mw

Here are the histograms:

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])<abs(i-i1)+abs(j-j1) then return `no solutions` fi;
od; od; od; od;
L:=sort(select(e->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]<MS,seq(L[-1]+k,k=0..MS-L[-1]-1),NULL)];
  

OneStepLeft:=proc(A1::listlist)
local s, M, m, k, T;
uses ListTools;
s:=Search(a, 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 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(i<m, x[i+1, j, k+1], NULL)
                                                       , ifelse(j>1, x[i, j-1, k+1], NULL)
                                                       , ifelse(j<n, x[i, j+1, k+1], NULL)
                                                      )
                                 , i = 1..m)
                             , j = 1..n)
                         , k = 1 .. m*n-1));

    one := &or(seq(seq(x[i,j,1], i=1..m), j=1..n));   


    single := &not(&or(seq(seq(seq(seq(x[i,j,k] &and x[i,j,kk], kk = k+1..m*n), k = 1..m*n-1)
                                , i = 1..m), j = 1..n)));

    sol := Logic:-Satisfy(&and(initial, adjacent, one, single));
    
    if sol = NULL then
        error "no solution";
    end if;
if inline then
        S := M;
     else
        S := Matrix(m,n);
    end if;

    for eq in sol do
        (var, val) := op(eq);
        if val then
            S[op(1..2, var)] := op(3,var);
        end if;
    end do;
    S;
end proc:

           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);

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 5, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 9})

 

Matrix(3, 3, {(1, 1) = 3, (1, 2) = 4, (1, 3) = 5, (2, 1) = 2, (2, 2) = 7, (2, 3) = 6, (3, 1) = 1, (3, 2) = 8, (3, 3) = 9})

 

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 5, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 8, (3, 3) = 0})

 

Matrix(%id = 18446746210121682686), Matrix(%id = 18446746210121682806), Matrix(%id = 18446746210121674750), Matrix(%id = 18446746210121674870)

(1)


Comparison with Numbrix procedure. The example is taken from
http://rosettacode.org/wiki/Solve_a_Numbrix_puzzle 

 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));

Matrix(9, 9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 46, (2, 4) = 45, (2, 5) = 0, (2, 6) = 55, (2, 7) = 74, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 38, (3, 3) = 0, (3, 4) = 0, (3, 5) = 43, (3, 6) = 0, (3, 7) = 0, (3, 8) = 78, (3, 9) = 0, (4, 1) = 0, (4, 2) = 35, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 71, (4, 9) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 33, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 59, (5, 8) = 0, (5, 9) = 0, (6, 1) = 0, (6, 2) = 17, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 67, (6, 9) = 0, (7, 1) = 0, (7, 2) = 18, (7, 3) = 0, (7, 4) = 0, (7, 5) = 11, (7, 6) = 0, (7, 7) = 0, (7, 8) = 64, (7, 9) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 24, (8, 4) = 21, (8, 5) = 0, (8, 6) = 1, (8, 7) = 2, (8, 8) = 0, (8, 9) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 0})

 

memory used=7.85MiB, alloc change=-3.01MiB, cpu time=172.00ms, real time=212.00ms, gc time=93.75ms

 

Matrix(9, 9, {(1, 1) = 49, (1, 2) = 50, (1, 3) = 51, (1, 4) = 52, (1, 5) = 53, (1, 6) = 54, (1, 7) = 75, (1, 8) = 76, (1, 9) = 81, (2, 1) = 48, (2, 2) = 47, (2, 3) = 46, (2, 4) = 45, (2, 5) = 44, (2, 6) = 55, (2, 7) = 74, (2, 8) = 77, (2, 9) = 80, (3, 1) = 37, (3, 2) = 38, (3, 3) = 39, (3, 4) = 40, (3, 5) = 43, (3, 6) = 56, (3, 7) = 73, (3, 8) = 78, (3, 9) = 79, (4, 1) = 36, (4, 2) = 35, (4, 3) = 34, (4, 4) = 41, (4, 5) = 42, (4, 6) = 57, (4, 7) = 72, (4, 8) = 71, (4, 9) = 70, (5, 1) = 31, (5, 2) = 32, (5, 3) = 33, (5, 4) = 14, (5, 5) = 13, (5, 6) = 58, (5, 7) = 59, (5, 8) = 68, (5, 9) = 69, (6, 1) = 30, (6, 2) = 17, (6, 3) = 16, (6, 4) = 15, (6, 5) = 12, (6, 6) = 61, (6, 7) = 60, (6, 8) = 67, (6, 9) = 66, (7, 1) = 29, (7, 2) = 18, (7, 3) = 19, (7, 4) = 20, (7, 5) = 11, (7, 6) = 62, (7, 7) = 63, (7, 8) = 64, (7, 9) = 65, (8, 1) = 28, (8, 2) = 25, (8, 3) = 24, (8, 4) = 21, (8, 5) = 10, (8, 6) = 1, (8, 7) = 2, (8, 8) = 3, (8, 9) = 4, (9, 1) = 27, (9, 2) = 26, (9, 3) = 23, (9, 4) = 22, (9, 5) = 9, (9, 6) = 8, (9, 7) = 7, (9, 8) = 6, (9, 9) = 5})

 

memory used=1.21GiB, alloc change=307.02MiB, cpu time=37.00s, real time=31.88s, gc time=9.30s

 

Matrix(%id = 18446746210094669942)

(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)):

Matrix(5, 5, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 25})

 

memory used=0.94GiB, alloc change=-22.96MiB, cpu time=12.72s, real time=11.42s, gc time=2.28s

 

104

 

memory used=34.74MiB, alloc change=0 bytes, cpu time=781.00ms, real time=783.00ms, gc time=0ns

 

 


 

Download NumbrixPuzzle.mw

For no particular reason at all, these are parametric equations that print "Maplesoft" in handwritten cursive script when plotted

restart:
X := -2.05*sin(-2.70 + 2.45*t) - 3.36*sin(1.12 + 1.43*t) - 4.82*sin(-2.19 + 2.03*t) - 2.02*sin(1.36 + 2.31*t) - 2.41*sin(1.08 + 2.59*t) - 14.2*sin(1.51 + 0.185*t) - 5.25*sin(-2.04 + 1.85*t) - 2.81*sin(0.984 + 2.36*t) - 3.01*sin(-2.04 + 1.80*t) - 1.80*sin(-2.61 + 2.73*t) - 0.712*sin(-3.94 + 1.89*t) - 6.90*sin(-1.90 + 1.52*t) - 0.600*sin(-3.39 + 2.26*t) - 0.631*sin(-4.65 + 2.68*t) - 3.10*sin(-2.22 + 2.17*t) - 2.95*sin(1.38 + 1.25*t) - 1.43*sin(0.383 + 2.40*t) - 8.25*sin(-1.66 + 0.323*t) - 1.39*sin(-3.08 + 2.63*t) - 0.743*sin(-2.43 + 0.647*t) - 6.25*sin(-1.73 + 0.832*t) - 273.*sin(-1.58 + 0.0462*t) - 4.58*sin(-2.00 + 1.48*t) - 5.70*sin(-1.80 + 1.20*t) - 2.30*sin(1.42 + 0.462*t) - 3.24*sin(1.51 + 0.277*t) - 16.0*sin(-1.64 + 0.231*t) - 1.58*sin(0.779 + 1.71*t) - 0.571*sin(-2.08 + 0.970*t) - 8.85*sin(-1.88 + 1.34*t) - 1.10*sin(-2.24 + 2.08*t) - 1.49*sin(-2.27 + 1.02*t) - 2.19*sin(-1.70 + 1.94*t) - 4.47*sin(-2.06 + 1.57*t) - 2.08*sin(-2.02 + 1.06*t) - 5.70*sin(-1.86 + 1.62*t) - 2.26*sin(-1.66 + 1.16*t) - 3.95*sin(-1.98 + 1.29*t) - 0.928*sin(-2.08 + 1.76*t) - 2.98*sin(1.36 + 1.11*t) - 0.390*sin(-2.33 + 2.22*t) - 3.81*sin(1.01 + 2.54*t) - 0.613*sin(-1.43 + 1.66*t) - 19.7*sin(-1.60 + 0.138*t) - 0.524*sin(-2.87 + 0.414*t) - 2.15*sin(-4.63 + 0.694*t) - 0.782*sin(-1.56 + 2.49*t) - 5.27*sin(-1.81 + 1.38*t) - 5.18*sin(1.51 + 0.0923*t) - 6.83*sin(1.37 + 0.923*t) - 0.814*sin(-1.72 + 0.600*t) - 2.98*sin(-1.82 + 0.738*t) - 5.49*sin(1.44 + 0.509*t) - 3.90*sin(-1.76 + 0.785*t) - 0.546*sin(-2.18 + 0.876*t) - 1.92*sin(0.755 + 1.98*t) - 8.16*sin(1.38 + 0.553*t) - 0.504*sin(-1.56 + 0.371*t) - 3.43*sin(1.14 + 2.12*t):
Y := -1.05*sin(-3.81 + 2.68*t) - 7.72*sin(-4.59 + 0.231*t) - 6.38*sin(1.37 + 1.11*t) - 4.24*sin(-2.36 + 2.31*t) - 7.06*sin(1.18 + 1.80*t) - 4.60*sin(1.28 + 2.03*t) - 0.626*sin(-0.285 + 2.45*t) - 0.738*sin(-1.89 + 2.26*t) - 1.45*sin(-1.73 + 1.57*t) - 2.30*sin(-4.51 + 2.59*t) - 9.58*sin(-2.07 + 1.71*t) - 0.792*sin(-0.578 + 0.647*t) - 4.55*sin(1.49 + 1.25*t) - 14.0*sin(-2.13 + 1.62*t) - 1.02*sin(0.410 + 0.277*t) - 19.2*sin(-1.54 + 0.0462*t) - 17.3*sin(-1.86 + 1.20*t) - 1.96*sin(-0.845 + 2.63*t) - 0.754*sin(-0.0904 + 2.73*t) - 4.74*sin(1.11 + 1.48*t) - 1.79*sin(0.860 + 2.17*t) - 25.2*sin(-1.77 + 0.832*t) - 3.88*sin(1.30 + 0.462*t) - 20.8*sin(-1.66 + 0.323*t) - 17.6*sin(1.20 + 1.29*t) - 4.83*sin(0.169 + 2.36*t) - 10.8*sin(-2.01 + 1.85*t) - 8.69*sin(-2.17 + 2.22*t) - 5.48*sin(-1.69 + 1.34*t) - 18.1*sin(1.18 + 1.43*t) - 4.71*sin(0.728 + 2.08*t) - 1.15*sin(-3.44 + 1.52*t) - 2.53*sin(-2.61 + 2.54*t) - 5.48*sin(-2.02 + 1.94*t) - 4.67*sin(1.30 + 1.66*t) - 9.10*sin(1.37 + 0.970*t) - 6.45*sin(1.31 + 1.02*t) - 5.18*sin(-2.09 + 1.76*t) - 18.3*sin(-1.77 + 1.06*t) - 27.3*sin(1.31 + 1.16*t) - 2.83*sin(-3.01 + 2.40*t) - 2.93*sin(-1.70 + 0.138*t) - 4.17*sin(-2.06 + 2.12*t) - 1.60*sin(-4.25 + 1.38*t) - 2.69*sin(-1.89 + 0.371*t) - 7.92*sin(-1.78 + 0.600*t) - 19.6*sin(-1.79 + 0.738*t) - 22.6*sin(1.48 + 0.509*t) - 13.5*sin(1.21 + 0.923*t) - 5.53*sin(-1.64 + 0.0923*t) - 1.20*sin(0.145 + 2.49*t) - 3.15*sin(-1.57 + 0.414*t) - 1.74*sin(0.655 + 1.98*t) - 3.98*sin(-2.14 + 0.876*t) - 11.3*sin(-1.82 + 0.694*t) - 10.4*sin(0.987 + 1.89*t) - 8.39*sin(-1.53 + 0.185*t) - 27.8*sin(-1.76 + 0.785*t) - 9.39*sin(1.38 + 0.553*t):
plot([X, Y, t = 0 .. 68], scaling = constrained, axes = boxed);

Hare in the forest

The rocket flies

  

Быльнов_raketa_letit.mws

 

Plotting the function of a complex variable

Plotting_the_function_of_a_complex_variable.mws

 

Animated 3-D cascade of dolls

 

3d_matryoshkas_en.mws

 

With this application developed entirely in Maple using native syntax and embedded components for science and engineering students. Just replace your data and you're done.

Pearson_Coeficient.mw

Lenin Araujo Castillo

Ambassador of Maple

 

Foucault’s Pendulum Exploration Using MAPLE18

https://www.ias.ac.in/describe/article/reso/024/06/0653-0659

In this article, we develop the traditional differential equation for Foucault’s pendulum from physical situation and solve it from
standard form. The sublimation of boundary condition eliminates the constants and choice of the local parameters (latitude, pendulum specifications) offers an equation that can be used for a plot followed by animation using MAPLE. The fundamental conceptual components involved in preparing differential equation viz; (i) rotating coordinate system, (ii) rotation of the plane of oscillation and its dependence on the latitude, (iii) effective gravity with latitude, etc., are discussed in detail. The accurate calculations offer quantities up to the sixth decimal point which are used for plotting and animation. This study offers a hands-on experience. Present article offers a know-how to devise a Foucault’s pendulum just by plugging in the latitude of reader’s choice. Students can develop a miniature working model/project of the pendulum.

Exercises solved online with Maple exclusively in space. I attach the explanation links on my YouTube channel.

Part # 01

https://www.youtube.com/watch?v=8Aa2xzU8LwQ

Part # 02

https://www.youtube.com/watch?v=qyGT28CeSz4

Part # 03

https://www.youtube.com/watch?v=yf8rjSPbv5g

Part # 04

https://www.youtube.com/watch?v=FwHPW7ncZTg

Part # 05

https://www.youtube.com/watch?v=bm3frpukb0I

Link for download the file:

Vector_Exercises-Force_in_space.mw

Lenin AC

Ambassador of Maple

 

 

 


 

Solving a Numbrix Puzzle with Logic

Background

 

 

Parade magazine, a filler in the local Sunday newspaper, contains a Numbrix puzzle, the object of which is to find a serpentine path of consecutive integers, 1 through 81, through a nine by nine grid.  The puzzle typically specifies the content of every other border cell.

 

The Maple Logic  package has a procedure, Satisfy , that can be used to solve this puzzle.  Satisfy is a SAT-solver; given a boolean expression it attempts to find a set of equations of the form {x__1 = b__1, x__2 = b__2, () .. ()}, where x__i are the boolean variables in the given expression and b__i are boolean values (true or false) that satisfy the expression (cause it to evaluate true).

 

A general technique to solve this and other puzzles with Satisfy is to select boolean-values variables that encode the state of the puzzle (a trial solution, whether valid or not), generate a boolean-expression of these variables that is satisfied (true) if and only if the variables are given values that correspond to a solution, pass this expression to Satisfy, then translate the returned set of boolean values (if any) to the puzzle solution.

Setup

 

Assign a matrix that defines the grid and the initial position.  Use zeros to indicate the cells that need values. To make it easy to inspect the expressions, a small 2 x 3 matrix is used for this demo---a full size example is given at the end.

M := Matrix(2,3, {(1,1) = 1, (1,3) = 5});

Matrix(2, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 5, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0})

(2.1)

 

Extract the dimensions of the Matrix

(m,n) := upperbound(M);

2, 3

(2.2)

Boolean Variables

 

Let the boolean variable x[i,j,k] mean that cell (i,j) has value k. For example, x[2,3,6] is true when cell (2,3) contains 6, otherwise it is false. There are (m*n)^2 boolean variables.

Initial Position

 

The initial position can be expressed as the following and-clause.

initial := &and(seq(seq(ifelse(M[i,j] = 0, NULL, x[i,j,M[i,j]]), i = 1..m), j = 1..n));

`&and`(x[1, 1, 1], x[1, 3, 5])

(4.1)

Adjacent Cells

 

The requirement that an interior cell with value k is adjacent to the cell with value k+1 can be expressed as the implication
   

   x[i,j,k] &implies &or(x[i-1,j,k+1], x[i+1,j,k+1], x[i,j-1,k+1], x[i,j+1,k+1])

 

Extending this to handle all cells results in the following boolean expression.

adjacent := &and(seq(seq(seq(
         x[i,j,k] &implies &or(NULL
                               , ifelse(1<i, x[i-1, j, k+1], NULL)
                               , ifelse(i<m, x[i+1, j, k+1], NULL)
                               , ifelse(1<j, x[i, j-1, k+1], NULL)
                               , ifelse(j<n, x[i, j+1, k+1], NULL)
                               )
                            , i = 1..m)
                        , j = 1..n)
                    , k = 1 .. m*n-1));

`&and`(`&implies`(x[1, 1, 1], `&or`(x[2, 1, 2], x[1, 2, 2])), `&implies`(x[2, 1, 1], `&or`(x[1, 1, 2], x[2, 2, 2])), `&implies`(x[1, 2, 1], `&or`(x[2, 2, 2], x[1, 1, 2], x[1, 3, 2])), `&implies`(x[2, 2, 1], `&or`(x[1, 2, 2], x[2, 1, 2], x[2, 3, 2])), `&implies`(x[1, 3, 1], `&or`(x[2, 3, 2], x[1, 2, 2])), `&implies`(x[2, 3, 1], `&or`(x[1, 3, 2], x[2, 2, 2])), `&implies`(x[1, 1, 2], `&or`(x[2, 1, 3], x[1, 2, 3])), `&implies`(x[2, 1, 2], `&or`(x[1, 1, 3], x[2, 2, 3])), `&implies`(x[1, 2, 2], `&or`(x[2, 2, 3], x[1, 1, 3], x[1, 3, 3])), `&implies`(x[2, 2, 2], `&or`(x[1, 2, 3], x[2, 1, 3], x[2, 3, 3])), `&implies`(x[1, 3, 2], `&or`(x[2, 3, 3], x[1, 2, 3])), `&implies`(x[2, 3, 2], `&or`(x[1, 3, 3], x[2, 2, 3])), `&implies`(x[1, 1, 3], `&or`(x[2, 1, 4], x[1, 2, 4])), `&implies`(x[2, 1, 3], `&or`(x[1, 1, 4], x[2, 2, 4])), `&implies`(x[1, 2, 3], `&or`(x[2, 2, 4], x[1, 1, 4], x[1, 3, 4])), `&implies`(x[2, 2, 3], `&or`(x[1, 2, 4], x[2, 1, 4], x[2, 3, 4])), `&implies`(x[1, 3, 3], `&or`(x[2, 3, 4], x[1, 2, 4])), `&implies`(x[2, 3, 3], `&or`(x[1, 3, 4], x[2, 2, 4])), `&implies`(x[1, 1, 4], `&or`(x[2, 1, 5], x[1, 2, 5])), `&implies`(x[2, 1, 4], `&or`(x[1, 1, 5], x[2, 2, 5])), `&implies`(x[1, 2, 4], `&or`(x[2, 2, 5], x[1, 1, 5], x[1, 3, 5])), `&implies`(x[2, 2, 4], `&or`(x[1, 2, 5], x[2, 1, 5], x[2, 3, 5])), `&implies`(x[1, 3, 4], `&or`(x[2, 3, 5], x[1, 2, 5])), `&implies`(x[2, 3, 4], `&or`(x[1, 3, 5], x[2, 2, 5])), `&implies`(x[1, 1, 5], `&or`(x[2, 1, 6], x[1, 2, 6])), `&implies`(x[2, 1, 5], `&or`(x[1, 1, 6], x[2, 2, 6])), `&implies`(x[1, 2, 5], `&or`(x[2, 2, 6], x[1, 1, 6], x[1, 3, 6])), `&implies`(x[2, 2, 5], `&or`(x[1, 2, 6], x[2, 1, 6], x[2, 3, 6])), `&implies`(x[1, 3, 5], `&or`(x[2, 3, 6], x[1, 2, 6])), `&implies`(x[2, 3, 5], `&or`(x[1, 3, 6], x[2, 2, 6])))

(5.1)

 

All Values Used

 

The following expression is true when each integer k, from 1 to m*n, is assigned to one or more cells.

allvals := &and(seq(seq(&or(seq(x[i,j,k], k=1..m*n)), i=1..m), j=1..n));

`&and`(`&or`(x[1, 1, 1], x[1, 1, 2], x[1, 1, 3], x[1, 1, 4], x[1, 1, 5], x[1, 1, 6]), `&or`(x[2, 1, 1], x[2, 1, 2], x[2, 1, 3], x[2, 1, 4], x[2, 1, 5], x[2, 1, 6]), `&or`(x[1, 2, 1], x[1, 2, 2], x[1, 2, 3], x[1, 2, 4], x[1, 2, 5], x[1, 2, 6]), `&or`(x[2, 2, 1], x[2, 2, 2], x[2, 2, 3], x[2, 2, 4], x[2, 2, 5], x[2, 2, 6]), `&or`(x[1, 3, 1], x[1, 3, 2], x[1, 3, 3], x[1, 3, 4], x[1, 3, 5], x[1, 3, 6]), `&or`(x[2, 3, 1], x[2, 3, 2], x[2, 3, 3], x[2, 3, 4], x[2, 3, 5], x[2, 3, 6]))

(6.1)

Single Value

 

The following expression is satisfied when each cell has no more than one value.

 single := &not &or(seq(seq(seq(seq(x[i,j,k] &and x[i,j,kk], kk = k+1..m*n), k = 1..m*n-1), i = 1..m), j = 1..n));

`&not`(`&or`(`&and`(x[1, 1, 1], x[1, 1, 2]), `&and`(x[1, 1, 1], x[1, 1, 3]), `&and`(x[1, 1, 1], x[1, 1, 4]), `&and`(x[1, 1, 1], x[1, 1, 5]), `&and`(x[1, 1, 1], x[1, 1, 6]), `&and`(x[1, 1, 2], x[1, 1, 3]), `&and`(x[1, 1, 2], x[1, 1, 4]), `&and`(x[1, 1, 2], x[1, 1, 5]), `&and`(x[1, 1, 2], x[1, 1, 6]), `&and`(x[1, 1, 3], x[1, 1, 4]), `&and`(x[1, 1, 3], x[1, 1, 5]), `&and`(x[1, 1, 3], x[1, 1, 6]), `&and`(x[1, 1, 4], x[1, 1, 5]), `&and`(x[1, 1, 4], x[1, 1, 6]), `&and`(x[1, 1, 5], x[1, 1, 6]), `&and`(x[2, 1, 1], x[2, 1, 2]), `&and`(x[2, 1, 1], x[2, 1, 3]), `&and`(x[2, 1, 1], x[2, 1, 4]), `&and`(x[2, 1, 1], x[2, 1, 5]), `&and`(x[2, 1, 1], x[2, 1, 6]), `&and`(x[2, 1, 2], x[2, 1, 3]), `&and`(x[2, 1, 2], x[2, 1, 4]), `&and`(x[2, 1, 2], x[2, 1, 5]), `&and`(x[2, 1, 2], x[2, 1, 6]), `&and`(x[2, 1, 3], x[2, 1, 4]), `&and`(x[2, 1, 3], x[2, 1, 5]), `&and`(x[2, 1, 3], x[2, 1, 6]), `&and`(x[2, 1, 4], x[2, 1, 5]), `&and`(x[2, 1, 4], x[2, 1, 6]), `&and`(x[2, 1, 5], x[2, 1, 6]), `&and`(x[1, 2, 1], x[1, 2, 2]), `&and`(x[1, 2, 1], x[1, 2, 3]), `&and`(x[1, 2, 1], x[1, 2, 4]), `&and`(x[1, 2, 1], x[1, 2, 5]), `&and`(x[1, 2, 1], x[1, 2, 6]), `&and`(x[1, 2, 2], x[1, 2, 3]), `&and`(x[1, 2, 2], x[1, 2, 4]), `&and`(x[1, 2, 2], x[1, 2, 5]), `&and`(x[1, 2, 2], x[1, 2, 6]), `&and`(x[1, 2, 3], x[1, 2, 4]), `&and`(x[1, 2, 3], x[1, 2, 5]), `&and`(x[1, 2, 3], x[1, 2, 6]), `&and`(x[1, 2, 4], x[1, 2, 5]), `&and`(x[1, 2, 4], x[1, 2, 6]), `&and`(x[1, 2, 5], x[1, 2, 6]), `&and`(x[2, 2, 1], x[2, 2, 2]), `&and`(x[2, 2, 1], x[2, 2, 3]), `&and`(x[2, 2, 1], x[2, 2, 4]), `&and`(x[2, 2, 1], x[2, 2, 5]), `&and`(x[2, 2, 1], x[2, 2, 6]), `&and`(x[2, 2, 2], x[2, 2, 3]), `&and`(x[2, 2, 2], x[2, 2, 4]), `&and`(x[2, 2, 2], x[2, 2, 5]), `&and`(x[2, 2, 2], x[2, 2, 6]), `&and`(x[2, 2, 3], x[2, 2, 4]), `&and`(x[2, 2, 3], x[2, 2, 5]), `&and`(x[2, 2, 3], x[2, 2, 6]), `&and`(x[2, 2, 4], x[2, 2, 5]), `&and`(x[2, 2, 4], x[2, 2, 6]), `&and`(x[2, 2, 5], x[2, 2, 6]), `&and`(x[1, 3, 1], x[1, 3, 2]), `&and`(x[1, 3, 1], x[1, 3, 3]), `&and`(x[1, 3, 1], x[1, 3, 4]), `&and`(x[1, 3, 1], x[1, 3, 5]), `&and`(x[1, 3, 1], x[1, 3, 6]), `&and`(x[1, 3, 2], x[1, 3, 3]), `&and`(x[1, 3, 2], x[1, 3, 4]), `&and`(x[1, 3, 2], x[1, 3, 5]), `&and`(x[1, 3, 2], x[1, 3, 6]), `&and`(x[1, 3, 3], x[1, 3, 4]), `&and`(x[1, 3, 3], x[1, 3, 5]), `&and`(x[1, 3, 3], x[1, 3, 6]), `&and`(x[1, 3, 4], x[1, 3, 5]), `&and`(x[1, 3, 4], x[1, 3, 6]), `&and`(x[1, 3, 5], x[1, 3, 6]), `&and`(x[2, 3, 1], x[2, 3, 2]), `&and`(x[2, 3, 1], x[2, 3, 3]), `&and`(x[2, 3, 1], x[2, 3, 4]), `&and`(x[2, 3, 1], x[2, 3, 5]), `&and`(x[2, 3, 1], x[2, 3, 6]), `&and`(x[2, 3, 2], x[2, 3, 3]), `&and`(x[2, 3, 2], x[2, 3, 4]), `&and`(x[2, 3, 2], x[2, 3, 5]), `&and`(x[2, 3, 2], x[2, 3, 6]), `&and`(x[2, 3, 3], x[2, 3, 4]), `&and`(x[2, 3, 3], x[2, 3, 5]), `&and`(x[2, 3, 3], x[2, 3, 6]), `&and`(x[2, 3, 4], x[2, 3, 5]), `&and`(x[2, 3, 4], x[2, 3, 6]), `&and`(x[2, 3, 5], x[2, 3, 6])))

(7.1)

Solution

 

Combine the boolean expressions into a a single and-clause and pass it to Satisfy.

sol := Logic:-Satisfy(&and(initial, adjacent, allvals, single));

{x[1, 1, 1] = true, x[1, 1, 2] = false, x[1, 1, 3] = false, x[1, 1, 4] = false, x[1, 1, 5] = false, x[1, 1, 6] = false, x[1, 2, 1] = false, x[1, 2, 2] = false, x[1, 2, 3] = false, x[1, 2, 4] = false, x[1, 2, 5] = false, x[1, 2, 6] = true, x[1, 3, 1] = false, x[1, 3, 2] = false, x[1, 3, 3] = false, x[1, 3, 4] = false, x[1, 3, 5] = true, x[1, 3, 6] = false, x[2, 1, 1] = false, x[2, 1, 2] = true, x[2, 1, 3] = false, x[2, 1, 4] = false, x[2, 1, 5] = false, x[2, 1, 6] = false, x[2, 2, 1] = false, x[2, 2, 2] = false, x[2, 2, 3] = true, x[2, 2, 4] = false, x[2, 2, 5] = false, x[2, 2, 6] = false, x[2, 3, 1] = false, x[2, 3, 2] = false, x[2, 3, 3] = false, x[2, 3, 4] = true, x[2, 3, 5] = false, x[2, 3, 6] = false}

(8.1)

Select the equations whose right size is true

sol := select(rhs, sol);

{x[1, 1, 1] = true, x[1, 2, 6] = true, x[1, 3, 5] = true, x[2, 1, 2] = true, x[2, 2, 3] = true, x[2, 3, 4] = true}

(8.2)

Extract the lhs of the true equations

vars := map(lhs, sol);

{x[1, 1, 1], x[1, 2, 6], x[1, 3, 5], x[2, 1, 2], x[2, 2, 3], x[2, 3, 4]}

(8.3)

Extract the result from the indices of the vars and assign to a new Matrix

S := Matrix(m,n):

for v in vars do S[op(1..2,v)] := op(3,v); end do:

S;

Matrix(2, 3, {(1, 1) = 1, (1, 2) = 6, (1, 3) = 5, (2, 1) = 2, (2, 2) = 3, (2, 3) = 4})

(8.4)

Procedure

 

We can now combine the manual steps into a procedure that takes an initialized Matrix and fills in a solution.

Numbrix := proc( M :: ~Matrix, { inline :: truefalse := false } )

Example

 

Create the initial position for a 9 x 9 Numbrix and solve it.

P := Matrix(9, {(1,1)=11, (1,3)=7, (1,5)=3, (1,7)=81, (1,9)=77, (3,9)=75, (5,9)=65, (7,9)=55, (9,9)=53
               , (9,7)=47, (9,5)=41, (9,3)=39, (9,1)=37, (7,1)=21, (5,1)=17, (3,1)=13});

Matrix(9, 9, {(1, 1) = 11, (1, 2) = 0, (1, 3) = 7, (1, 4) = 0, (1, 5) = 3, (1, 6) = 0, (1, 7) = 81, (1, 8) = 0, (1, 9) = 77, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 13, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 75, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (5, 1) = 17, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (5, 9) = 65, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (6, 9) = 0, (7, 1) = 21, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (7, 9) = 55, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0, (8, 9) = 0, (9, 1) = 37, (9, 2) = 0, (9, 3) = 39, (9, 4) = 0, (9, 5) = 41, (9, 6) = 0, (9, 7) = 47, (9, 8) = 0, (9, 9) = 53})

(10.1)

CodeTools:-Usage(Numbrix(P));

memory used=0.77GiB, alloc change=220.03MiB, cpu time=15.55s, real time=12.78s, gc time=3.85s

 

Matrix(9, 9, {(1, 1) = 11, (1, 2) = 10, (1, 3) = 7, (1, 4) = 81, (1, 5) = 3, (1, 6) = 4, (1, 7) = 81, (1, 8) = 78, (1, 9) = 77, (2, 1) = 12, (2, 2) = 9, (2, 3) = 8, (2, 4) = 7, (2, 5) = 6, (2, 6) = 5, (2, 7) = 80, (2, 8) = 79, (2, 9) = 76, (3, 1) = 13, (3, 2) = 14, (3, 3) = 27, (3, 4) = 28, (3, 5) = 71, (3, 6) = 72, (3, 7) = 73, (3, 8) = 74, (3, 9) = 75, (4, 1) = 16, (4, 2) = 15, (4, 3) = 26, (4, 4) = 29, (4, 5) = 70, (4, 6) = 69, (4, 7) = 68, (4, 8) = 67, (4, 9) = 66, (5, 1) = 17, (5, 2) = 18, (5, 3) = 25, (5, 4) = 30, (5, 5) = 61, (5, 6) = 62, (5, 7) = 63, (5, 8) = 64, (5, 9) = 65, (6, 1) = 20, (6, 2) = 19, (6, 3) = 24, (6, 4) = 31, (6, 5) = 60, (6, 6) = 59, (6, 7) = 58, (6, 8) = 57, (6, 9) = 56, (7, 1) = 21, (7, 2) = 22, (7, 3) = 23, (7, 4) = 32, (7, 5) = 43, (7, 6) = 44, (7, 7) = 49, (7, 8) = 50, (7, 9) = 55, (8, 1) = 36, (8, 2) = 35, (8, 3) = 34, (8, 4) = 33, (8, 5) = 42, (8, 6) = 45, (8, 7) = 48, (8, 8) = 51, (8, 9) = 54, (9, 1) = 37, (9, 2) = 38, (9, 3) = 39, (9, 4) = 40, (9, 5) = 41, (9, 6) = 46, (9, 7) = 47, (9, 8) = 52, (9, 9) = 53})

(10.2)

 

numbrix.mw

1 2 3 4 5 6 7 Last Page 1 of 54