## 10455 Reputation

6 years, 124 days

## MaplePrimes Activity

These are Posts that have been published by vv

## A useful exercise for Maple programmers...

Maple

The Maple help contains a nice example of a puzzle solver named alphametic,  see  ?Iterator,Permute.
The goal is to determine the distinct digits represented by letters which satisfy a given  equation. The provided solved example is:

"16540*781 = 12904836 + 12904"
i.e.  m=1, a=6, p=5 etc.

It's worth studying the program (written as a module) because it includes a few useful techniques.
I would suggest the following exercise for a Maple (young) programmer.

1.  When solving the "classical" puzzle  "FORTY+TEN+TEN=SIXTY", instead of the correct answer "29786+850+850=31486",   you will see
"2978Y+850+850=3148Y".
Try to find what's going on and correct the code.

2. The solutions given by alphametic include numbers beginning with a zero.
Execute e.g. .
Modify the code to produce only standard numbers.

## Great Hall Floor and Maple...

Maple

I have recently visited the Queen's House at Greenwich  (see wiki),  an  important building in British architectural history (17th century).
I was impressed by the Great Hall Floor, whose central geometric decoration seems to be generated by a Maple program :-)

Here is my code for this. I hope you will like it.

```restart;
with(plots): with(plottools):
n:=32: m:=3:#    n:=64: m:=7:

a[0], b[0] := exp(-Pi*I/n), exp(Pi*I/n):
c[0]:=b[0]+(a[0]-b[0])/sqrt(2)*exp(I*Pi/4):
for k to m+1 do
c[k]:=a[k-1]+b[k-1]-c[k-1];
b[k]:=c[k]*(1+exp(2*Pi*I/n))-b[k-1];
a[k]:=conjugate(b[k])  od:
b[-1]:=c[0]*(1+exp(2*Pi*I/n))-b[0]:
a[-1]:=conjugate(b[-1]):
c[-1]:=a[-1]+b[-1]-c[0]:
seq( map[inplace](evalf@[Re,Im], w), w=[a,b,c] ):
Q:=polygonplot([seq([a[k],c[k],b[k],c[k+1]],k=0..m-1), [c[m],a[m],b[m]], [a[-1],b[-1],c[0]]]):
display(seq(rotate(Q, 2*k*Pi/n, [0,0]),k=0..n-1), disk([0,0],c[m][1]/3), axes=none, size=[800,800], color=black);
```

## A limit from an undergraduate competitio...

Maple

At a recent undegraduate competition the students had to compute the following limit

 > Limit( n * Diff( (exp(x)-1)/x, x\$n), n=infinity ) assuming x<>0;
 (1)

Maple is able to compute the symbolic n-fold derivative and I hoped that the limit will be computed at once.

Unfortunately it is not so easy.
Maybe someone finds a more more straightforward way.

 > restart;
 > f := n * diff( (exp(x)-1)/x, x\$n );
 (2)
 > limit(%, n=infinity);
 (3)
 > simplify(%) assuming x>0;
 (4)

So, Maple cannot compute directly the limit.

 > convert(f, Int) assuming n::posint;
 (5)
 > J:=simplify(%)  assuming n::posint;
 (6)
 > L:=convert(J, Int) assuming n::posint;
 (7)
 > L:=subs(_k1=u, L);
 (8)

Now it should be easy, but Maple needs help.

 > with(IntegrationTools):
 > L1:=Change(L, u^n = t, t) assuming n::posint;
 (9)
 > limit(L1, n=infinity);  # OK
 (10)
 > ####################################################################

Note that the limit can also be computed using an integration by parts, but Maple refuses to finalize:

 > Parts(L, exp(u*x)) assuming n::posint;
 (11)
 > simplify(%);
 (12)
 > limit(%, n=infinity);
 (13)
 > value(%);  # we are almost back!
 (14)
 >

## add, floats, and Kahan sum...

Maple

add, floats, and Kahan sum

I found an intresting fact about the Maple command add for floating point values.
It seems that add in this case uses a summation algorithm in order to reduce the numerical error.
It is probably the Kahan summation algorithm (see wiki), but I wonder why this fact is not documented.

Here is a simple Maple procedure describing and implementing the algorithm.

 > restart;
 > Digits:=15;
 (1)
 > KahanSum := proc(f::procedure, ab::range)   local S,c,y,t, i;      # https://en.wikipedia.org/wiki/Kahan_summation_algorithm S := 0.0;              # S = result (final sum: add(f(n), n=a..b)) c := 0.0;              # c = compensation for lost low-order bits. for i from lhs(ab) to rhs(ab) do     y := f(i) - c;          t := S + y;                   c := (t - S) - y;             S := t;                   od;                          return S end proc:

Now, a numerical example.

 > f:= n ->  evalf(1/(n+1/n^3+1) - 1/(n+1+1/(n+1)^3+1));
 (2)
 > n := 50000; K := KahanSum(f, 1..n);
 (3)
 > A := add(f(k),k=1..n);
 (4)
 > s:=0.0:  for i to n do s:=s+f(i) od: 's' = s;
 (5)
 > exact:=( 1/3 - 1/(n+1+1/(n+1)^3+1) );
 (6)
 > evalf( [errK = K-exact, errA = A-exact, err_for=s-exact] );
 (7)
 > evalf[20]( [errK = K-exact, errA = A-exact, err_for=s-exact] );
 (8)
 >

## Bug corrected in inttrans[hilbert]...

Maple

There is a bug in inttrans:-hilbert:

 > restart;
 > inttrans:-hilbert(sin(a)*sin(t+b), t, s); # should be: sin(a)*cos(s+b);   expand(%);
 (1)
 > ########## correction ##############
 > `inttrans/expandc` := proc(expr, t) local xpr, j, econst, op1, op2;       xpr := expr;             for j in indets(xpr,specfunc(`+`,exp)) do           econst := select(type,op(j),('freeof')(t));           if 0 < nops(econst) and econst <> 0 then               xpr := subs(j = ('exp')(econst)*combine(j/('exp')(econst),exp),xpr)           end if       end do;       for j in indets(xpr,{('cos')(linear(t)), ('sin')(linear(t))}) do           if type(op(j),`+`) then               op1:=select(has, op(j),t); ##               op2:=op(j)-op1;            ##               #op1 := op(1,op(j));               #op2 := op(2,op(j));               if op(0,j) = sin then                   xpr := subs(j = cos(op2)*sin(op1)+sin(op2)*cos(op1),xpr)               else                   xpr := subs(j = cos(op1)*cos(op2)-sin(op1)*sin(op2),xpr)               end if           end if       end do;       return xpr end proc:
 > #######################################
 > inttrans:-hilbert(sin(a)*sin(t+b), t, s); expand(%);
 (2)
 >

 1 2 3 4 5 6 Page 3 of 6
﻿