Carl Love

Carl Love

27291 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

@michaelvio 

The vast majority of the slowness that you were experiencing was due to your use of int (symbolic integration) versus Int (numeric integration, when it's used in the context that you were using).

It's true that using a larger epsilon and a smaller numpoints will save some time, but that savings is infinitesimal in this case. I recommend that you remove both options so that those things will be done at their default values (although epsilon= 1e-6 is sufficient for plotting purposes).

You should also set Digits:= 15 right after the restart at the top of the worksheet. The extra precision makes a huge difference in the final plot, and the extra time needed isn't even noticeable.

To construct the function J1 as you just described, all that you need to do is

J1:= unapply(subs(x0= 1/t0, J), t0);

This is your last plot at the default Digits:= 10:

And this is using Digits:= 15:

It's a similar issue as with your recent Question regarding applying eval to Int or Intateval uses rules, usually to ensure the mathematical validity of the proposed substitution. The rule in this case is the procedure `eval/limit`. To create a more-sophisticated rule, read the procedure (showstat(`eval/limit`), it's a fairly easy read with only a moderate amount of Maple knowledge, and it's only 42 lines (Maple 2023)) and either write an overload or, perhaps, modify the procedure directly (ToInertFromInert).

If you'd like it to follow rules additional to what's already in the procedure, and you can state them for me with a moderate amount of mathematical precision and generality, I may be able to suggest a way to code it. 

Let's consider your examples individually:

1. e:= A+B*limit(f(x), x= infinity);  eval(e, limit(f(x), x= infinity)= 1);
This case doesn't involve changing anything inside any limit. Thus, the rule isn't used (verified by using trace(`eval/limit`)), and you got your expected result.

2. e:= A+B*limit(2*f(x), x= infinity);  eval(e, limit(f(x), x= infinity)= 1);
In this case you're asking it to evaluate an expression containing a limit by using information about another limit (albeit, a very closely related one); the eval rule will be used. I think that the mathematical rule (i.e., abstract, not Maple) that you're implicitly invoking here is always valid, although I'm not 100% sure about all weird "corner cases", domains other than 2-sided real, and multivariate (more than one bound variable (e.g., limit(..., (x,y)= ...)) cases. That rule, stated generally, is that factors that don't depend on the bound variable(s) can be factored out of the limit. Command expand will do this, so a quick solution is 
eval(expand(e), limit(f(x), x= infinity)= 1)

[to be continued]

The correct command is is

is(sqrt(3) < 3);

Presumably it just expands things to a sufficient number of decimal places to make sure.

You wrote:

  • Should one then use is(n>=m) instead of evalb(n>=m) always?

No! That's inefficient. evalb is builtin (kernel code rather than Maple). is is a highly symbolic command written in Maple. In the vast majority of cases, you know a priori that both sides are type numeric, and is isn't needed.

  • I can also always apply evalf() on each side before. But this seems like an overkill to me.

I think that's the approach that is takes, and it might increase Digits if the user's set value is not high enough to make absolutely sure of the inequality.

I said that I could probably come up with something better for large n and < < n, and here it is.

To make the algorithm fast, there are three things to avoid:

  1. finding all the large prime factors of n when that is difficult to do,
  2. finding all the divisors and then filtering them when the number of divisors is large,
  3. checking every positive integer < m when m is large.

The command ifactors(n, easy) finds only the prime factors that are easy to find.

Small Divisors of Integers

A fast algorithm for finding the divisors less than some upper bound of a large integer

 

Text and code author: Carl Love <carl.j.love@gmail.com> 2024-Sept-27

 

restart
:

interface(prompt= "")
:

To understand my SmallDivisors algorithm, you need to understand the output of the integer factorization command ifactors(..., easy). So, for example, I'll build an integer n that has an easy-to-factor part and a hard-to-factor part and apply the 'easy' algorithm to it. I believe (although I'm not absolutely sure) that the 'easy' algorithm will find any prime factor < 2^22. It will often also find nontrivial prime factors larger than that.

p1:= prevprime(2^25); 2^25-p1;

33554393

39

p2:= nextprime(p1);

33554467

p3:= nextprime(p2);

33554473

n:= (2*3*5*7)^32*(p1*p2*p3)^2;

