Carl Love

Carl Love

28010 Reputation

25 Badges

12 years, 290 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

For polynomials in C[x], the distinguishing of roots is done by the index option of RootOf. See help pages ?RootOf and ?RootOf,indexed. An ordering of the complex numbers based on argument and magnitude is used for this purpose. The ordering can be changed (this is described on ?RootOf,indexed). Since it looks like you're a new Maple user, I wouldn't recommend changing the ordering on your own just yet.

Here is an example:

The RootOf command always has a single bound variable, which it changes to _Z regardless of your initial specification. In the case of an insoluble irrecducible polynomial in C[x] with exact coefficients, the allvalues command separates to roots by index.

Rts:= [allvalues](RootOf(x^5+x^2+1, x));

[RootOf(_Z^5+_Z^2+1, index = 1), RootOf(_Z^5+_Z^2+1, index = 2), RootOf(_Z^5+_Z^2+1, index = 3), RootOf(_Z^5+_Z^2+1, index = 4), RootOf(_Z^5+_Z^2+1, index = 5)]

NRts:= <evalf(Rts)>:
DataFrame(<NRts | argument~(NRts) | abs~(NRts)>, columns= [Root, Arg, Abs]);

DataFrame(Matrix(5, 3, {(1, 1) = .7515192324+.7846159210*I, (1, 2) = .8069402580, (1, 3) = 1.086463667, (2, 1) = -.1545896767+.8280741332*I, (2, 2) = 1.755357607, (2, 3) = .8423803999, (3, 1) = -1.193859111, (3, 2) = 3.141592654, (3, 3) = 1.193859111, (4, 1) = -.1545896767-.8280741332*I, (4, 2) = -1.755357607, (4, 3) = .8423803999, (5, 1) = .7515192324-.7846159210*I, (5, 2) = -.8069402580, (5, 3) = 1.086463667}), rows = [1, 2, 3, 4, 5], columns = [Root, Arg, Abs])

 

Download IndexedRootOf.mw

This works:

thaw(PDEtools:-Solve(subs(Ys=~ (Fz:= freeze~(Ys)), sys), Fz));

The variable names that have session-dependent assumptions are the indices of the table `property/OrigName`. Thus the following supplies the information that you're looking for, although not exactly in the format that you showed:

about(indices(`property/OrigName`, 'nolist'));

You can learn where the information is stored by reading the code of about, which is fairly short and straightforward:

showstat(about);

@C_R The fact that map(t-> t*~2, eqs) "works" is just a red herring, and has nothing to do with the elementwise operator, since map(t-> t*2, eqs) also works, producing identical results. Indeed, applying any elementwise operator to an equation (or indeed to any structure other than a setlisttable, rtable, or expression sequence) will give identical results to the non-elementwise form of the operator.

The anomaly with `*` must be due to the lack of a final eval internally, the same thing that you noticed. I can't yet explain the anomaly with `^`, but it might be related to either or both of these:

  • Of the 5 operators you tested, only `^` is non-associative.
  • Of those 5, only `^` requires a specific number of arguments, namely 2.

I don't think that the behavior of the standard arithmetic operators on non-algebraic structures is documented, so I can't call what you discovered a bug; it's just an anomaly. You can change the behavior like this:

AdjustArithOps:= module()
export
    `*`:= eval@:-`*`,

    `^`:= proc(a,b,$)
        try :-`^`(a,b) 
        catch "non-algebraic base": 
            if a::atomic then error else map(thisproc, a, b) fi
        catch "non-algebraic exponent": 
            if b::atomic then error else map2(thisproc, a, b) fi
        end try
    end proc
;
end module
:
eqs:= [a=b, c=d]:
use AdjustArithOps in 
    #Now you can do what you tried in your Question, or you can 
    #simply use the elementwise forms:
    map(`^`, eqs, 2); eqs^~2; map(`*`, eqs, 2); eqs*~2
