When I qualified my suggestion with "if efficiency matters" and "or critical code" I had in mind other scenarios like pure computation, as opposed to this task of printing and viewing. Sorry if I wasn't clear enough about that.

acer

When I qualified my suggestion with "if efficiency matters" and "or critical code" I had in mind other scenarios like pure computation, as opposed to this task of printing and viewing. Sorry if I wasn't clear enough about that.

acer

@Robert Israel That's like what you answered in another Question posted by the same member.

@Robert Israel That's like what you answered in another Question posted by the same member.

@mzh It would take away too much from Maple if the constructing such things produces an error. For one thing, such constructions may have different interpretations, and thus should be "OK" to exist as intermediate objects.

Consider, as just one example

myequation := 32 = 64;

                            32 = 64

evalb(myequation); # of course

                             false

myequation mod 2;

                             0 = 0

evalb(myequation mod 2);

                              true

So, construction of `myequation` might get done as part of some involved computation involving several values, and only later might the intermediate objects (equations) get tested modulo some value. It wouldn't be possible to code that, if the construction immediately produced an error. In general, there are many more such examples.

@mzh It would take away too much from Maple if the constructing such things produces an error. For one thing, such constructions may have different interpretations, and thus should be "OK" to exist as intermediate objects.

Consider, as just one example

myequation := 32 = 64;

                            32 = 64

evalb(myequation); # of course

                             false

myequation mod 2;

                             0 = 0

evalb(myequation mod 2);

                              true

So, construction of `myequation` might get done as part of some involved computation involving several values, and only later might the intermediate objects (equations) get tested modulo some value. It wouldn't be possible to code that, if the construction immediately produced an error. In general, there are many more such examples.

Another way to handle this is to ensure that the arguments are row-column. For this column-column example one could transpose the first Vector. This causes LinearAlgebra:-Multiply to do the work, without the conjugation,

<a, b, c> . <d, e, f>;
                        _     _     _  
                        a d + b e + c f

<a, b, c>^%T . <d, e, f>;

                        a d + b e + c f

<a|b|c> . <d,e,f>; # same orientations as above

                        a d + b e + c f

Transposition can be done inplace (highly memory and speed efficient) on a Vector, either by using the LinearAlgebra:-Transpose command or using the `subtype` parameter of the `rtable_options` command

Other straightforward ways to do it include making the assumption that the relevant names are purely real (or to use a command like `evalc` which does that in a blanketing way),

evalc(<a, b, c> . <d, e, f>);

                        a d + b e + c f

<a, b, c> . <d, e, f> assuming real;

                        a d + b e + c f

<a, b, c> . <d, e, f> assuming a::real, b::real, c::real;

                        a d + b e + c f

acer

I overlooked the addition of Jacobi based svd as new functions dgejsv and dgesvj in version 3.2 or LAPACK. (See section 9.5 of these release v.3.2 notes.)

If you were to compile CLAPACK 3.2.1, and its accompanying reference BLAS, as a pair of .dll libraries with those entry point functions exported, then you should be able to use them directly from Maple using the `define-external` command. (I might try building this one first, perhaps. The only pre-built libraries for this that I know of are .lib, and you'd need dynamic .dll instead. So you'd need to build it yourself.)

acer

I overlooked the addition of Jacobi based svd as new functions dgejsv and dgesvj in version 3.2 or LAPACK. (See section 9.5 of these release v.3.2 notes.)

If you were to compile CLAPACK 3.2.1, and its accompanying reference BLAS, as a pair of .dll libraries with those entry point functions exported, then you should be able to use them directly from Maple using the `define-external` command. (I might try building this one first, perhaps. The only pre-built libraries for this that I know of are .lib, and you'd need dynamic .dll instead. So you'd need to build it yourself.)

acer

Very nice explanation, with more sophistication indeed than in my answer.

And thanks for `style=point`, which looks much better! I was trying to steer him towards bounding the error (or "last" term), but didn't want to give away what might be a homework question too easily.

