Carl Love

## 26862 Reputation

11 years, 304 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## interface(screenwidth)...

The output width of both lprint and showstat is controlled by interface(screenwidth). Its default value is 79, which I recall was a typical actual screenwidth in the 1980s. Modern screens are much wider.

## Horner's is effective in Maple...

There are many stumbling blocks that can occur when trying to time Maple commands, and I think that you've been fooled by one. For example, it's invalid to compare two measurements at least one of which has been rounded down to 0. The same is true for any measurements, even outside of computers.

The following example shows that using Horner's rule is much faster:

```restart:
#degree, number of terms, and number of test-data points:
deg:= 2^6:  nterms:= deg+1:  nData:= 2^12:
p:= x^deg + randpoly(x, degree= deg-1, terms= nterms-1):
(P,Ph):= unapply~([p, convert(p, horner)], x)[]:
#Test data is rational numbers in fraction form:
Data:= `~`[`/`]('rtable(1..nData, random(), subtype= Vector[row])' \$ 2):
gc(); Rh:= CodeTools:-Usage(map(Ph, Data)):
memory used=454.41MiB, alloc change=64.00MiB,
cpu time=359.00ms, real time=702.00ms, gc time=93.75ms

gc(); R:= CodeTools:-Usage(map(P, Data)):
memory used=1.06GiB, alloc change=0 bytes,
cpu time=3.28s, real time=7.22s, gc time=140.62ms

max(abs~(R-Rh)); #Check that results are identical.
0

```

## Updated, simpler Factrix...

There have been many innovations in Maple since Robert Israel wrote Factrix (circa 24 years ago). The following works for any rtable (i.e., vector, matrix, array), of any number of dimensions, whose entries are all exact explicit rational numbers (integers or fractions). It factors out the largest possible constant such that the remaining rtable is all integer.

Factrix:= (M::rtable)->
(C-> `if`(C=0, M, C %* (M/~C)))(((igcd@op@numer~)/(ilcm@op@denom~))({seq}(M)))
:
Factrix(b);
#b is the matrix from the original question.

This is not meant to be a complete replacement for the original Factrix because this only handles rational numbers.

## try ... catch, max[index], densityplot...

The first step is to adjust the functions that you want to plot so that they return an extended_numeric value for all numeric inputs from the plot intervals (alpha= 0..1, delta= 0..1). There's no need for you to understand the subtle distinction between numeric and extended_numeric. You just need to know that undefined is an extended_numeric value that all plot commands can handle gracefully. Currently, diffrP1 and diffrP2 return errors over a large portion of the plot intervals. Most plot commands can handle those errors and can convert them to undefined, but this doesn't work very well when you're trying to indicate the max of the 3 functions.

Here is how to trap the errors and convert them to undefined:

for k to 3 do
DiffrP||k:= subs(
_D= diffrP||k,
proc(a,d) try _D(a,d) catch: undefined end try end proc
)
od:

It is simply defining 3 new functions that call the original 3, "catch" any errors, and return undefined in those cases.

This next command plots those 3 functions in side-by-side plots. This is not strictly necessary for your final plot, but it'll help you understand it.

plots:-display(
plots:-densityplot~(
`<|>`(DiffrP||(1..3)), 0..1, 0..1, title=~ ["diffrP"||(1..3)],
labels= [alpha, delta], labeldirections= ["horizontal", "vertical"],
)
);

I've omitted a direct display of those plots here.

Finally, we create a new function that returns either 12, or 3 to indicate which of the functions returns the maximum value and excludes undefined from consideration. A simple way to write that new function is

max[index, defined]@[DiffrP||(1..3)]

And we plot it like this:

plots:-densityplot(
max[index, defined]@[DiffrP||(1..3)], 0..1, 0..1,
labeldirections= ["horizontal", "vertical"], labels= [alpha, delta],
);

So, the third function is usually the maximum, and there are small regions where the first or second are the maximum.

The colorbar on the right of the plot was a feature introduced in Maple 2023. Without having the bar, you just need to know that red = 3, white = 2, blue = 1.

## applyrule is much older than seq...