end use; 
                       [ 2    2   2    2]
                       [a  = b , c  = d ]

                       [ 2    2   2    2]
                       [a  = b , c  = d ]

                     [2 a = 2 b, 2 c = 2 d]

                     [2 a = 2 b, 2 c = 2 d]

 

Yes, your computations can be greatly simplified by using elementwise operations (see help page ?elementwise).

restart

C := x^2/A^2+y^2/B^2-1

x^2/A^2+y^2/B^2-1

L := y = k*x+b

y = k*x+b

#option "explicit" not necessary, but might help for more complicated equations:
x||(1..2):= solve(subs(L,C), x, explicit);

-(b*k*A-(A^2*B^2*k^2+B^4-B^2*b^2)^(1/2))*A/(A^2*k^2+B^2)

-(b*k*A+(A^2*B^2*k^2+B^4-B^2*b^2)^(1/2))*A/(A^2*k^2+B^2)

L||(1..2):= subs~([y=~ y||(1..2)], [x=~ x||(1..2)], L)[];

y1 = -k*(b*k*A-(A^2*B^2*k^2+B^4-B^2*b^2)^(1/2))*A/(A^2*k^2+B^2)+b

y2 = -k*(b*k*A+(A^2*B^2*k^2+B^4-B^2*b^2)^(1/2))*A/(A^2*k^2+B^2)+b

y1_plus_y2, y1_by_y2 := simplify(`~`[rhs]([L1+L2, L1*L2]))[]

2*B^2*b/(A^2*k^2+B^2), B^2*(-A^2*k^2+b^2)/(A^2*k^2+B^2)

NULL

Download Elementwise.mw

The possibility of a programmatic or automatic solution to this has been mentioned, but no code actually given. This is easier to code than one may expect.

restart
:

IntofSum2SumofInt:= IofS->
    evalindets(
        IofS,
        And(specfunc({Int,int}), specfunc({sum,Sum}) &under curry(op,1))
            &under combine,
        e-> local Op:= rcurry(op, combine(e))@`[]`;
            Op(1,0)(Op(0)(Op(1,1), Op(2..)), Op(1, 2..))
    )
:

inte_eq:= int(
    -(sum(B[i]*beta[i]*exp(-beta[i]*(t-tau))*tau, i= 1..n)),
    tau= 0..t
);

int(-(sum(B[i]*beta[i]*exp(-beta[i]*(t-tau))*tau, i = 1 .. n)), tau = 0 .. t)

IntofSum2SumofInt(inte_eq);

sum(-B[i]*(beta[i]*t+exp(-beta[i]*t)-1)/beta[i], i = 1 .. n)

 

Download IntofSum.mw

I don't know why the prettyprintred integrals and sums didn't appear above, but you can see them in the worksheet.

You could replace the extra factor of tau with 1/tau in the above, but the result is trivial because the integrals diverge at tau=0 (unless the Bs and betas make all the integrands zero).

Using assign is usually a bad idea. It can usually be replaced by eval, as is the case here. Your code can be shortened to this:

L:= [3,4,x,x^2,x^3]:
((X,x::algebraic)-> local n, a;
    `if`(patmatch(X, x^n::'nonunit'(anything), a), eval(n,a), NULL)
)~(L,x);
                             [2, 3]

You wanted to make the procedure shorter, yet you seem to go out of your way to make your code unnecessarily verbose by using RETURNmap instead of ~, etc.

Even shorter is

op~(2, select(patmatch, L, x^n::'nonunit'(anything)));
 

Yes, the way to do it is to create the types in a ModuleLoad. Here is an example:

M:= module()
local
    AddType:= proc(T::name, H::{type, procedure})
    uses TT= TypeTools;
        if TT:-Exists(T) then TT:-RemoveType(T) fi;
        TT:-AddType(T,H)
    end proc,

    ModuleLoad:= proc()
        AddType(Point, (e, n::{2,3})-> type(e, [algebraic$n]));
        AddType(Line, (e, n::{2,3})-> type(e, [Point(n)$2]))
    end proc
