vv

12453 Reputation

19 Badges

9 years, 284 days

MaplePrimes Activity


These are answers submitted by vv

The modp1 arithmetic is very fast in maple (but especially when the modulus is a hardware integer).
The usual mod operations use modp1, so they are fast too.
Using modp1 directly, the gain is modest:

restart;
m := 8009;
N := ceil(m/2) + 2;
read "nn.txt":
st := time[real]():
r := modpol(a*a, f_t, t, 2^N):
t1 := time[real]() - st;

N2:=2^N:
A:=modp1(ConvertIn(a,t),N2):
F:=modp1(ConvertIn(f_t,t),N2):

st := time[real]():
R:=modp1(Rem(Multiply(A,A),F), N2):
t2 := time[real]() - st;

rr:=modp1(ConvertOut(R,t),N2):
r-rr; # 0 
t2/t1;

st := time[real]():
modp1(Rem(Power(A,2),F), N2):
t3 := time[real]() - st;
t3/t1;

So, probably in your case that's all Maple can do.

For a general approach, you can use the package QuadraticInt (for working in the ring Z[sqrt(d)]); see the Application Center.

restart;

S:=Sum((-1)^(n+1)/(n*(n+1)),n = 1 .. infinity):
S=value(S);

Sum((-1)^(n+1)/(n*(n+1)), n = 1 .. infinity) = -1+2*ln(2)

(1)

# For a proof, let's change the variable in Sum

'S' = subs((n = 1.. infinity) = (k=0..infinity), n=k+1, S);

S = Sum((-1)^(k+2)/((k+1)*(k+2)), k = 0 .. infinity)

(2)

F:=int(ln(1+x),x);

ln(1+x)*(1+x)-1-x

(3)

Sx:=convert(F, FPS);

-1+Sum((-1)^(k+2)*x^(k+2)/((k+2)*(k+1)), k = 0 .. infinity)

(4)

# So,

S = eval(F + 1, x = 1);

Sum((-1)^(n+1)/(n*(n+1)), n = 1 .. infinity) = -1+2*ln(2)

(5)


 

 

Your Bst is the parametric representation of a surface.
But it is in the plane z = 0. So you obviously have a single "contour", but this "contour" is actually a 2-dimensional set in R^2 (it has interior points).

isolve is not able to manage inequalities in this case. Select the desired solutions later:

isolve(29 = x^2 + y^2):
select(u -> (eval(x,u)>=1 and eval(y,u)>=1), [%]);

 

It does not work because eval needs subexpressions (i.e. operands of some level).
[E.g.  a*b  is not subexpression in a*b*c;  for such substitution there is algsubs, but in Physics almost everything is redefined, so this probably fails too.]

But here you can simply use:

eval((2), Dagger(A[i]) = -A[i]);

 

 

If you just want symbolic computations, you may use the Physics package:
with(Physics):
Setup(noncommutativeprefix = {A,B});

Now, the adjoint (Hermitian conjugate) ^* will act (almost) as a generic involution (displayed: Dagger).

(A^(-1) * B)^*; # ==>  B^*  * (A^*)^(-1)

n:=10;  u:=5.2;
V:=LinearAlgebra:-RandomVector(n, generator=rand(0.0 .. 10.0));

n := 10

 

u := 5.2

 

Vector[column](%id = 18446744074617251110)

(1)

s:=add(b[i]*V[i],i=1..n);

2.075595849*b[1]+8.602833750*b[2]+6.237874345*b[3]+3.783826445*b[4]+6.378864096*b[5]+6.726260294*b[6]+.5091763968*b[7]+3.694547905*b[8]+7.678663330*b[9]+8.996029595*b[10]

(2)

sol:=Optimization:-LPSolve(s-u, [u<=s], assume=binary, depthlimit=10000);

[.570143754, [b[1] = 1, b[2] = 0, b[3] = 0, b[4] = 0, b[5] = 0, b[6] = 0, b[7] = 0, b[8] = 1, b[9] = 0, b[10] = 0]]

(3)

remove(`=`, eval([op(s)],sol[2]),0);

[2.075595849, 3.694547905]

(4)


Download mimiz.mw

A faster version.

F := proc(N::posint)
  local V:=Vector(4), p:=2;
  while (p:=nextprime(p))<N do V[(irem(p,8)+1)/2]++ od;
  [seq(V)]
end proc:

 

A := Matrix(2, 3, [[1, 2, 3], [3, 1, 2]]); B := Matrix(2, 3, [[1, x, 3], [y, 1, z]]);
solve(Equate(A,B), {x,y,z});

Note that solve needs a list or set for its first argument, so solve(A=B, {x,y,z}); will not work.
The second argument is optional, so, solve(Equate(A,B)); is enough.

simplify ... assuming ... should be enough.
Unfortunately an extra step is needed, for example:

arctan(cos(beta)/sin(beta));
simplify(convert(%, tan)) assuming beta>0, beta<Pi/2;

   Pi/2 - beta

 

