## Alec Mihailovs

Dr. Aleksandrs Mihailovs

## 4470 Reputation

19 years, 353 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

## Social Networks and Content at Maplesoft.com

I received my Ph.D. from the University of Pennsylvania in 1998 and I have been teaching since then at SUNY Oneonta for 1 year, at Shepherd University for 5 years, at Tennessee Tech for 2 years, at Lane College for 1 year, and this year I taught at the University of Massachusetts Lowell. My research interests include Representation Theory and Combinatorics.

## Hypercube Line Picking...

See Eric Weisstein's page Hypercube Line Picking in the Mathworld.

Alec

1. That happens with other people, too, including me - see that thread. I usually copy everything before clicking "submit", then seeing an empty comment, click "Edit" under it, and paste there what I just copied. That doesn't give you points on this MaplePrimes version - you get 5 points for every vote up of your message, or if somebody adds it to his/her favorites. No points for just posting. I could delete those your messages if you don't have that option (under the message).

2. If you can get a 32-bit version of Classic worksheet, it might work with your otherwise 64-bit Maple installation - see Dave Linder's post explaining the necessary steps in Linux.

3. If you can't get Classic, there is a possibility of using Command line Maple. It can be run in PowerShell, by the way.

Alec

## closed form...

For symbolical minimization, you don't necessarily need the closed form of this integral. You need the integral of derivative of F(x,y) with respect to y.

Alec

## Classic Maple...

Use Classic Maple (like most of us, MaplePrimes, do) - no problems with copying and pasting.

Alec

## By hand...

My guess is that the best way in such case is to do that by hand.

For example (not directly including the determinat computation, but including a sparse LUDecomposition which could be used for the determinant calculation) the following seems to involve Maple in a very long process,

```B:=Vector(10^6,i->i,datatype=float[8],storage=sparse):
A:=LinearAlgebra:-DiagonalMatrix(B,datatype=float[8],storage=sparse):
x:=LinearAlgebra:-LinearSolve(A,B,method=SparseLU);
```

A simple hand-written procedure would first check if the number of 1s is less than the matrix size - then the determinant is 0, then check that each row and each column has at least one 1 - if not, the determinant is 0, then check if the matrix is lower triangular or upper triangular - then the determinant is 1 if all diagonal entries are 1 and 0 otherwise. Then the case of permutation matrices could be done - if the number of entries equals the matrix size, then the determinant is ±1 depending on the permutation sign. If that doesn't work - then an attempt of doing a sparse LU decomposition could be done. Different methods are known depending on the matrix structure. That could be compiled (using Compiler:-Compile).

For smaller sizes, Linbox might be useful, but I don't think it would be helpful for the sizes you would like to use.

Alec

## convert FromMma...

You could use converter from Mma to Maple,

```convert("f[#, a, #, b] &[x, y]",FromMma);

unapply(f(_Z1, a, _Z1, b), _Z1)(x, y)

%;

f(x, a, x, b)
```

Alec

## help pages...

I put a link to the corresponding help page in the answer. It has several links to other pages. Some examples (but not many) could be found on this site, especially after its search function will be improved.

About books or other texts - it is in section 6.3 and the end of section 6.5 in the Programming Guide, p. 218 about the seq - but it looks as if it just repeats the help pages.

Alec

## For example...

For example,

```Pk:=Array(indets(plot(f(x),x=0..2*Pi),listlist)[],datatype=float[8]);

[ 1..200 x 1..2 2-D Array ]
Pk := [ Data Type: float[8]     ]
[ Storage: rectangular    ]
[ Order: Fortran_order    ]

plots:-pointplot(Pk);
```

and you could achieve a similar effect using style=point in the plot command,

`plot(f(x),x=0..2*Pi,style=point);`

Array dimensions can be found as

```ArrayDims(Pk);

1 .. 200, 1 .. 2

op([2,1,2],Pk);

200
```

The step is not constant,

```[min,max](seq(Pk[i+1,1]-Pk[i,1],i=1..199));

[0.0287234985352272, 0.0344642955948600]
```

Alec

## Simple...

```diff(a(z(tau)),tau);

/ d         \
D(a)(z(tau)) |---- z(tau)|
\dtau       /

convert(%,diff);

/ d       \|            / d         \
|--- a(t1)||            |---- z(tau)|
\dt1      /|t1 = z(tau) \dtau       /

convert(%,D);

D(a)(z(tau)) D(z)(tau)
```

So one can chose the form that (s)he prefers. Derivatives can be displayed in various formats using ?PDEtools,declare (as well as typesetting rules in Standard Maple with extended typesetting.)

Alec

## _passed and seq modifier...