;
export
    CreateLine:= overload([
        proc(Pts::Line(2), $) option overload; "2D line" end proc,
        proc(Pts::Line(3), $) option overload; "3D line" end proc,
        ()-> "Unrecognized input format"
    ])
;
    ModuleLoad()
end module
:
M:-CreateLine([[1,2,3],[4,5,6]]);
                           "3D line"

M:-CreateLine([[1,2,3],[4,5]]);
                  "Unrecognized input format"

 

You almost had it correct with your subs. To completely eliminate the integral, the diff(eq2, x) should also be a target of the substitution: 

restart
:
eq:= diff(y(x),x) = y(x) + exp(x)*exp(-3*x)/2 + int(exp(x+t)*y(t), t= 0..x):
IC:= y(0)=1: 

eq2:= g(x) = int(exp(x+t)*y(t), t= 0..x):
IC_2:= eval(eq2, x= 0):  #Generate 2nd IC automatically.

#Eliminate the integral:
sys:= subs((rhs=lhs)(eq2), {eq, diff(eq2,x), IC, IC_2});
        / d                               
sys := { --- g(x) = g(x) + exp(2 x) y(x), 
        \ dx                              

   d                1                                            
  --- y(x) = y(x) + - exp(x) exp(-3 x) + g(x), g(0) = 0, y(0) = 1
   dx               2                                            

  \ 
   }
  / 

Now sys can be solved symbolically with dsolve(sys) or numerically with dsolve(sys, numeric).

The idea of factoring out the exp(x) will also work, but the above is simpler: Just make a substitution for the entire original integral.

This code does essentially the same thing as @dharr 's but does it with the index-free and functional syntax that I prefer. The entire computation---the construction of the Array---is done by a single line of code. The remaining two lines are simply to put that Array into a neatly displayed Table.

restart
:

Given: Two  leagues (lists of teams). It isn't necessary that the teams or the leagues have the same number of members, althoughh that

happens to be true in this case.

new:= [
    ["Pablo G", "Lucy Z", "Matt C"], ["John S", "Rosalie B", "Morgan B"],
    ["Daniel N", "Cohen M", "Steve S"], ["Thomas S", "Peter M", "Tony W"],
    ["Susanne A", "Stephen J", "Jorah S"], ["Sapan B", "Rishab S", "Jesse B"],
    ["Dean B", "Jesse W", "Hak Y"]
]:
past:= [
    ["Jesse W", "Daniel N", "Lucy Z"], ["John S", "Jesse B", "Rosalie B"],
    ["Jorah S", "Peter M", "Thomas S"], ["Stephen J", "Pablo G", "Matt C"],
    ["Rishab S", "Cohen M", "Steve S"], ["Tony W", "Dean B", "Morgan B"],
    ["Susanne A", "Hak Y", "Sapan B"]
]:

Goal: Find all intersections of new teams with past teams.

#Array of the intersections of the pairs of the Cartesian product of the 2 leagues.
#The raw Array is more suitable than the displayed Table for programmatic access:
new_Xinter_past:= ArrayTools:-GeneralOuterProduct({op}~(new), op@`intersect`, {op}~(past));

"[[["Lucy Z",,,"Matt C","Pablo G",,,],[,"John S","Rosalie B",,,,"Morgan B",],["Daniel N",,,,"Cohen M","Steve S",,],[,,"Peter M","Thomas S",,,"Tony W",],[,,"Jorah S","Stephen J",,,"Susanne A"],[,"Jesse B",,,"Rishab S",,"Sapan B"],["Jesse W",,,,,"Dean B","Hak Y"]]]"