This reminds me of a usenet thread from 2010, which revolved around confusion due to Maple's use of double-precision `sin` which is (on some platforms) not always the same as what the operating system's runtime provides. Axel provided a fascinating reference, during that discussion, BTW.

Any readers interested in the range reduction from values of much higher multiples of Pi into (0,Pi/4) could start with Mueller's great book, Elementary Functions, Algorithms, and Implementation. (Axel, we've discussed that book once before, I think...)

Posted code in reponses here by me and Axel contain conversion to Horner's form, and therein lies a whole other deep area. Which reminds me of something else I had sheets on for potential blogging: Horner's form for Matrix polynomial evaluation. And then there is Estrin's form for use with Threads...

acer

Very nice explanation, with more sophistication indeed than in my answer.

And thanks for `style=point`, which looks much better! I was trying to steer him towards bounding the error (or "last" term), but didn't want to give away what might be a homework question too easily.

This reminds me of a usenet thread from 2010, which revolved around confusion due to Maple's use of double-precision `sin` which is (on some platforms) not always the same as what the operating system's runtime provides. Axel provided a fascinating reference, during that discussion, BTW.

Any readers interested in the range reduction from values of much higher multiples of Pi into (0,Pi/4) could start with Mueller's great book, Elementary Functions, Algorithms, and Implementation. (Axel, we've discussed that book once before, I think...)

Posted code in reponses here by me and Axel contain conversion to Horner's form, and therein lies a whole other deep area. Which reminds me of something else I had sheets on for potential blogging: Horner's form for Matrix polynomial evaluation. And then there is Estrin's form for use with Threads...

acer

@Pavel Kriz By computing the log-base-10 of the error, one can see how many decimal places of the result are accurate.

I chose 16 since evalhf and compiled double-precision C is close to 15-16 decimal places of width for the mantissa. A natural question is: what 16 decimal digit formula can be written out to C or Fortran and then used to get roughly 16 correct decimal digits for all x in (0,Pi/4).

BTW, I don't know if this is assigned work. Sometimes, when computing a series approximation one can consider bounding the error and in turn finding the number of terms to ensure that,

@Pavel Kriz By computing the log-base-10 of the error, one can see how many decimal places of the result are accurate.

I chose 16 since evalhf and compiled double-precision C is close to 15-16 decimal places of width for the mantissa. A natural question is: what 16 decimal digit formula can be written out to C or Fortran and then used to get roughly 16 correct decimal digits for all x in (0,Pi/4).

BTW, I don't know if this is assigned work. Sometimes, when computing a series approximation one can consider bounding the error and in turn finding the number of terms to ensure that,

@alfred Yes, of course, my x=stage2 (or x=stage2[1] with the 'location' option) is not right at all, when computing stage1. It could be repaired by isolating G=g for x, evaluating that with G=stage2, and then using that for the x value to compute stage1. One reason I won't write that out right now is that in general you may not be able to isolate G=g for x (either explicitly, or uniquely).

Here's another shot. The first is a different symbolic approach (even though I just said why symbolics might not do, in general).

And here's what might be a correction to the earlier numeric approach, keeping in mind that what was wrong before was using y=v(x) making the same kind of oops.

restart:

f := (.5+x-y)^2;

g := (.5-x)^2*(.7-y)^2*(.3-y)^2;

                                               2
                             f := (0.5 + x - y) 

                                  2          2          2
                    g := (0.5 - x)  (0.7 - y)  (0.3 - y) 

stage2:=maximize(maximize(g,x=0..1),y=0..1,location);

        stage2 := 0.01102500000, {[{y = 0.}, 0.01102500000], [{y = 1.}, 0.01102500000]}

seq(T[1], T in stage2[2]);

                             {y = 0.}, {y = 1.}

seq(maximize(eval(f,op(T[1])),x=0..1,location), T in stage2[2]);

           2.250000000, {[{x = 1.}, 2.250000000]}, 0.2500000000, 

             {[{x = 0.}, 0.2500000000], [{x = 1.}, 0.2500000000]}