For example,

```f:=()->add(i,i=_passed):

f(1,2,3,4);

10
```

which, in this example, could be written more simple as

`f:=`+`:`

The number of parameters passed is _npassed.

The seq modifier can be used, too, see ?using_parameters .

Also, x_1 and x[1] are two different things. And assume is not the best way to do the assuming. It is better to write assuming n>0 where you need it. That simplifies some things and you don't have to worry about displaying n as n~.

Alec

## normalization...

First, the Fourier series starts with a[0]/2, so you have to correct your a[0],

`a[0]:= 2*op(1, Coefficients);`

Then, Z should be calculated with full normalization,

`Z := FourierTransform(<seq(evalf(f(i*dt)),i=0..N-1)>, normalization=full);`

After that, the approximations of a[i] and b[i] for real f can be defined as

```a0:=2*Z[1];
A:=Array(1..trunc((N-1)/2),datatype=float[8]):
ArrayTools:-Copy(trunc((N-1)/2),2*Re(Z),1,A);
B:=Array(1..trunc((N-1)/2),datatype=float[8]):
ArrayTools:-Copy(trunc((N-1)/2),-2*Im(Z),1,B):
```

For not necessarily real f, a0 can be defined the same, but the a[i] and b[i] coefficients are approximated for i from 1 to floor((N-1)/2) by

```a[i] = Z[i+1]+Z[-i];

b[i] = (Z[i+1]-Z[-i])*I;
```

Alec

## Jacobian...

Jacobian produces a Matrix, and you need a Vector, which can be obtained from it as

`Dwx := Jacobian(x(p,w),[w])[..,1];`