#Display as a Table. The middle operator (()-> ...) of the @-construct is
#for the elegant display of Table entries, without extraneous punctuation:
(Tabulate @ (()-> if nargs=0 then else [args, if nargs=1 then "" else fi] fi)~ @ DataFrame)(
    new_Xinter_past, rows= ["new "||($nops(new))], columns= ["past "||($nops(past))]
):

Download PairwiseIntersection.mw

The Table didn't transfer to MaplePrimes, but you can see it in the worksheet. It's identical to @dharr's both in content and displayed form.

You ( @Christopher2222 ) are a Moderator, and thus you can read the deleted content. If you do this, you will see that the vast majority of Moderator-deleted (as opposed to author-deleted) content is spam.

You are wrong about lost knowledge: It is the failure to delete inappropriate content that will cause the greater loss of knowledge, because good writers will stop using MaplePrimes, and it will become irrelevant. Do you think periodicals shouldn't have editors and peer review?

The defined coordinate systems are stored in an unexported table of an  undocumented package, as given in the title. Here's a procedure to access it:

Coords:= proc(C::{name,function}, V::list(name))
uses PC= Plot:-CoordinateSystems;  #undocumented package
local
    opq:= kernelopts(opaquemodules= false),    
    r:= PC:-CoordinateSystemsTable[String(`if`(C::name, C, op(0,C)))](
        V[],`if`(C::name, [][], op(C))
    )
;
    kernelopts(opaquemodules= opq);
    r
end proc
:
Coords(spherical_physics, [r, theta, phi]);
      [r*cos(phi)*sin(theta), r*sin(theta)*sin(phi), r*cos(theta)]

Since the package is undocumented, don't rely too much on the syntax remaining exactly the same as the way I just accessed it in the procedure.

By using a Newton's method iteration with pure-integer builtin (i.e., GMP) arithmetic, I can find roots of integers significantly faster than with any method that uses evalf. While most of Maple's integer-arithmetic commands whose names begin with i are builtin, that isn't true of iroot: It's written in Maple, and it uses evalf.

My Newton procedure is the second procedure listed in the table methods in the worksheet below. The Newton iteration for finding the cube root r of is simply

r-> (2*r^3 + n)/(3*r^2)

restart
:

#Make garbage collection single threaded, because multi-threaded makes
#"real-time" comparisons difficult:
kernelopts(gcmaxthreads= 1)
:

#
# The procedures to be tested and compared
# ----------------------------------------
methods:= table([
    #iroot with floor correction:
    Iroot= (n-> local r:= iroot(n,3); `if`(r^3 > n, r-1, r)),

    #Integerized Newton's method:
    Newton= proc(n)
    local r1:= integermul2exp(1, trunc(ilog2(1+n)/3 + 1/2)), r, r2;
        do r1:= iquo(2*(r2:= (r:= r1)^2)*r + n, 3*r2) until r-r1 in {-1,0,1};
        min(r,r1)
    end proc,

    #Two of @dharr's methods:
    root@evalf= (n-> floor(root(evalf(n),3))),
    evalf@root= (n-> floor(evalf(root(n,3)))),

    #iroot's own method, minimalized:
    evalf@`^`= (n-> trunc(evalf(n^(1/3))))
]):
meth_names:= [indices](methods, nolist);

[`@`(evalf, `^`), `@`(root, evalf), Iroot, Newton, `@`(evalf, root)]

# Generate random test cases:
#
nWords:= 3:  iterations:= 2^16:
R:= rand(2^((nWords-1)*64)..2^(nWords*64)-1):
#Arbitrary key for consistent test data on different runs and platforms.
#Change the key to generate different random test data.
randomize(23):
testdata:= table(
    [meth_names[], accuracy]=~ ['['R'()$iterations]' $ 1+numelems(methods)]
):

#Set sufficient floating-point precision (see showstat(iroot), line 39)
#(only needed for methods that directly use evalf):
d:= length(2^(nWords*64)-1):  Digits:= iquo(d,3)+5+length(d);

26

