Alec Mihailovs

Dr. Aleksandrs Mihailovs

4415 Reputation

20 Badges

17 years, 32 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

Social Networks and Content at Maplesoft.com

Maple Application Center

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.

MaplePrimes Activity


These are answers submitted by Alec Mihailovs

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

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

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

Alec

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

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

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,

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

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

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

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

In Maple, sum is for symbolic summation. For just adding numbers, add should be used instead.

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

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

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