diff and D[2] should work assuming that you executed with(VectorCalculus) somewhere earlier, if you specify the (Maple's) dot product explicitely,

`Dwx.x(p,w)^%T;`

Alec

## add vs. sum and option remember...

Also, recursive procedures like this one work much faster with option remember,

```y:=proc(n::nonnegint) option remember;
if n=0 then a else add(y(k),k=0..n-1) fi
end:

y(4);

8 a

y(1000);

5357543035931336604742125245300009052807024058527668037218751941\
8517552556246806124659918940784792906379733645877657341259\
3572642846157021799228878734928740196728388741211549271053\
7302531185570938977091076523237491790970633699383779582771\
9730385314572855982388432710838302149158263121934186028340\
34688 a

ifactor(op(1,%));
999
(2)
```

Alec

## selecting elements with positive imagina...

For example,

```f:=(b,q)->select(x->Im(evalf(eval(x,[beta=b,Q=q])))>0,t1):

f(1,1);

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  - 8 I Pi  beta Q|
[-|------------------------------------------------|   , I,
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  + 8 I Pi  beta Q|
|------------------------------------------------|   ]
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /

f(-1,1);

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  - 8 I Pi  beta Q|
[|------------------------------------------------|   ,
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  + 8 I Pi  beta Q|
-|------------------------------------------------|   , I]
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /

f(0,1);

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  - 8 I Pi  beta Q|
[|------------------------------------------------|   , I,
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /

/     2       2         2  2   2         2       \1/2
|-4 Pi  - beta  + 4 beta  Q  Pi  + 8 I Pi  beta Q|
|------------------------------------------------|   ]
|     2       2         2              2  2   2  |
\ 4 Pi  + beta  - 4 beta  Q Pi + 4 beta  Q  Pi   /
```

For the purpose of evaluating an integral using residues, it is better to select the indices of the poles in the list instead of selecting the poles themselves,

```F:=(b,q)->select(x->Im(evalf(eval(t1[x],[beta=b,Q=q])))>0,[\$1..nops(t1)]):

F(1,1);

[2, 4, 6]

F(-1,1);

[1, 3, 4]

F(0,1);

[1, 4, 6]
```

Alec

## MinimalPolynomial...

In this example, minimal polynomials could be found as follows,

```alias(a=RootOf(_Z^4+_Z+1)):
minpoly:=x->(Sqrfree(PolynomialTools:-MinimalPolynomial(Expand(x) mod 2,4)) mod 2)[2,1,1]:
```

For example,

```minpoly(a^3);
2     3     4
1 + _X + _X  + _X  + _X

Expand(eval(%,_X=a^3)) mod 2;
0
minpoly(a^1000000);
2
1 + _X + _X

Expand(eval(%,_X=a^1000000)) mod 2;

0
```

A more efficient version is

```Minpoly:=proc(x,y:=_X)(Sqrfree(evala(Norm(Expand(y-x) mod 2))) mod 2)[2,1,1]end:

CodeTools:-Usage(seq(Minpoly(a^k,x),k=1..15)):
memory used=0.91MiB, alloc change=0.69MiB, cpu time=0.03s, real time=0.05s

CodeTools:-Usage(seq(minpoly(a^k),k=1..15)):
memory used=3.09MiB, alloc change=2.75MiB, cpu time=0.06s, real time=0.11s

CodeTools:-Usage(seq(Expand(minpol(k)) mod 2,k=1..15)):
memory used=1.51MiB, alloc change=1.19MiB, cpu time=0.05s, real time=0.06s
```

It works as

```Minpoly(a^2+a, x);
2
x + 1 + x
```

The following is more efficient,

```MinPoly:=proc(x,y:=_X)
local A, i;
A:=Matrix(4,datatype=integer[4]);
for i to 4 do
A[i]:=frontend(PolynomialTools:-CoefficientVector,[Expand(x*a^(i-1)) mod 2,a]) od:
(Sqrfree(LinearAlgebra:-Modular:-CharacteristicPolynomial(2,A,y)) mod 2)[2,1,1]
end:

CodeTools:-Usage(seq(MinPoly(a^k,x),k=1..15));
memory used=369.42KiB, alloc change=255.95KiB, cpu time=0.00s, real time=0.00s
```

The difference can be seen better in larger examples,

```alias(a=RootOf(T^100+T^97+T^96+T^93+T^91+T^89+T^87+T^86+T^82+T^81+T^71+T^70+T^67+T^61+
T^60+T^57+T^54+T^53+T^52+T^49+T^48+T^45+T^44+T^42+T^39+T^36+T^33+T^32+T^31+T^29+T^28+T^27+
T^26+T^24+T^23+T^22+T^18+T^17+T^16+T^14+T^13+T^12+T^10+T^8+T^7+T^6+T^3+T+1)):

MinPoly:=proc(x,y:=_X)
local A, i;
A:=Matrix(100,datatype=integer[4]);
for i to 100 do
A[i]:=frontend(PolynomialTools:-CoefficientVector,[Expand(x*a^(i-1)) mod 2,a]) od:
(Sqrfree(LinearAlgebra:-Modular:-CharacteristicPolynomial(2,A,y)) mod 2)[2,1,1]
end:

CodeTools:-Usage(seq(MinPoly(a^k,x),k=1..1000)):
memory used=1.15GiB, alloc change=28.56MiB, cpu time=50.78s, real time=51.30s

CodeTools:-Usage(seq(Minpoly(a^k,x),k=1..1000)):
memory used=9.21GiB, alloc change=0 bytes, cpu time=14.13m, real time=14.43m
```

minpoly is not good for such large degrees,

```minpoly:=x->(Sqrfree(PolynomialTools:-MinimalPolynomial(Expand(x) mod 2,100)) mod 2)[2,1,1]:
minpoly(a);

18     16     13     8     4     2
_X   + _X   + _X   + _X  + _X  + _X  + _X
```

which is obviously wrong.

Now, minpol doesn't seem to run in this case, and a modified version of it which I wrote, was slower than MinPoly and used more memory.

Alec Mihailovs

PS MinPoly can be made even more efficient, working directly in GF,

```alias(a=RootOf(T^100+T^97+T^96+T^93+T^91+T^89+T^87+T^86+T^82+T^81+T^71+T^70+T^67+T^61+
T^60+T^57+T^54+T^53+T^52+T^49+T^48+T^45+T^44+T^42+T^39+T^36+T^33+T^32+T^31+T^29+T^28+T^27+
T^26+T^24+T^23+T^22+T^18+T^17+T^16+T^14+T^13+T^12+T^10+T^8+T^7+T^6+T^3+T+1)):

F:=GF(2,100,op(a)):
z:=F:-input(2):

MinPolyGF:=proc(x,y:=_X)
local A, i;
A:=Matrix(100,[seq([op(F:-`*`(x,F:-`^`(z,i)))],i=0..99)],datatype=integer[4]);
(Sqrfree(LinearAlgebra:-Modular:-CharacteristicPolynomial(2,A,y)) mod 2)[2,1,1]
end:

CodeTools:-Usage(seq(MinPolyGF(F:-`^`(z,k),x),k=1..1000)):
memory used=0.79GiB, alloc change=27.81MiB, cpu time=19.10s, real time=19.23s
```

Another advantage in using MinPolyGF is that

```n:=(2^100-1)/3:
use F in MinPolyGF(z^n,x) end;
2
x  + x + 1
```

while

```MinPoly(a^n,x);
Error, (in mod/Normal/algnum) modp1: degree too large
```

Alec

 1 2 3 4 5 6 7 Last Page 2 of 76
﻿