# The tests:
#
#Accuracy test: If all methods agree, then all will be listed in a single set:
entries(ListTools:-Classify(j-> methods[j]~(testdata[accuracy]), meth_names), nolist);

{Iroot, Newton, `@`(evalf, `^`), `@`(evalf, root), `@`(root, evalf)}

#Efficiency test:
for j in meth_names do
    f:= methods[j];
    printf("%a:\n", j);
    CodeTools:-Usage(seq[reduce= ()](f(n), n= testdata[j]));
    printf("\n")
od:

evalf@`^`:
memory used=1.59GiB, alloc change=-16.00MiB, cpu time=4.27s, real time=7.89s, gc time=750.00ms

root@evalf:
memory used=1.79GiB, alloc change=0 bytes, cpu time=3.69s, real time=9.42s, gc time=578.12ms

Iroot:
memory used=485.51MiB, alloc change=0 bytes, cpu time=953.00ms, real time=1.96s, gc time=218.75ms

Newton:
memory used=211.56MiB, alloc change=0 bytes, cpu time=750.00ms, real time=1.42s, gc time=203.12ms

evalf@root:
memory used=3.62GiB, alloc change=0 bytes, cpu time=13.97s, real time=25.77s, gc time=1.73s
 

Newton's method is the winner.

 

The Newton-method procedure for any particular root be generated by this procedure:

NewtonIroot:= (k::And(posint, Not(1)))-> subs(
    _k= k,
    proc(n::nonnegint)
    local r1:= integermul2exp(1, trunc(ilog2(1+n)/_k + 1/2)), r, r2;
        do r1:= iquo((_k-1)*(r2:= (r:= r1)^(_k-1))*r + n, _k*r2)
        until r-r1 in {-1,0,1};
        min(r,r1)
    end proc
):

showstat((Iroot5:= NewtonIroot(5)));


proc(n::nonnegint)
local r1, r, r2;
   1   r1 := integermul2exp(1,trunc(1/5*ilog2(1+n)+1/2));
   2   do
   3       r1 := iquo(4*(r2 := (r := r1)^4)*r+n,5*r2)
       until r-r1 in {-1, 0, 1};
   4   min(r,r1)
end proc
 

r:= rand();

441235504072

n:= (r+1)^5-1:  Iroot5(n);

441235504072

 

Download Iroot.mw

The votes awarded after deletion do apply to the author's reputation. It has always been true (at least for the 10 years or so that I've been using MaplePrimes) that votes awarded to an item which is subsequently deleted continue to apply to the author's reputation.

I totally agree with acer that items Deleted As Spam should not be retained. There are now several Posts where the majority of Replies have been deleted as spam yet remain as clutter.

The syntax of Maple's hypergeom always has exactly three arguments. The first two arguments are arbitrary-length lists of algebraic expressions, and the third is a single algebraic expression rather than a list. If it's a 2F1 function, then the first list has 2 elements, and the second has 1.

The biggest and ugliest (IMO) differences between Maple and Mathematica syntax are that Maple's lists are enclosed in square brackets ([...]) rather than curly braces ({...}) and that Maple's function arguments are always in parentheses ((...)) rather than square brackets. (I consider the Mathematica to be uglier.)

In another Answer, @janhardo shows the correct syntax of a simple Maple hypergeom function (of class 2F1): 

hypergeom([1,1], [2], -z)

Your example could be entered as

-r^2/3*((1-beta)/`&gamma;`)^n*
    hypergeom(
        [n, n/(1-beta)], [(beta-1-n)/(beta-1)], -r^(3/n*(1-beta))*chi^2/`&gamma;`
    )

I used `&gamma;` instead of gamma because gamma is by default a pre-defined constant (similar to Pi). If you mean the constant, then just use gamma. If it's intended to be a variable or parameter, there are several workarounds, one of which is `&gamma;`.

2 3 4 5 6 7 8 Last Page 4 of 394