Carl Love

Carl Love

28100 Reputation

25 Badges

13 years, 103 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@vv Has the posted code been changed? I can find no use of indices in either of the posted worksheets. I searched with the GUI's Find command (from the Edit menu). I searched both inside and outside of the Code Edit Regions.

@YasH Your test with randomize is very strong evidence that the variability of the order of the returned eigenvectors is not due to the code's use of some type of randomization. So, the variability must be due to some address-based storage issue, such as what's discussed about the indices command above. I use indices in nearly all of my programs where the order of the returned results doesn't matter because it's an extremely efficient command, so I'd guess that it's widely used in the library code also.

@vv Yes, indices are stored by some address scheme, and hence the order is session dependent. If you want the results of indices or entries to be returned in a consistent order, then use option indexorder:

indices(MyTable, indexorder);

and likewise for entries. The order will then be the same as if the indices (not the entries) were an ordinary Maple set.

@YasH That's a good question, and, yes, it's related. I don't know the answer, but it's not too difficult to figure out with a few experiments. For example, after a restart, do the four always come in the same order? If yes, that's a sign that random numbers are used in the algorithm. If no, that's a sign that the ordering is based on the addresses of certain objects, most likely vectors.

The idea that random numbers might be used in the algorithm can be further investigated by checking if the pre-determined sequence of random numbers changes after calling Eigenvectors.

@vv Vote up. My tests below show that a simple for-loop approach, taking the remainder at each step, is fastest, at least for word-sized integers.

In a post a few months ago, I showed that the divide-and-conquer approach was much faster than plain mul or a linear for-loop for plain multiplication. Apparently, plain mul uses a linear simple-loop approach. This probably prevents the use of the much faster FFT-based multiplication algorithm (which is only used when both factors are supra-word-length integers).

Taking the remaindering into account, the divide-and-conquer approach is still much faster than mul, but it doesn't beat the simple for loop. ModMul.mw

restart:


#Divide-and-conquer multiplication; no mod
DandCMul:= proc(X)
local n:= nops(X), h;
     if n < 4 then `*`(X[])
     else
          h:= iquo(n,2);
          thisproc(X[..h])*thisproc(X[h+1..])
     end if
end proc:


#Divide-and-conquer and do mod at each step
DandCModMul:= proc(X, p)
local n:= nops(X), h, r;
     if n < 4 then r:= `*`(X[])
     else
          h:= iquo(n,2);
          r:= thisproc(X[..h],p)*thisproc(X[h+1..],p)
     end if;
     irem(r,p)
end proc:
 


#This should be similiar to `fold` but a little faster and less garbage.
LinearModMul:= proc(X, p)
local x, r:= 1;
     for x in X do r:= irem(r*x, p) end do
end proc:
 

 

Test with word-length integers and modulus.

 

R:= rand(1..2^64):
L:= ['R()'$2^17]:

p:= nextprime(R()):


gc4:= proc() to 4 do gc() od end proc:


gc4(): r1:= CodeTools:-Usage(irem(mul(x, x= L), p)):

memory used=62.57GiB, alloc change=23.49MiB, cpu time=68.28s, real time=62.95s, gc time=49.62s

gc4(): r2:= CodeTools:-Usage(DandCModMul(L,p)):

memory used=64.66MiB, alloc change=0 bytes, cpu time=594.00ms, real time=607.00ms, gc time=93.75ms

gc4(): r3:= CodeTools:-Usage(LinearModMul(L, p)):

memory used=15.52MiB, alloc change=0 bytes, cpu time=250.00ms, real time=202.00ms, gc time=109.38ms

gc4(): r4:= CodeTools:-Usage(irem(DandCMul(L),p)):

memory used=86.11MiB, alloc change=0 bytes, cpu time=641.00ms, real time=636.00ms, gc time=93.75ms

evalb~([seq](r1=r||k, k= 2..4));

[true, true, true]

 

 


Download ModMul.mw

@vv Yes, my tests show that you are correct about this.

@taro The semicolon at the end of the added line needs to be a comma, and the procedure needs to be surrounded by parentheses:

Change this:

numer_expand := a->expand(numer(a))/denom(a);

to this:

numer_expand := (a->expand(numer(a))/denom(a)),

@sigl1982 Presumably your original code had some condition by which it would decide when it was time to stop going to Step8, or else it'd be an infinite loop. Since you didn't show that condition, I just made up the name condition0 for it. It's not the same as continue or not continue.

@sigl1982 A simple control variable (I called it continue below) should do it:

while condition0 do
      #Set up parameters.
      continue:= true;
      while condition and continue do
          for s in S while continue do
              #  do anything
              if condition2 then   
                 #  do anything
                 if condition3 then
                       #  do anything
                       continue:= false;   
                 end if;
              end if;  
          end do;    
      end do;
end do;

Please explain in more detail how you can do it manually. Maybe I'd be able to use that info to figure out how to do it with Maple.

@Kitonum But this is not even close to the OP's proposed result.

@sand15athome The following may be already totally obvious to you. But, just in case you don't know...

PLOT isn't a command of the Maple language; it's just an inert data structure that's interpreted by the GUI. (So, perhaps it can be considered a command to the GUI, but to the Maple kernel it's inert.) All commands that produce static (i.e., non-animated) 2-D plots have this data structure as their output. As long as you're aware of all that, I think that it's fine to write your own procedures that produce PLOT output; indeed, I often find it necessary.

