Carl Love

Carl Love

28085 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I changed your example procedure to one where the order of the arguments matters so that the example is more generally applicable.

myproc:= (a,b,c)-> a/b-c:

 

The `~` command below applies the procedure in "batch" form.

(myproc@op) ~ ([[d,e,f], [g,h,i]]);

Because of operator precedences, all the parentheses in the above command are necessary, unfortunately. Also, the arguments must be in lists [] rather than sets {} because the order matters. I think this command is fairly close to your desired syntax.

The above command is equivalent to the more traditional form

map(myproc@op, [[d,e,f], [g,h,i]]);

 

The op is used to convert lists, such as [d,e,f], to sequences, such as d, e, f; since myproc expects a sequence of arguments.
                      

Example 1: Factor the first 29 odd 27-bit integers. First, the "old" way:

restart;
N:= 2^(2^7-1)+1:
B:= [seq](N+k, k= 0..2^9-1, 2):
CodeTools:-Usage(map(ifactor, B)):

memory used=1.12GiB, alloc change=212.68MiB, cpu time=18.76s, real time=21.10s

...and now using package Threads:

restart;
N:= 2^(2^7-1)+1:
B:= [seq](N+k, k= 0..2^9-1, 2):
CodeTools:-Usage(Threads:-Map(ifactor, B)):

memory used=1.23GiB, alloc change=0.53GiB, cpu time=33.50s, real time=8.97s

 

Example 2: Generate 214 random 27-bit primes. The old way:

restart;
N:= 2^(2^7-1):
Gen:= RandomTools:-MersenneTwister:-NewGenerator(range= N..2*N):
seeds:= [seq](Gen(), k= 1..2^14):
Primes:= CodeTools:-Usage(map(nextprime, seeds)):

memory used=2.76GiB, alloc change=327.83MiB, cpu time=27.50s, real time=27.50s

Primes[-1];  #Inspect last prime in list...

            197727613027066792031910829881298775459

...and using Threads:

restart;
N:= 2^(2^7-1):
Gen:= RandomTools:-MersenneTwister:-NewGenerator(range= N..2*N):
seeds:= [seq](Gen(), k= 1..2^14):
Primes:= CodeTools:-Usage(Threads:-Map(nextprime, seeds)):

memory used=3.74GiB, alloc change=1.56GiB, cpu time=77.12s, real time=11.00s

Primes[-1];  #...just to verify that the two methods are doing the same computation.

            197727613027066792031910829881298775459

Records are mutable, as the following shows:

restart;
A:= Record(a=1):
B:= Record(a=1):
T[A]:= 1: T[B]:= 2:
eval(T);

   table([B = 2, A = 1])

Compare that with what happens with an immutable structure, in this case an unevaluated function call (named `record`, just for fun):

restart;
A:= record(a=1):
B:= record(a=1):
T[A]:= 1: T[B]:= 2:
eval(T);

   table([record(a = 1) = 2])

But that doesn't mean that they absolutely can't be used with remember tables. If the records being compared are address identical, not just copies which "just happen to be" equal, then the remember table will recognize them. If just a few cases slip through unrecognized because they just happen to be equal, it won't ruin your recursion. On the other hand, if your algorithm relies on making copies of records (not just copying the address, or pointer), the remember table won't work.

Cache tables, remember tables, and just-plain tables all function like the tables in the examples above. The only difference that cache tables have is that the amount of memory that they use can be easily controlled.

restart;
`convert/range`:= proc(ineqs::{set,list})
# Convert a pair of simple inequalities to form x= a..b.
   local ab:= indets(ineqs, numeric);
   indets(ineqs, name)[]= min(ab)..max(ab)
end proc:
n:= 1..3:
s1:= (x-1)*(y-x):
s2:= (7-y)*(1-x):
s3:= (x-y)*(y-7):
feas:= [solve](({s||n} >=~ 0) union {-2<=x,x<=3, 0<=y,y<=11}, {x,y});

      [{y = x, 0 <= x, x <= 1}, {y = 7, 1 <= x, x <= 3},
        {x = 1, 1 <= y, y <= 7}]