292090706291643821582093796327445337477364073600065563995291265004658240987511107950432900000000000000000000000000000000

The hard-to-factor part will be p2*p3.

p2*p3;

1125902456980891

1+ilog10(%); #Count its base-10 digits.

16

ifactors(n, easy);

[1, [[2, 32], [3, 32], [5, 32], [7, 32], [33554393, 2], [_c16_1, 2]]]

So, the above says that n = 1 * 2^32* 3^32 * 5^32 * 7^32 * (2^25 - 39)^2 * "a 16-digit composite"^2.

 

Here's the procedure:

(*
Fast algorithm to return the divisors of n less than m.

This procedure may return a syntax error if entered as 2D Input.
Please use 1D input (aka, Maple Input or "plaintext").
*)
SmallDivisors:= proc(n::And(posint, Not(1)), m::And(posint, Not(1)))
description "divisors of n less than m";
option `Author: Carl Love <carl.j.love@gmail.com> 2024-Sept-28`;
local
    F:= ifactors(n, 'easy')[2],
    PF:= remove(type, F, [name,posint]),
    P:= op~(1,PF), #the prime factors
    #Get max possible exponent for each prime:
    E:= min~(op~(2,PF), trunc~(evalhf~(ln(m)/~(ln~(P))))),
    r:= {1}, p, k, D, d
;
    if ilog2(m) > 21 and F[-1][1]::name and E[-1] <> 0 then
        WARNING("Some divisors < %1 may be missing", m)
    fi;
    for k,p in P do
        D:= r[];
        r union= {to E[k] do D:= for d in D while d*p < m do d*p od od}    
    until E[k]=0
end proc
:

Use it on Kitonum's example:

d1:= CodeTools:-Usage(SmallDivisors(10^100+1, 10^6)): d1;

memory used=301.66KiB, alloc change=0 bytes, cpu time=0ns, real time=101.00ms, gc time=0ns

{1, 73, 137, 401, 1201, 1601, 10001, 29273, 54937, 87673, 116873, 164537, 219337, 481601, 642001}

Compare with Kitonum's method, which obviouly has time complexity O(m):

Brute:= (n,m)-> local k; {seq(`if`(irem(n,k)=0, k, NULL), k= 1..m-1)}:

d2:= CodeTools:-Usage(Brute(10^100+1, 10^6)): d2;

memory used=146.08MiB, alloc change=52.63MiB, cpu time=47.00ms, real time=618.00ms, gc time=0ns

{1, 73, 137, 401, 1201, 1601, 10001, 29273, 54937, 87673, 116873, 164537, 219337, 481601, 642001}

Compare using a larger upper bound m:

d1:= CodeTools:-Usage(SmallDivisors(10^99+1, 10^7)): d1;

memory used=214.85KiB, alloc change=0 bytes, cpu time=0ns, real time=6.00ms, gc time=0ns

{1, 7, 11, 13, 19, 23, 77, 91, 121, 133, 143, 161, 209, 247, 253, 299, 437, 847, 1001, 1463, 1573, 1729, 1771, 2093, 2299, 2717, 2783, 3059, 3289, 4093, 4807, 5681, 8779, 11011, 16093, 19019, 19481, 23023, 28651, 29887, 33649, 36179, 39767, 45023, 52579, 52877, 53209, 61453, 62491, 77767, 94139, 96569, 114127, 166801, 201917, 209209, 253253, 315161, 368053, 370139, 372463, 437437, 495253, 544369, 578369, 585299, 658973, 675983, 683527, 687401, 798889, 855437, 999001, 1010971, 1035529, 1062259, 1167607, 1209317, 1223807, 1255397, 1413419, 1788641, 1834811, 2168413, 2221087, 2624921, 3466771, 3836423, 4048583, 4097093, 4784689, 4811807, 5988059, 6362059, 6438289, 6993007, 7076797, 7248703, 7435813, 7518797, 8465219, 8566649, 8787779, 9409807}

d2:= CodeTools:-Usage(Brute(10^99+1, 10^7)): d2;

memory used=1.22GiB, alloc change=96.66MiB, cpu time=1.98s, real time=7.04s, gc time=1.53s