@TomM Good work investigating this.

I made your ideas into a procedure. I use randomly selected complex evaluation points selected over a square centered at the origin of user-specified dimension (defaults to 4x4).

restart:

DetUnifloat:= proc(
     M::Matrix(square, polynom(float)),
     {radius::positive:= 2},
     {discardthreshold::positive:= Float(1, 3-Digits)}
)
description
     "A replacement for LinearAlgebra:-Determinant(..., method= unifloat)"
;
local
     R:= rand(evalf(-radius..radius)),
     x:= indets(M, And(name, Not(constant))),
     n:= op([1,1], M),
     maxdeg, mindeg, i, j, E, P, MC;
;
     if nops(x) > 1 then
          error "Only one variable allowed"
     end if;
     x:= `if`(nops(x) > 0, x[], 'x');
     maxdeg:= min(
          add(max(degree~(M[i,..])), i= 1..n),
          add(max(degree~(M[..,j])), j= 1..n)
     );
     mindeg:= max(
          add(min(ldegree~(M[i,..])), i= 1..n),
          add(min(ldegree~(M[..,j])), j= 1..n)
     );
     E:= ['R()' + 'R()'*I $ maxdeg+1];
     P:= CurveFitting:-PolynomialInterpolation(
          E, [seq](LinearAlgebra:-Determinant(eval(M, x= e)), e= E), x
     );
     MC:= max(abs~([coeffs](P,x)));
     select(
          t-> degree(t,x) >= mindeg,
          sort(MC*simplify(fnormal(P/MC, Digits, discardthreshold), 'zero'))
     )
end proc:
     

data_deg_ex := Matrix(6, 6, [[180*sigma2^2, -17980.4676072929*sigma2, 60*sigma2^2, -1792.74593023400*sigma2, -597.004097753022*sigma2, -60*sigma2], [-17980.4676072929*sigma2, 60*sigma2^2, -17980.4676072929*sigma2, -597.004097753022*sigma2, -597.581976744666*sigma2, 5993.48920243096], [60*sigma2^2, -17980.4676072929*sigma2, 180*sigma2^2, -597.581976744666*sigma2, -1791.01229325907*sigma2, -60*sigma2], [-1792.74593023400*sigma2, -597.004097753022*sigma2, -597.581976744666*sigma2, -60*sigma2, 5993.48920243096, 597.581976744666], [-597.004097753022*sigma2, -597.581976744666*sigma2, -1791.01229325907*sigma2, 5993.48920243096, -60*sigma2, 597.004097753022], [-60*sigma2, 5993.48920243096, -60*sigma2, 597.581976744666, 597.004097753022, 60]]);

data_deg_ex := Matrix(6, 6, {(1, 1) = 180*sigma2^2, (1, 2) = -17980.4676072929*sigma2, (1, 3) = 60*sigma2^2, (1, 4) = -1792.74593023400*sigma2, (1, 5) = -597.004097753022*sigma2, (1, 6) = -60*sigma2, (2, 1) = -17980.4676072929*sigma2, (2, 2) = 60*sigma2^2, (2, 3) = -17980.4676072929*sigma2, (2, 4) = -597.004097753022*sigma2, (2, 5) = -597.581976744666*sigma2, (2, 6) = 5993.48920243096, (3, 1) = 60*sigma2^2, (3, 2) = -17980.4676072929*sigma2, (3, 3) = 180*sigma2^2, (3, 4) = -597.581976744666*sigma2, (3, 5) = -1791.01229325907*sigma2, (3, 6) = -60*sigma2, (4, 1) = -1792.74593023400*sigma2, (4, 2) = -597.004097753022*sigma2, (4, 3) = -597.581976744666*sigma2, (4, 4) = -60*sigma2, (4, 5) = 5993.48920243096, (4, 6) = 597.581976744666, (5, 1) = -597.004097753022*sigma2, (5, 2) = -597.581976744666*sigma2, (5, 3) = -1791.01229325907*sigma2, (5, 4) = 5993.48920243096, (5, 5) = -60*sigma2, (5, 6) = 597.004097753022, (6, 1) = -60*sigma2, (6, 2) = 5993.48920243096, (6, 3) = -60*sigma2, (6, 4) = 597.581976744666, (6, 5) = 597.004097753022, (6, 6) = 60})

DetUnifloat(evalf(data_deg_ex));

HFloat(1.8662406249298306e11)*sigma2^8+HFloat(1.4795511563262647e14)*sigma2^7+HFloat(4.0379217068100504e16)*sigma2^6+HFloat(4.393918722504788e18)*sigma2^5+HFloat(1.6722907322306802e20)*sigma2^4

LinearAlgebra:-Determinant(data_deg_ex, method= minor);

186624000000*sigma2^8+147955115634640.*sigma2^7+0.403792170686782e17*sigma2^6+0.439391872250447e19*sigma2^5+0.167229073223072e21*sigma2^4

 


Download unifloat.mw

@tomleslie I think that you may be right---that unlike dsolve, pdsolve's solution procedures don't contain the (correct) partial derivatives. However, Maple says that the requested derivative is 0---it doesn't produce an empty plot or give an error message. It also gives 0 if you use the pds:-value method.

@vv And the value of that Stieltjes integral is not 0?

First 379 380 381 382 383 384 385 Last Page 381 of 709