z:= add(sqrt(s||i), i= n):
m:= 1..nops(feas):
zr||m:= seq(eval(z, indets(feas[i], `=`)), i= m):
seq(maximize(zr||i, convert(indets(feas[i], `<=`), range)), i= m);

                      (1/2)   (1/2)   (1/2)
                     7     , 8     , 9   

So the answer is sqrt(9), or 3.

In your example, you don't need to create the Array until you know how many entries it will need. Presumably blocks is a list with the same number of entries as you want in the Array. That number can be accessed as nops(blocks). So right after the line where you create blocks, you can put the line

c:= Array(1..nops(blocks));

You can eliminate the for loop by including an initializer in the above statement, and putting it after the definitions for e and m:

c:= Array(1..nops(blocks), i-> fme(convert(blocks[i], bytes), e, m));

and not printing c until after the procedure returns.

Here's my solution:

> restart;
    invtau:= proc(n::posint)
   # Inverse tau: Returns the smallest positive integer r such
   # tau(r) = n, where tau(r) is the number of positive integer
   # divisors of r (see ?numtheory,tau).
       local pair, r:= 1, r_prime:= 1, exponent;
       for pair in ListTools:-Reverse(sort(ifactors(n)[2])) do
           exponent:= pair[1] - 1;
           to pair[2] do
               r_prime:= nextprime(r_prime);
               r:= r*r_prime^exponent
           od
       od;
       r
  end proc:

> st:= time():
> [seq]([n,invtau(n)], n= 1..100);


[[1, 1], [2, 2], [3, 4], [4, 6], [5, 16], [6, 12], [7, 64], [8, 30], [9, 36], [10, 48], [11, 1024], [12, 60], [13, 4096], [14, 192], [15, 144], [16, 210], [17, 65536], [18, 180], [19, 262144], [20, 240], [21, 576], [22, 3072], [23, 4194304], [24, 420], [25, 1296], [26, 12288], [27, 900], [28, 960], [29, 268435456], [30, 720], [31, 1073741824], [32, 2310], [33, 9216], [34, 196608], [35, 5184], [36, 1260], [37, 68719476736], [38, 786432], [39, 36864], [40, 1680], [41, 1099511627776], [42, 2880], [43, 4398046511104], [44, 15360], [45, 3600], [46, 12582912], [47, 70368744177664], [48, 4620], [49, 46656], [50, 6480], [51, 589824], [52, 61440],
[53, 4503599627370496], [54, 6300], [55, 82944], [56, 6720], [57, 2359296], [58, 805306368], [59, 288230376151711744], [60, 5040], [61, 1152921504606846976], [62, 3221225472], [63, 14400], [64, 30030], [65, 331776], [66, 46080], [67, 73786976294838206464], [68, 983040], [69, 37748736], [70, 25920], [71, 1180591620717411303424], [72, 13860], [73, 4722366482869645213696], [74, 206158430208], [75, 32400], [76, 3932160], [77, 746496], [78, 184320], [79, 302231454903657293676544], [80, 18480], [81, 44100], [82, 3298534883328], [83, 4835703278458516698824704], [84, 20160], [85, 5308416], [86, 13194139533312], [87, 2415919104], [88, 107520], [89, 309485009821345068724781056], [90, 25200], [91, 2985984], [92, 62914560], [93, 9663676416], [94, 211106232532992], [95, 21233664], [96, 60060], [97, 79228162514264337593543950336], [98, 233280], [99, 230400], [100, 45360]]

> time()-st;
                             0.031

A:= x^n+y^3.5:
                     
indets(A, anything^float);
                             / 3.5\
                            { y    }
                             \     /
indets(A, anything^symbol);
                                n
                             { x  }
                             
indets(A, `^`);
                             n   3.5
                          { x , y    }
                                   

There is a way to solve this without resorting to commands outside of dsolve. The help page ?dsolve,numeric_bvp,advanced gives advice on how to apply dsolve's options to resolve some error messages given by `dsolve/numeric/bvp`. For messages related to Newton iteration convergence, it suggests the 'continuation' option, which is defined on help page ?dsolve,numeric,bvp (with two examples of its use given on the former help page).

To use the method, one introduces a parameter λ into the problem such that λ can vary continuously from 0 to 1, with λ=0 yielding a relatively easier BVP and λ=1 giving the desired BVP. The parameter can appear anywhere in the problem: in the ODEs, or the boundary conditions, or both. There is some subtlety required in deciding where and how to place λ. I tried this two ways with the problem at hand, and they both worked.

My first idea was to have λ=0 eliminate the convolution of f(x) and g(x), so I applied it to the g(x) term in the first ODE and to the term containing f(x) in the second ODE.

restart;
ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 ) (x) + lambda*g(x) = 0;
ode2:= (D@@2)(g)(x) + lambda*3*10*(f*D(g))(x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = 1, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color= black);

I recommend plotting the interval 5.7 .. 6 separately because there is an interesting fluctuation there which is not visible at the g-axis scale of the overall plot.

My second idea was to use λ=0 to "level" the boundary conditions so that everything starts and ends at 0 on the vertical axis. In this case, that simply means changing g(0)=1 to g(0)=λ.

restart;
ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 + g ) (x) = 0;
ode2:= ( (D@@2)(g) + 3*10*f*D(g) ) (x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = lambda, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color = black);

Same plot, of course.

You are mixing new and old matrix commands. The new commands all begin with capital letters and tend to be spelled in full rather than abbreviated. Use Matrix and Vector instead of array. Use Adjoint instead of adjoint. Use Determinant instead of Det. See ?Matrix ?Vector ?Determinant etc.

I don't know why this unanswered question from seven years ago is appearing near the top of the "Active Conversations" list. Perhaps an answer was recently deleted and that counts as activity.

Anyway, I am guessing that you are trying to compute the exact determinant of a matrix over the integers. In that case you should use LinearAlgebra:-Determinant(A, method= integer). If your matrix is over some other ring, let me know.

After using the command kernelopts(opaquemodules= false), you can use the command showstat to list the code of any procedure in a module (or package), both exported and local procedures.

So, from your code, I'm guessing that for each integer i from 1 to 100, you want to find the smallest n such that tau(n) = i.  In other words, you want a partial inverse of the tau function. (Note to other readers: tau(n) is the number of positive integer divisors of n.) But your method won't work because for i = 97, n will have 29 decimal digits, and the whole age of the universe wouldn't be enough time to get to that number if you work by counting up from 1. I don't want to say what that 29-digit number is, because it would give too much of a hint.

Hint: You are correct that the prime factorization of n is the key to computing tau(n). Indeed, the primes themselves don't matter; it can be computed just from the exponents in the factorization. Take a look at showstat(numtheory:-tau). It's a very short procedure; only the last actual line of code matters here; the other lines are just there to quickly handle the trivial cases. Once you know what you want the exponents to be, just pick the prime bases that minimize the product. 

By adding sqrt(1+2*Pi), sqrt(1+4*Pi), sqrt(2*Pi), etc., to the bases used by identify (by using the extension option), I've identified several of the numerical solutions.  If cpu time becomes an issue, reduce the degree used by identify to 2 from its default 6. See ?identify

I'd like to take a look at this. What does it say when it crashes? What are the numbers for "Memory" and "Time" in the bottom edge, right side of your Maple window? What values do you get if you execute the commands kernelopts(cpulimit) and kernelopts(datalimit) in a fresh Maple session?

I took a look at your code, and I think that you're wasting cpu time by checking for small and medium factors with igcd, because isprime already does those checks. Take a look at showstat(isprime).

You set N:= (b-a)/h, with floating point a, b, h; so, N is floating point. But when you access N as the upper limit of a seq, it becomes integer. You should say N:= round((b-a)/h),

First 392 393 394 395 Page 394 of 395