{1, 7, 11, 13, 19, 23, 77, 91, 121, 133, 143, 161, 209, 247, 253, 299, 437, 847, 1001, 1463, 1573, 1729, 1771, 2093, 2299, 2717, 2783, 3059, 3289, 4093, 4807, 5681, 8779, 11011, 16093, 19019, 19481, 23023, 28651, 29887, 33649, 36179, 39767, 45023, 52579, 52877, 53209, 61453, 62491, 77767, 94139, 96569, 114127, 166801, 201917, 209209, 253253, 315161, 368053, 370139, 372463, 437437, 495253, 544369, 578369, 585299, 658973, 675983, 683527, 687401, 798889, 855437, 999001, 1010971, 1035529, 1062259, 1167607, 1209317, 1223807, 1255397, 1413419, 1788641, 1834811, 2168413, 2221087, 2624921, 3466771, 3836423, 4048583, 4097093, 4784689, 4811807, 5988059, 6362059, 6438289, 6993007, 7076797, 7248703, 7435813, 7518797, 8465219, 8566649, 8787779, 9409807}

evalb(d1 = d2);

true

Now let's consider a number with a great many small divisors. The n1 below has 3^14 (~4.8 million) divisors.

n1:= mul(ithprime(k), k= 1..14)^2;

171158644061440576710668800200900

d3:= CodeTools:-Usage(SmallDivisors(n1, 10^9)): nops(d3);

memory used=10.56MiB, alloc change=4.27MiB, cpu time=0ns, real time=79.00ms, gc time=0ns

133547

d4:= CodeTools:-Usage(select(`<`, NumberTheory:-Divisors(n1), 10^9)):

memory used=0.86GiB, alloc change=206.21MiB, cpu time=3.44s, real time=6.41s, gc time=2.89s

evalb(d3=d4);

true

SmallDivisors has comparable performance with NumberTheory:-Divisors when all divisors are requested:

d5:= CodeTools:-Usage(SmallDivisors(n1, n1+1)): nops(d5)=3^14;

memory used=443.29MiB, alloc change=118.64MiB, cpu time=2.00s, real time=4.19s, gc time=1.70s

4782969 = 4782969

d6:= CodeTools:-Usage(NumberTheory:-Divisors(n1)): nops(d6);

memory used=322.98MiB, alloc change=-2.08MiB, cpu time=3.33s, real time=3.77s, gc time=3.16s

4782969

The following number has 15! (~1.3 trillion) divisors, so it would be infeasible to apply NumberTheory:-Divisors to it due to memory limitation. Howver, SmallDivisors can easily find all divisors less than 1 trillion.

n2:= mul(ithprime(i)^(15-i), i= 1..14);

222046261244808869776040900239313953120666458401987626347805929809038681543479663134628000000000000

d7:= CodeTools:-Usage(SmallDivisors(n2, 10^12)):
nops(d7), d7[-1], ifactor(d7[-1]);

memory used=456.89MiB, alloc change=27.00MiB, cpu time=2.03s, real time=3.45s, gc time=1.86s

6055806, 999999663336, ``(2)^3*``(3)^5*``(17)*``(23)*``(31)^2*``(37)^2

 

 

 

 

Download SmallDivisors.mw

This may be suitable for you:

select(`<`, Divisors(n), m)

Surely it's more efficient than what you were doing. For an easy-to-factor n >> 2^64 and m << sqrt(n), I could probably come up with something even better.

My guess is that you actually want the simplified product of an arbitrary number of factors of T[r]*P[r] --

T[1]*P[1]*...*T[n]*P[n]

-- because otherwise your Question would be exceedingly trivial. But first, you should stop trying to use an index (your r) as a parameter. Maple does allow the use of indices as parameters, but it's a somewhat more advanced topic than your current level of Maple knowledge. So, for now, let's just make r an ordinary parameter:

T:= r-> exp(r+1);  P:= r-> exp(r-1);

Your product can be represented as an inert product with the Product command:

Product(T(r)*P(r), r= 1..3);

             

Now, here's the non-trivial (but still conceptually easy) part: Make the upper bound, 3, arbitrary, and make the product non-inert by using product (lowercase):

product(T(r)*P(r), r= 1..n);

             

Of course that can be greatly simplified:

simplify(%);

             


By the way, can anybody give me an English word that's actually used in mathematical prose that completes this formal two-level analogy?:

