Kitonum

16053 Reputation

24 Badges

12 years, 199 days

MaplePrimes Activity


These are Posts that have been published by Kitonum

This post is related to the this thread

The recursive procedure  PosIntSolve  finds the number of non-negative or positive solutions of any linear Diophantine equation

a1*x1+a2*x2+ ... +aN*xN = n  with positive coefficients a1, a2, ... aN .  

Formal parameters: L is the list of coefficients of the left part, n  is the right part,  s (optional) is nonneg (by default) for nonnegint solutions and  pos  for positive solutions.

The basic ideas:

1) If you make shifts of the all unknowns by the formulas  x1'=x1-1,  x2'=x2-1, ... , xN'=xN-1  then  the number of positive solutions of the first equation equals the number of non-negative solutions of the second equation.

2) The recurrence formula (penultimate line of the procedure) can easily be proved by induction.

 

The code of the procedure:

restart;

PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)

local N, n0;

option remember;

if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;

N:=nops(L);

if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;

add(PosIntSolve(subsop(1=NULL,L),n0-k*L[1]), k=0..floor(n0/L[1]));

end proc:

 

Examples of use.

 

Finding of the all positive solutions of equation 30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

st:=time():

PosIntSolve([30,75,110,85,255,160,15,12,120], 8000, pos);

time()-st;

                                       13971409380

                                             2.125

 

To test the procedure, solve (separately for non-negative and positive solutions) the simple equation  2*x1+7*x2+3*x3=2000  in two ways (by the  procedure and brute force method):

ts:=time():

PosIntSolve([2,7,3], 2000);

PosIntSolve([2,7,3], 2000, pos);

time()-ts;

                47905

                47334

                 0.281

 

ts:=time():

k:=0:

for x from 0 to 2000/2 do

for y from 0 to floor((2000-2*x)/7) do

for z from 0 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k; 

k:=0:

for x from 1 to 2000/2 do

for y from 1 to floor((2000-2*x)/7) do

for z from 1 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k;

time()-ts; 

                   47905

                   47334

                   50.063

 

Another example - the solution of the famous problem: how many ways can be exchanged $ 1 using the coins of smaller denomination.

PosIntSolve([1,5,10,25,50],100);

                        292

 

 Number-of-solutions.mw

 

 Edit.  The code has been slightly edited 

This post is related to the question. There were  proposed two ways of finding the volume of the cutted part of a sphere in the form of a wedge.  Here the procedure is presented that shows the rest of the sphere. Parameters procedure: R - radius of the sphere, H1 - the distance the first cutting plane to the plane  xOy,  H2 -  the distance the second cutting plane to the plane  zOy. Necessary conditions:  R>0,  H1>=0,  H2>=0,  H1^2+H2^2<R^2 . For clarity, different surfaces are painted in different colors.

restart;

Pic := proc (R::positive, H1::nonnegative, H2::nonnegative)

local A, B, C, E, F;

if R^2 <= H1^2+H2^2 then error "Should be H1^(2)+H2^(2)<R^(2)" end if;

A := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = arctan(sqrt(-H1^2-H2^2+R^2), H2) .. 2*Pi-arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. Pi, color = green);

B := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. arccos(sqrt(R^2-H2^2-H2^2*tan(phi)^2)/R), color = green);

C := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = arccos(H1/R) .. Pi, color = green);

E := plot3d([r*cos(phi), r*sin(phi), H1], phi = -arccos(H2/sqrt(R^2-H1^2)) .. arccos(H2/sqrt(R^2-H1^2)), r = H2/cos(phi) .. sqrt(R^2-H1^2), color = blue);

F := plot3d([H2, r*cos(phi), r*sin(phi)], phi = arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)) .. Pi-arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)), r = H1/sin(phi) .. sqrt(R^2-H2^2), color = gold);

plots[display](A, B, C, E, F, axes = none, view = [-1.5 .. 1.5, -1.5 .. 1.5, -1.5 .. 1.5], scaling = constrained, lightmodel = light4, orientation = [60, 80]);

end proc:

 

Example of use:

Pic(1,  0.5,  0.3);

                             

 

 

The work contains two procedures called  SetPartition  and  NumbPart .

The first procedure  SetPartition  generates all the partitions of a set  S  into disjoint subsets of specific sizes defined by a list  Q. The third optional parameter is the symbol  `ordered`  or  `nonordered`  (by default) . It determines whether the order of subsets in a partition is significant.

The second procedure  NumbPart  returns the number of all partitions of the sizes indicated by  Q .

 

Codes of these procedures:

SetPartition:=proc(S::set, Q::list(posint), K::symbol:=nonordered)  # Procedure finds all partitions of the set  S  into subsets of the specific size given Q.

local L, P, T, S1, S2, n, i, j, m, k, M;

uses ListTools,combinat;

if `+`(op(Q))<>nops(S) then error "Should be `+`(op(Q))=nops(S)" fi;

L:=convert(Q,set);

T:=[seq([L[i],Occurrences(L[i],Q)], i=1..nops(L))];

P:=`if`(K=ordered,convert(T,list),convert(T,set));

S1:=S;  S2:={`if`(K=ordered,[],{})}; n:=nops(P);

for i to n do

m:=P[i,1]; k:=P[i,2];

for j to k do

S2:={seq(seq(`if`(K=ordered,[op(s),t],{op(s),t}), t=choose(S1 minus `union`(op(s)),m)), s=S2)};

od; od;

if K=ordered then {map(op@permute,S2)[]} else S2 fi;