When a graph is obtained via a "constructor" such as CycleGraph or CompleteGraph,
some known properties (such as "chromatic_number" or "clique_number")  will be included as attributes.
But the commands DeleteEdge, AddEdge  etc  simply keep these attributes in the new graph, which is wrong.
For example, CliqueNumber and CromaticNumber will use these attributes instead of (re)computing them.

restart;
with(GraphTheory):
G := CycleGraph(3); # or CompleteGraph(3)
H := DeleteEdge(G, {1,2}, inplace=false):
H1:=Graph(Vertices(H), Edges(H)):  # identical to H but without attributes

GRAPHLN(undirected, unweighted, [1, 2, 3], Array(%id = 18446744074327621750), `GRAPHLN/table/1`, 0)

(1)

GetGraphAttribute(G);

["chromatic_number" = 3, "clique_cover_number" = 1, "clique_number" = 3, "independence_number" = 1]

(2)

GetGraphAttribute(H); # copied from G (incorrect!)

["chromatic_number" = 3, "clique_cover_number" = 1, "clique_number" = 3, "independence_number" = 1]

(3)

GetGraphAttribute(H1);

[]

(4)

CliqueNumber(H), CliqueNumber(H1);

3, 2

(5)

ChromaticNumber(H), ChromaticNumber(H1);

3, 2

(6)

Yes, it's a bug. With the col parameter the result is correct.

restart;
with(GraphTheory):
G1:=AddEdge(CycleGraph(4),{{1,3}}):
ChromaticNumber(G1);       # 2   wrong
ChromaticNumber(G1,'col'); # 3   correct
col;                       # [[2, 4], [3], [1]]

Edit. Strangely, all the methods give a wrong result:

for met in [hybrid, optimal, brelaz, dsatur, greedy, welshpowell, sat] do
  ChromaticNumber(G1, 'method'=met);
od;

Actually it seems that the bug is in GraphTheory:-GraphInfo:-GetAttrib
called by GraphTheory:-GetGraphAttribute.

For details and examples see  ?parallelism

P.S. It has happened to me many times that a Google search including a keyword "maple" is more successful than a search in Maple, even if the information is present in the help system.

 

I wrote a couple of years ago a procedure for this (multiple integrals over domains defined by polynomial inequalities).
https://www.maplesoft.com/applications/view.aspx?SID=154350

Actually, in that code the order of integration was chosen automatically.
Here is a version where the order is set by the user.

MultiInt:=proc(f::algebraic, rel::set(relation), Lvars::list(name))
local S,  Dom,u, vars:={Lvars[]};
S:=solve(rel,Lvars);
if _SolutionsMayBeLost=true then print("solve says 'SolutionsMayBeLost'") fi;
if type(S,list) then S:=remove(hastype,S,`=`) else error "solve cannot decompose" fi;
Dom:=proc(r)
  local z,a,b, found:={},dom:=NULL,ok ;
  while found<>vars do   
     for z in vars minus found do
     ok:=false;
     a:=lhs~(select(t ->rhs(t)=z and type(lhs(t), freeof(vars minus found)), r));
     b:=rhs~(select(t ->lhs(t)=z and type(rhs(t), freeof(vars minus found)), r));   
     if nops(a)=1 and nops(b)=1 then
        found:=found union {z}; dom:=z=a[]..b[],dom; ok:=true; break;
     fi;
     od;
     if not ok then return FAIL fi
  od;  
  Int(f,dom)
end proc:  # Dom
add(Dom(u), u=S)
end proc:

 

In your case the domain is defined by the inequalities:

dom := {(4*x-x^2)<=(y^2), y^2<=4*x, x<=4, x>=0, y>=0};

{0 <= x, 0 <= y, x <= 4, y^2 <= 4*x, -x^2+4*x <= y^2}

(1)

MultiInt( f(x,y), dom, [x,y]);

Int(f(x, y), y = (-x^2+4*x)^(1/2) .. 2*x^(1/2), x = 0 .. 4)

(2)

MultiInt( f(x,y), dom, [y,x]);

Int(f(x, y), x = (1/4)*y^2 .. 2-(-y^2+4)^(1/2), y = 0 .. 2)+Int(f(x, y), x = 2+(-y^2+4)^(1/2) .. 4, y = 0 .. 2)+Int(f(x, y), x = (1/4)*y^2 .. 4, y = 2 .. 4)

(3)

# Take now f = 1

MultiInt( 1,dom, [x,y]);

Int(1, y = (-x^2+4*x)^(1/2) .. 2*x^(1/2), x = 0 .. 4)

(4)

value(%);

32/3-2*Pi

(5)

MultiInt( 1, dom, [y,x]);

Int(1, x = (1/4)*y^2 .. 2-(-y^2+4)^(1/2), y = 0 .. 2)+Int(1, x = 2+(-y^2+4)^(1/2) .. 4, y = 0 .. 2)+Int(1, x = (1/4)*y^2 .. 4, y = 2 .. 4)

(6)

value(%);

32/3-2*Pi

(7)

 

First 31 32 33 34 35 36 37 Last Page 33 of 111