(Addition is to summation as multiplication is to __?__)
as in Maple
(add is to sum as mul is to product).

That word, if it exists, would replace "products" in this Answer's Title.

Their are 9 vertices that I'll order top-down, left-to-right. That makes 18 variables. By making the lower right vertex the origin and the right edge the y-axis, we can make 5 variables 0. By reflecting the given image over the y-axis, we can constrain all points to the 1st quadrant. I generate 13 equations from the given information in your diagram, then use fsolve. The 13 equations are:

  • Perimeters of each of the 7 primary triangles
  • Perimeters of the 2 composite triangles
  • Collinearity of the left side (2 equations because it's 4 points)
  • Collinearity of the bottom side
  • Collinearity of the bottoms of the triangles with perimeters 9 and 12.

I'll leave it to you to color them.

restart
:

#9 points, each with x and y coordinates.
P:= Matrix((9,2), symbol= p)
:

#distance-between-points function:
d:= (i,j)-> sqrt((p[i,1]-p[j,1])^2 + (p[i,2]-p[j,2])^2)
:

#perimeter-of-triangle:
T:= (i,j,k)-> d(i,j)+d(j,k)+d(k,i)
:
#co-linearity constraint generator
CoL:= (i,j,k)-> d(i,j)+d(j,k) = d(i,k)
:
#the primary triangles:

Ts:= [[1,2,3], [2,4,5], [2,3,5], [3,5,6], [4,7,8], [4,6,8], [6,8,9]]
:
#their perimeters:
Ps:= [20, 9, 13, 12, 10, 15, 11]
:

#Make lower right corner the origin, and right edge the y-axis:
eqs0:= {p[9,1], p[9,2], p[1,1], p[3,1], p[6,1]}=~ 0
:
eqs:= eval(
    {
        ((T@op)~(Ts) =~ Ps)[],
        CoL(1,2,4), CoL(2,4,7), CoL(4,5,6), CoL(7,8,9),
        T(1,4,6) = 20+9+12 - 13, T(1,7,9) = 20+9+12+10+11 -13-15
    },
    eqs0
):

#1st quadrant:
Q1:= ({seq}(P) minus lhs~(eqs0))=~ 0..infinity:

P:= eval(P, fsolve(eqs, Q1) union eqs0);

Matrix(9, 2, {(1, 1) = 0, (1, 2) = 14.90060263, (2, 1) = 3.382184235, (2, 2) = 10.23285104, (3, 1) = 0, (3, 2) = 6.046658910, (4, 1) = 5.188020843, (4, 2) = 7.740618787, (5, 1) = 3.710858752, (5, 2) = 6.351898354, (6, 1) = 0, (6, 2) = 2.863276643, (7, 1) = 5.954748645, (7, 2) = 6.682472994, (8, 1) = 3.216507499, (8, 2) = 3.609582404, (9, 1) = 0, (9, 2) = 0})

Poly:= T->
    plottools:-polygon(
        convert(P[T], listlist), linestyle= solid, color= white
    )
:
Label:= k-> local i;
    plots:-textplot([seq(add(P[Ts[k][i]], i= 1..3)/3), Ps[k]])
:

plottools:-rotate(
    plottools:-reflect(
        plots:-display(
            Poly~(Ts)[], seq(Label(k), k= 1..nops(Ts)),
            axes= none, scaling= constrained, gridlines= false
        ),
        [[0,0],[0,1]]  #reflect over y-axis
    ),
    Pi/6
);

;

 

Download TrianglesByPerimeters.mw

Your initial expression V__phi has a pair of square brackets [ ]. These need to be replaced by round brackets ( ) (also known as parentheses). In Maple (and most computer languages) square brackets can't be used as a "higher level" of parentheses like they are used in math books. Only parentheses can be used for algebraic groupings.

Making this change, the whole expression simplifies quickly.

You can use Edit ==> Find from the menus to find the square brackets.

The read command is for reading files of plaintext Maple code, not Maple worksheets. The file that you tried to read, "EKHAD.mw", is a worksheet. You should open it like you'd open any other worksheet: through the File ==> Open menu dialog.

The first read command that you gave has already loaded the code. Just type

ezra();

I believe that many of the algorithms from A=B have been directly implemented in Maple already.

The difference between eval (in its form used in this Question) and subs is that eval follows rules that are often specified in Maple code. The rule in this case is the procedure `eval/Int`. Presumably this procedure has been updated since Maple 2015. After a restart (to clear remember tables), do

trace(`eval/Int`);

and rerun your problem.

This should work for all your cases:

SeriesInt:= e->
    subsindets(
        e, specfunc({Int, int, Intat, intat}),
        J-> op(0, J)(mtaylor(op(1, J), op([2,1], J)= 0, 5), op(2.., J))
    )
;

The suggestion by @dharr will work in most cases. But it won't work when the subscript is an integer because then the subscript will be interpreted as an index into the list on its left. The command Typesetting:-Typeset can be used to turn mathematical expressions into inert objects that display exactly like the original expressions. So, this works:

restart:
`print/SubSup`:= (e, a, b)->
    [e][Typesetting:-Typeset(a)]^Typesetting:-Typeset(b)
:
SubSup:= (e::uneval, a, b)-> 'procname'(e,a,b)
:
SubSup(int(x^2, x), 2, 5);


 

The plot command is only for 2D. For 3D, use

plots:-spacecurve([t*sin(t), t*cos(t), (2*sqrt(2))/3*t^(3/2)], t= 0..Pi);

For 2D (assuming (just to keep this Answer simple) the orthogonal projection into the xy-plane), use

plot([t*sin(t), t*cos(t), t= 0..Pi]);

Note that in the first command, the range t= 0..Pi is outside the square brackets, and in the second command, it's inside.

(simplify@eval)(
    sqrt(a^2+b^2),
    solve({d-b = c-d, a/d = tan(45*Pi/180), a*(c-b)/2 = 60, a^2+c^2 = 19^2})
);

                               11

The variables sx, and y are clearly irrelevant because they only appear once in the problem.

Since the constraints are all easily made polynomial, this works also:

(simplify@@2)(
    sqrt(a^2+b^2), 
    {d-b = c-d, a/d = tan(45*Pi/180), a*(c-b)/2 = 60, a^2+c^2 = 19^2}
);

                               
11

This form of simplify requires only the constraints, not the objective, to be polynomial, or easily convertible to polynomial (as in a/d = ...). It's help page is ?simplify,siderels.

The extra simplify (as in simplify@@2) is unfortunately required to simplify sqrt(121) to 11. To avoid that awkwardness, you could do this:

(sqrt@simplify)(
    a^2+b^2, 
    {d-b = c-d, a/d = tan(45*Pi/180), a*(c-b)/2 = 60, a^2+c^2 = 19^2}
);

                               11

The other responders are making the mistake of assuming that all the variables are integer. The problem only requires AD to be integer. On a number line, let A, B, C, D 0, 7/3, 20/3, 21/3. Then AD = 7.

This is an excellent problem for testing any mathematical natural-language AI system. Solving it requires more than just middle-school Algebra I. This is a GRE-level problem (admittedly a fairly easy one). (GRE = Graduate Record Exam, a nationally standardized exam that one takes in the U.S. before applying to graduate school.)

@mmcdara In English, there is no ambiguity or lack of precision in this problem's statement. "Consecutive" means that they are distinct points, and are in the stated order. "Confounded" is not a word used in ordinary mathematical English; its only common usages are to describe mental states.

The diagram has two quadrilaterals that I'll call the "kite" and the "dart". They share two sides. I'll label these points as vertices:

  • A: point at the 60-degree angle
  • N: point where n is used
  • C: point in center where two lines cross
  • M: point where is used
  • X: point where x is used.

The kite is ANCM; the dart is XNCM. Now make two equations using that the sum of the angles of a quadrilateral is 360 degrees.

#quadrilateral ANCM:
Kite:= 60 + (180-2*n) + (180-40) + (180-2*m) = 360;
#       A       N           C          M
                 Kite := 560 - 2 m - 2 n = 360

#quadrilateral XNCM:
Dart:= x + n + (180+40) + m = 360;
#      X   N       C      M
                 Dart := x + m + 220 + n = 360

x = eval(x, solve({Kite, Dart}));
                             x = 40

 

1 2 3 4 5 6 7 Last Page 1 of 390