[The title refers to the seq parameter modifier, not the seq command.]

Addressing your question about why this inconsistency exists: I guess that a large part of the reason is that applyrule is many years older than parameter modifiers. I'd guess about 10 years.

## Plot x vs. y...

Another method is to plot the inverse function using interchanged axes. Maple's solve and plot commands make that much easier than it sounds.

plot([solve(y = x^(1/3), x), y, y= -2..2], labels= [x,y]);

If you know the inverse, you could do

plot([y^3, y, y= -2..2], labels= [x,y])

Some adjustments to the1st method are needed if the function isn't injective on the desired interval and solve returns multiple inverses.

## Sum...

Like this:

Sum((-1)^r*x^r, r= 0..infinity)

## selectremove...

My code below creates two procedures: LmR converts a set or list of equations and expressions to all expressions; Eq0 converts it to all equations. The key part is selectremove, which separates the input into two subgroups: the equations and the expressions. Then one of those subgroups is converted as desired, and the other is passed through unchanged.

```(LmR,Eq0):= subs~(
_M=~ [(E,A)-> ((lhs-rhs)~(E)[], A[]), (E,A)-> (E[], (A =~ 0)[])],
(S::{set,list}({`=`,thistype}(algebraic)))->
`if`(S::set,`{}`,`[]`)(_M(selectremove(type, S, `=`)))
)[]:
ex:= {a*x^2 + b*x = v, c*x^2 + d*x - w}:
LmR(ex);
Eq0(ex);
/   2               2          \
{ a x  + b x - v, c x  + d x - w }
\                              /

/   2               2              \
{ a x  + b x = v, c x  + d x - w = 0 }
\                                  /

```

## Add an objective variable to linearize...

Note the role of the added variable z below. It is the objective to be minimized, and it's constrained to be greater than or equal to all of the linear expressions from your original objective. Thus it's greater than or equal to their max. By using an extra variable in this way, the objective and constraints are all linear. Having max directly in the objective makes it nonlinear, which is much more complicated to solve. The Optimization package doesn't handle nonlinear integer programs.

Your should be an array of symbolic variables. Do not give it a datatype

 > restart :
 > num_profiles:= 4;  num_websites:= 3;  num_groups:= 2; X:= Matrix((num_profiles, num_websites), symbol= x); objs:= seq(add(X[i,j], j= 1..num_websites), i= 1..num_profiles); constraints:= {     seq(         1 = add(             add(X[i,j], i= k*num_profiles/num_groups + 1 .. (k+1)*num_profiles/num_groups),             k = 0..num_groups - 1         ),         j = 1..num_websites       ),     z >=~ objs }; sol:= Optimization:-Minimize(     z, constraints, integervariables= {z}, binaryvariables= {seq}(X) ); optimal_assignment:= eval(X, sol[2]); for i to num_profiles do     for j to num_websites do         if optimal_assignment[i, j] = 1 then             print("Profile ", i, " assigned to Website ", j)         fi     od od; print("Objective Value (Minimized Maximum Website Overlap): ", sol[1]);

## See help for SurfaceInt...

Surface is not a command per se; rather, it's just a token for parametirically specifying the surface of integration for the SurfaceInt command. The help page ?VectorCalculus,SurfaceInt has the complete instructions for specifyting a surface with Surface in the 4th paragraph of the 2nd bullet point of "Description".

## Calendar:-DayOfWeek...

For example, to do it for year 2024:

select(m-> Calendar:-DayOfWeek(2024, m, 13) = 6, [\$1..12]);
[9, 12]

So the 9th and 12th months have a Friday-the-13th, i.e., September and December. The 6 represents Friday (considered the 6th day of the week).

## Type suffixed...

Use

if indets(r, suffixed('_')) <> {} then . .

You could also change '_' to _Z, etc.

## Left-associative, and mod not distribute...

There are two anomalies happening in your example:

1. Left-associativity: All infix operators beginning with a single are left-associative (and those beginning with && are right-associative). So 2 &^ 2 &^ 2 &^ 2 is interpreted as ((2 &^ 2) &^ 2) &^ 2 rather than 2 &^ (2 &^ (2 &^ 2)). This happens in the kernel before mod sees the expression.