v := proc(X)
      if not type(X,numeric) then 'procname'(args);
      else Optimization:-NLPSolve(eval(g,x=X),y=0..1,maximize)[1];
      end if;
     end proc:

#plot(v(x),x=0..1);

st2num:=Optimization:-NLPSolve(v,0..1,maximize,method=branchandbound);

                    st2num := [0.0110249987483025, [0.]]

H:=proc(X)
    if not type(X,numeric) then 'procname'(args);
    else eval(f,{x=X,y=Optimization:-NLPSolve(v,0..1,maximize,
                                              method=branchandbound)[2][1]});
    end if;
   end proc:

#plot(H(x),x=0..1);  # takes a while

Optimization:-NLPSolve(H(x),{x>=0,1>=x},maximize,method=sqp);

Warning, no iterations performed as initial point satisfies first-order conditions
                       [2.25000000000000000, [x = 1.]]

Come to think of it... the definition of procedure `v` should probably use a global method just in case there are local maxima for some input X which don't happen to be the global maximum.

v := proc(X)
      if not type(X,numeric) then 'procname'(args);
      else Optimization:-NLPSolve(eval(g,x=X),y=0..1,maximize,
                                  method=branchandbound)[1];
      end if;
     end proc:

@alfred Yes, of course, my x=stage2 (or x=stage2[1] with the 'location' option) is not right at all, when computing stage1. It could be repaired by isolating G=g for x, evaluating that with G=stage2, and then using that for the x value to compute stage1. One reason I won't write that out right now is that in general you may not be able to isolate G=g for x (either explicitly, or uniquely).

Here's another shot. The first is a different symbolic approach (even though I just said why symbolics might not do, in general).

And here's what might be a correction to the earlier numeric approach, keeping in mind that what was wrong before was using y=v(x) making the same kind of oops.

restart:

f := (.5+x-y)^2;

g := (.5-x)^2*(.7-y)^2*(.3-y)^2;

                                               2
                             f := (0.5 + x - y) 

                                  2          2          2
                    g := (0.5 - x)  (0.7 - y)  (0.3 - y) 

stage2:=maximize(maximize(g,x=0..1),y=0..1,location);

        stage2 := 0.01102500000, {[{y = 0.}, 0.01102500000], [{y = 1.}, 0.01102500000]}

seq(T[1], T in stage2[2]);

                             {y = 0.}, {y = 1.}

seq(maximize(eval(f,op(T[1])),x=0..1,location), T in stage2[2]);

           2.250000000, {[{x = 1.}, 2.250000000]}, 0.2500000000, 

             {[{x = 0.}, 0.2500000000], [{x = 1.}, 0.2500000000]}

v := proc(X)
      if not type(X,numeric) then 'procname'(args);
      else Optimization:-NLPSolve(eval(g,x=X),y=0..1,maximize)[1];
      end if;
     end proc:

#plot(v(x),x=0..1);

st2num:=Optimization:-NLPSolve(v,0..1,maximize,method=branchandbound);

                    st2num := [0.0110249987483025, [0.]]

H:=proc(X)
    if not type(X,numeric) then 'procname'(args);
    else eval(f,{x=X,y=Optimization:-NLPSolve(v,0..1,maximize,
                                              method=branchandbound)[2][1]});
    end if;
   end proc:

#plot(H(x),x=0..1);  # takes a while

Optimization:-NLPSolve(H(x),{x>=0,1>=x},maximize,method=sqp);

Warning, no iterations performed as initial point satisfies first-order conditions
                       [2.25000000000000000, [x = 1.]]

Come to think of it... the definition of procedure `v` should probably use a global method just in case there are local maxima for some input X which don't happen to be the global maximum.

v := proc(X)
      if not type(X,numeric) then 'procname'(args);
      else Optimization:-NLPSolve(eval(g,x=X),y=0..1,maximize,
                                  method=branchandbound)[1];
      end if;
     end proc: