The original poster mentioned that the vector should be "not of a known length or known values." Your example does not match that abstract description.

acer

The original poster mentioned that the vector should be "not of a known length or known values." Your example does not match that abstract description.

acer

It occurred to me that someone might want such a fast hardware J0 to be utilized even in cases where the original `BesselJ` function was hard-bound into some other Maple Library routine.

By all that I mean the case where some other Library routine already had reference to BesselJ in it when it was read from source and saved into the Maple Library. In such a case it is the global name BesselJ which has been bound, prior to that other routine being saved. So it wouldn't change the behaviour if, say, I loaded a module which exported its own BesselJ. No, that would only affect the instances where I typed in BesselJ in my session, with that export's binding in effect. It wouldn't change the behaviour for those routines which were saved with reference to the global name in them already. In order to affect those routines I would instead have to redefine the global BesselJ name.

So below I show an example of that. I would advise that you don't run the code unless you understand what it does. I wouldn't want you inadvertantly to wreck the BesselJ in your installed Maple Library.

What the code below does is unprotect BesselJ, and in its place save a version that calls the fast external library when n=0 and x::numeric, and in all other cases call the original BesselJ. The original BesselJ can be referenced due to having copied it right at the start.

Since external calls may be session dependent, it should not be possible to properly save the call_external. So instead the Maple-level entry routine to the fast external function will have to redefine itself whenever it gets called the first time in a new session. This is brought about with a ModuleLoad action. See here for some comments on that.

> hwBJ := module()
> option package;
> export hwBesselJ0;
> local ModuleLoad;
>   ModuleLoad := proc()
>     # The OS libm.so may have double-precision Bessel J0.
>     unprotect(hwBesselJ0);
>     hwBesselJ0 := define_external('j0',
>        'r'::float[8],'RETURN'::float[8] ,
>        LIB="/usr/lib64/libm.so"):
>     protect(hwBesselJ0);
>   end proc;
> end module:
>
> __origBesselJ := copy(eval(BesselJ)):
> __origBesselJ(0,5.6);
                                 0.02697088469
>
> unprotect(BesselJ):
> BesselJ := proc(n,x)
>   if n=0 and Digits<=trunc(evalhf(Digits))
>    and type(x,'numeric') then
>     hwBJ:-hwBesselJ0(x);
>   else
>     __origBesselJ(n,x);
>   end if;
> end proc:
>
> march(create,"hwBJ.mla");
> libname:="hwBJ.mla",libname;
          libname := "hwBJ.mla", "/home/acer/maple/maple11.02/lib"

Now I can save all three parts: the module which sets up the call_external and then redefines its own export to use it, the copy of the original BesselJ, and my new modified BesselJ.

> savelib(hwBJ);
> savelib(__origBesselJ);
> savelib(BesselJ);

Now I can restart, and test it. I want it to use the hardware J0 when Digits is not high, and when the input is real and numeric, but otherwise I want the usual BesselJ behaviour. If I were being really careful I would have tested the hardware J0 to find out for what range of values it is accurate, and then restricted its use to that range

> restart:
> libname:="hwBJ.mla",libname;
          libname := "hwBJ.mla", "/home/acer/maple/maple11.02/lib"
 
>
> BesselJ(0,100.0);
                             0.0199858503042231184
 
> evalf(BesselJ(0,17));
                                 -0.1698542522
 
> BesselJ(0,y);
                                 BesselJ(0, y)
 
> BesselJ(1,100.0);
                                -0.07714535201
 
> Digits:=30: BesselJ(0,100.0); Digits:=10:
                       0.0199858503042231224242283909508
 
>
> evalhf(eval(BesselJ(0,17)));
                             -0.169854252151183549

At this point I should mention that I don't offhand know of places where BesselJ is actually hardcoded into Library sources. But these methods should be feasible for other special functions as well. And I do know that, for example, GAMMA is coded into all sorts of routines including some for floating-point calculations.

It might be fun to work out robust Maple platform-independent code to reliably replace a useful key set of special functions (when used at hardware or lower precision, and over ranges where they are accurate).

acer

It's not really like how the GMP (GnuMP) library is used, as that is quite tightly tied into the Maple kermel which is linked directly against it and since its data structures don't appear directly at the Maple Library level.

It is much closer to how the NAG library is used from the Maple Library, where a external library specialized in numerics is dynamically opened. The data may be visible at the Library level (hardware arrays and, now, even hardware scalars) and some Library routine decides (algorithmic switching, often invisible to the user) on usage according to problem type.

acer

I'm sorry that I don't understand the physical problem well enough to comment on your alternative method.

You asked how to get the operator body. You had this,

power:=theta -> 4*p[0]*BesselJ(1, k*r*sin(theta))^2/(k^2*r^2*sin(theta)^2);

Since the only place that theta is used in `power` is in the calls to `sin` (which returned unevaluated when passed a name) you could just call power(theta).

> power(theta);
                                                       2
                      4 p[0] BesselJ(1, k r sin(theta))
                      ----------------------------------
                               2  2           2
                              k  r  sin(theta)

Trying to be fancy..

> op(map(FromInert, subs({_Inert_PARAM(1)=_Inert_NAME("theta")},
> indets( ToInert(eval(power)), specfunc(anything,_Inert_STATSEQ) ))));
                                                       2
                      4 p[0] BesselJ(1, k r sin(theta))
                      ----------------------------------
                               2  2           2
                              k  r  sin(theta)

acer

I'm sorry that I don't understand the physical problem well enough to comment on your alternative method.

You asked how to get the operator body. You had this,

power:=theta -> 4*p[0]*BesselJ(1, k*r*sin(theta))^2/(k^2*r^2*sin(theta)^2);

Since the only place that theta is used in `power` is in the calls to `sin` (which returned unevaluated when passed a name) you could just call power(theta).

> power(theta);
                                                       2
                      4 p[0] BesselJ(1, k r sin(theta))
                      ----------------------------------
                               2  2           2
                              k  r  sin(theta)

Trying to be fancy..

> op(map(FromInert, subs({_Inert_PARAM(1)=_Inert_NAME("theta")},
> indets( ToInert(eval(power)), specfunc(anything,_Inert_STATSEQ) ))));
                                                       2
                      4 p[0] BesselJ(1, k r sin(theta))
                      ----------------------------------
                               2  2           2
                              k  r  sin(theta)

acer

At the author's special request due to licensing concerns, this comment has been replaced by

Will

At the author's special request due to licensing concerns, this comment has been replaced by

Will

Sure. What surprised me here was that both input objects started off a the table reference p[0], and to both I did the same context-menu action (or Format menu), and the result was two atomic identifiers which Maple did not see as the same name.

acer

Sure. What surprised me here was that both input objects started off a the table reference p[0], and to both I did the same context-menu action (or Format menu), and the result was two atomic identifiers which Maple did not see as the same name.

acer

The Maple input in your material is in a format that is called 2D Math. That's the format in which the input can appear as nicely marked up math with symbols and all that jazz.

In 2D Math, a name with a subscript (which in your example is entered as p_0) is actually the table reference p[0] under the hood.

If you select a section of input then you can use the menus (top-menu's Format -> Convert To, or right-click context-menu's Convert To) to toggle the selected between 2D Math and 1D input.

So, when you had p[0] assigned a value, then Maple got all confused when you tried to assign a procedure to p (which, incidentally, made use of p[0]) and then to plot it. Very often in Maple, it is not possible to use a name like p for one thing while also using its indexed entries for something else.

There's a way around this, if one really wants to use 2D Math. After typing in your p_0 (which is p[0] at that moment) you can select it with the mouse an use the right-click context menu and convert that subscripted item into its own unique name. That unique name, which in 2D Math would still appear as the subscripted p, would be fully distinct from p and so not suffer from a collision elsewhere. The context-menu choice would be 2D Math -> Convert To -> Atomic Identifier.

The idea is that, by toggling a subscripted name (identifier) as atomic (unique) it would be distinct from the unsubscripted same name.

There are some snags. The first, which is just bearable, is that you have to toggle every single instance of inputing p_0 in 2D Math. You could also use the `subliteral` palette object and get the same effect, but that would mean having to keep switching between keyboard and mouse action when entering math. The second is that table references for subscripted names is probably the wrong default, and can't be turned off as a preference. And the term "atomic identifier" is unintuitive and so not as good as "unique name". And the whole business is rather underdocumented, though it sure does come up a lot here on mapleprimes.

When I tried toggling the instances of p_0 in your uploaded worksheet (which were p[0]) I found that they produced different unique names. So the procedure which contained the p_0 would not be able to evaluate to an actual real number, since its p_0 was not assigned! By using the menus to convert both instances to 1D input, after I had toggled them as Atomic Identifiers, I found that they were different names underneath. The first one, where the assignment took place, was,

`#msub(mi("p",fontstyle = "normal"),mn("0"))`

The second one, inside the operator definition, was,

`#msub(mi("p"),mn("0"))`

And the assignment to the former would not get utilized when accessing the latter. For me, that was a new twist on this issue.

acer

The Maple input in your material is in a format that is called 2D Math. That's the format in which the input can appear as nicely marked up math with symbols and all that jazz.

In 2D Math, a name with a subscript (which in your example is entered as p_0) is actually the table reference p[0] under the hood.

If you select a section of input then you can use the menus (top-menu's Format -> Convert To, or right-click context-menu's Convert To) to toggle the selected between 2D Math and 1D input.

So, when you had p[0] assigned a value, then Maple got all confused when you tried to assign a procedure to p (which, incidentally, made use of p[0]) and then to plot it. Very often in Maple, it is not possible to use a name like p for one thing while also using its indexed entries for something else.

There's a way around this, if one really wants to use 2D Math. After typing in your p_0 (which is p[0] at that moment) you can select it with the mouse an use the right-click context menu and convert that subscripted item into its own unique name. That unique name, which in 2D Math would still appear as the subscripted p, would be fully distinct from p and so not suffer from a collision elsewhere. The context-menu choice would be 2D Math -> Convert To -> Atomic Identifier.

The idea is that, by toggling a subscripted name (identifier) as atomic (unique) it would be distinct from the unsubscripted same name.

There are some snags. The first, which is just bearable, is that you have to toggle every single instance of inputing p_0 in 2D Math. You could also use the `subliteral` palette object and get the same effect, but that would mean having to keep switching between keyboard and mouse action when entering math. The second is that table references for subscripted names is probably the wrong default, and can't be turned off as a preference. And the term "atomic identifier" is unintuitive and so not as good as "unique name". And the whole business is rather underdocumented, though it sure does come up a lot here on mapleprimes.

When I tried toggling the instances of p_0 in your uploaded worksheet (which were p[0]) I found that they produced different unique names. So the procedure which contained the p_0 would not be able to evaluate to an actual real number, since its p_0 was not assigned! By using the menus to convert both instances to 1D input, after I had toggled them as Atomic Identifiers, I found that they were different names underneath. The first one, where the assignment took place, was,

`#msub(mi("p",fontstyle = "normal"),mn("0"))`

The second one, inside the operator definition, was,

`#msub(mi("p"),mn("0"))`

And the assignment to the former would not get utilized when accessing the latter. For me, that was a new twist on this issue.

acer

If Maple always went to the expense of trying to do such simplifications then it would hurt the performance for anyone who didn't need or want it for a particular example.

How would one go about turning that simplification off?

And what if there were radicals in the coefficients of the cubic? It could become very expensive to try to determine whether the roots were purely real.

Would drawing the line at, say, rational coefficients appear inconsistent?

And what if one didn't want a result with trig in it, but wanted the explicit radicals only (when possible) from `solve`?

It may be difficult to justify an automatic and potentially expensive simplification which might not always suit everyone.

acer

If Maple always went to the expense of trying to do such simplifications then it would hurt the performance for anyone who didn't need or want it for a particular example.

How would one go about turning that simplification off?

And what if there were radicals in the coefficients of the cubic? It could become very expensive to try to determine whether the roots were purely real.

Would drawing the line at, say, rational coefficients appear inconsistent?

And what if one didn't want a result with trig in it, but wanted the explicit radicals only (when possible) from `solve`?

It may be difficult to justify an automatic and potentially expensive simplification which might not always suit everyone.

acer

I believe that you can insert a URL by selecting some existing text with the mouse and then, with that text highlighted, using the Insert/Edit Link button.

If you haven't yet written that part of the text then you can insert the link first and then edit it afterwards as source. By that I mean switch to source-mode using the "Source" button at the top of the Editor. Then edit by hand what lies between the <a..> and </a> (leaving the href part alone).

It's mildly inconvenient if the text portion cannot be entered or edited in the Insert/Edit Link pop-up.

That numeric integral is hard, I think, to get with good accuracy. Which is something that the author did not address. How accurate an answer did he want? How accurate an answer did he receive from Mathematica?  He appeared to have ignored Mathematica's printed advice and bumped up the MaxRecursion parameter but not chosen Method->Oscillatory. He did not mention that Maple's evalf(Int(..) has [documented] separate options to control both the working precision and an accuracy tolerance.

What I liked about the example is that it illustrates a dichotomy for software. Should it always just try and be super smart and magically handle everything? Or should it provide lots of powerful options and controls to assist with the nastier problems? How should it balance ease of use versus power? And how should it guide the user? Maple's evalf(Int()) can dump out a lot of userinfo messages, if infolevel[`evalf/int`] is set higher than 0. But I don't think that it will use Maple's WARNING() mechanism to offer succinct advice when it fails to converge and returns unevaluated.

An alternative view of the performance on that example is how the two systems behaved when the author threw a nasty problem at them -- likely deliberately in a naive way. On the one hand I think that Maple does have more documented controls than the author claimed to have found (he hinted that he tried stuff), but on the other hand Mathematica was helpful.

Hmm. Why isn't Warning a help alias for WARNING.

acer