2.Distributivity: mod doesn't distribute over exponentiation because it's usually not mathematically valid. In other words, a &^ b mod c isn't interpreted as (a mod c) ^ (b mod c) mod c. An example for which it isn't valid is a=2, b=4, c=3

## Just simplify...

I wouldn't use testeq at all. It's returned FAIL every time that I've ever tried to apply it to a transcendental function. But plain simplify or is works fine here. All of these return the information that you want:

is(expr=x), evalb(simplify(expr)=x), evalb(simplify(expr-x) = 0);

## Two better algorithms...

Here are two more-efficient algorithms for the sum of the proper divisors. The first of these (SPD2 below) is much faster than your original; the second (SPD3 below) is much faster than that.

I think that this will all work in your Maple 13; if not, let me know.

```restart
:
(* Sum of Proper Divisors (original version):
---------------------------------------------
This is your original procedure. The only changes I made were removing things
not directly related to getting the answer, such as print statements.
*)
SPD0:= proc(n)
local Sum1, b;
Sum1:= 1;
for b from 2 to ceil(n/2) do
if floor(n/b) = n/b then Sum1:= Sum1 + b fi
od;
Sum1
end proc
:
(* Sum of Proper Divisors (version 1):
--------------------------------------
This is the same algorithm as yours, recoded. The changes that I made were
to use operators that are more efficient for integer arithmetic: add, modp,
and iquo instead of `for`, ceil, `/`, and floor.
*)
SPD1:= proc(n::posint)
local d;
1 + add(`if`(modp(n,d) = 0, d, 0), d= 2..iquo(n,2))
end proc
:
(* Sum of Proper Divisors (version 2):
--------------------------------------
This algorithm uses this idea:
Let n be a positive integer. Let s = floor(sqrt(n)). For every divisor d of n,
q = n/d is also a divisor, and d <= s iff q >= s. If s^2 = n, we need to adjust
for the double counting of s.
*)
SPD2:= proc(n::posint)
local d, q, s;
s:= isqrt(n);  #same as floor(sqrt(n)) but faster
#irem gives the integer remainder and quotient in one step:
add(`if`(irem(n, d, 'q') = 0, d+q, 0), d= 2..s)  +  1 - `if`(s^2=n, s, 0)
end proc
:
(* Sum of Divisors of a Prime Power:
------------------------------------
For prime p, the divisors of p^e are clearly
{1, p, p^2, ..., p^e}.
The sum of those is a simple geometric sum:
sum(p^k, k= 0..e) = (p^(e+1) - 1)/(p - 1).
*)
SDPP:= (p::prime, e::posint)-> (p^(e+1) - 1)/(p-1)
:
(* Sum of Proper Divisors (version 3):
--------------------------------------
This uses the classic formula for the sum of divisors, then subtracts
the number itself.

This uses ifactors, which returns the same information as ifactor, but in a
list of lists format:
ifactors(n) = [-1 or +1, [[p1,e1], [p2,e2], ..., [pk,ek]]],
such that the prime factorization is
n = p1^e1 * p2^e2 * ... * pk^ek.
*)
SPD3:= proc(n::posint)
local P;
mul(SDPP(P[]), P= ifactors(n)[2]) - n
end proc
:
#Test all 4 on your example:
SPD||(0..3)(945);
975, 975, 975, 975

#Construct a much larger random example:
seq('rand(2..99)(), nextprime(rand(10^5..10^7)())', k= 1..9);
94, 8633311, 97, 9909901, 60, 1584343, 70, 3888067, 97, 3079871,
91, 5954371, 19, 8014679, 57, 682901, 84, 4397777

n:= `*`(%);
n := 7153933542604628643414930949912947237856549377795642829298491576067369227200

CodeTools:-Usage( SPD3(n) );
memory used=9.62MiB, alloc change=0 bytes,
cpu time=94.00ms, real time=160.00ms, gc time=0ns

28515792740732679749844128611744262891259901154773757373417700467974788852800

```

 5 6 7 8 9 10 11 Last Page 7 of 386
﻿