end proc:

 

NumbPart:=proc(Q::list(posint), K::symbol:=nonordered)  # Procedure finds the number of all partitions of a set into subsets of the specific size given  Q

local L, T, P, n, S, N;

uses ListTools;

L:=convert(Q,set);

T:=[seq([L[i],Occurrences(L[i],Q)], i=1..nops(L))];

P:=`if`(K=ordered,convert(T,list),convert(T,set));

n:=nops(P);  N:=add(P[i,2], i=1..n);

S:=add(P[i,1]*P[i,2],i=1..n)!/mul(P[i,1]!^P[i,2],i=1..n)/mul(P[i,2]!,i=1..n);

if K=nonordered then return S else  S*N! fi;

end proc:

 

Examples of use:

SetPartition({a,b,c,d,e}, [1,1,3]);  nops(%);  # Nonordered partitions and their number

SetPartition({a,b,c,d,e}, [1,1,3], ordered);  nops(%);  # Ordered partitions and their number

 

 

Here's a more interesting example. 5 fruits  {apple, pear, orange, mango, peach}  must be put on three plates so that  on each of two plates there are 2  fruits, and there is one fruit  on one plate. Two variants to solve: 1) plates are indistinguishable and 2) different plates. In how many ways can this be done?

SetPartition({apple,pear,orange,mango,peach}, [1,2,2]);  nops(%);  # plates are indistinguishable

 

NumbPart([1,2,2], ordered);  # Number of ways for  different plates

                                                              90

 

Another example - how many ways can be divided  100 objects into 10 equal parts?

NumbPart([10$10]);

    64954656894649578274066349293466217242333450230560675312538868633528911487364888307200

 

 SetPartition.mws

The procedure  Partition  significantly generalizes the standard procedure  combinat[partition]  in several ways. The user specifies the number of parts of the partition, and can also set different limitations on parts partition.

Required parameters:  n - a nonnegative integer, - a positive integer or a range (k  specifies the number of parts of the partition). The parameter  res  is the optional parameter (by default  res is  ). If  res  is a number, all elements of  k-tuples must be greater than or equal  res .  If  res  is a range  a .. b ,   all elements of  k-tuples must be greater than or equal  a  and  less than or equal  b . The optional parameter  S  - set, which includes elements of the partition. By default  S = {$ 0.. n} .

The code of the procedure:

Partition:=proc(n::nonnegint, k::{posint,range}, res::{range, nonnegint} := 1, S::set:={$0..n})  # Generates a list of all partitions of an integer n into k parts

local k_Partition, n1, k1, L;

 

k_Partition := proc (n, k::posint, res, S)

local m, M, a, b, S1, It, L0;

m:=S[1]; M:=S[-1];

if res::nonnegint then a := max(res,m); b := min(n-(k-1)*a,M)  else a := max(lhs(res),m); b := min(rhs(res),M) fi;

S1:={$a..b} intersect S;

if b < a or b*k < n or a*k > n  then return [ ] fi;

It := proc (L)

local m, j, P, R, i, N;

m := nops(L[1]); j := k-m; N := 0;

for i to nops(L) do

R := n-`+`(op(L[i]));

if R <= b*j and a*j <= R then N := N+1;

P[N] := [seq([op(L[i]), s], s = {$ max(a, R-b*(j-1)) .. min(R, b)} intersect select(t->t>=L[i,-1],S1) )] fi;

od;

[seq(op(P[s]), s = 1 .. N)];

end proc;

if k=1 then [[b]] else (It@@(k-1))(map(t->[t],S1))  fi;

end proc;

 

if k::posint then return k_Partition(n,k,res,S) else n1:=0;

for k1 from lhs(k) to rhs(k) do

n1:=n1+1; L[n1]:=k_Partition(n,k1,res,S)

od;

L:=convert(L,list);

[seq(op(L[i]), i=1..n1)] fi;

 

end proc:

 

Examples of use:

Partition(15, 3);

 

 

Partition(15, 3..5, 1..5);  # The number of parts from 3 to 5, and each summand from 1 to 5

 

 

Partition(15, 5, {seq(2*n-1, n=1..8)});  # 5 summands and all are odd numbers 

 

 

A more interesting example.
There are  k banknotes in possible denominations of 5, 10, 20, 50, 100 dollars. At what number of banknotes  k  the number of variants of exchange  $140  will be maximum?

n:=0:

for k from 3 to 28 do

n:=n+1: V[n]:=[k, nops(Partition(140, k, {5,10,20,50,100}))];

od:

V:=convert(V, list);

max(seq(V[i,2], i=1..nops(V)));

select(t->t[2]=8, V);

 

Here are these variants:

Partition(140, 10, {5,10,20,50,100});

Partition(140, 13, {5,10,20,50,100});

 

 Partition.mws

 

 

 

A heart shape in 3d:

 

 

The code of the animation:

A := plots[animate](plot3d, [[16*sin(t)^3*cos(s), 16*sin(t)^3*sin(s), 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)], t = 0 .. u, s = 0 .. 2*Pi, color = red, style = surface, axes = none], u = 0 .. Pi, frames = 100):

B := plots[animate](plot3d, [[16*sin(t)^3*cos(s), 16*sin(t)^3*sin(s), 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)], t = u .. Pi, s = 0 .. 2*Pi, color = "LightBlue", style = surface, axes = none], u = 0 .. Pi, frames = 100):

plots[display](A, B);

 

Edited. The direction of painting changed.

 

4 5 6 7 8 9 10 Page 6 of 11