@Kitonum 

What version did you use? In Maple 2018/2019 the code does not work due to a well known problem in fsolve.

A workaround for f(1):

F:= s -> int(sqrt(16 + 4*sin(2*t)^2 + 9*cos(3*t)^2), t = 0 .. s)-1:
Optimization:-NLPSolve('F'(s)^2, s=0..10);

[2.88626817600000*10^(-17), [s = 0.203398786591181]]

@tomleslie 

In v426 sol is 0 (the integral in [-Pi,Pi]  is 0), so, incorrect.

@Carl Love 

Nice!
BTW, you mean this?  Rabin-Shallit-Rand-alg-1986.pdf

@Rouben Rostamian  

You are right. When I wrote the comment It was late and I forgot the ilcm in computing the order of the permuataion (I did not use GroupTheory).  Here are the procedures.
 

f:=(n,k) -> piecewise(k<=floor((n+1)/2), n-2*k+2, 2*k-n-1):
P:= n -> [seq(f(n,k), k=1..n)]: # perm
IP:= n -> convert(convert(P(n),Perm)^(-1),list):  # inverse perm
C:= n -> convert(P(n),disjcyc):
N:= n -> ilcm(seq(nops(u), u=C(n))):  # Order
encrypt := s -> cat(convert(s,list)[IP(length(s))][]):
decrypt := s -> cat(convert(s,list)[P(length(s))][]):

But now the cipher seems to be less interesting :-)

@acer 

The method is equivalent to applying (formally) the standard Newton for sqrt(f).

@ActiveUser 

You should try to ask the question more clearly.

This cipher is more interesting than I thought.
Here are some simple facts I found using a few lines of code.
Obviously the cipher reduces to a permutation. Let's denote it by P[n] (n is the number of letters).

1.  For n>4,  Order(P[n]) >= n.
    Order(P[14]) = 14, that is why the Einstein trick works.

2.  The values n <= 1000  satisfying  Order(P[n]) = n are:

1, 2, 3, 5, 6, 9, 11, 14, 18, 23, 26, 29, 30, 33, 35, 39, 41, 50, 51, 53, 65, 69, 74, 81, 83, 86, 89, 90, 95, 98, 99, 105, 113, 119, 131, 134, 135, 146, 155, 158, 173, 174, 179, 183, 186, 189, 191, 194, 209, 210, 221, 230, 231, 233, 239, 243, 245, 251, 254, 261, 270, 273, 278, 281, 293, 299, 303, 306, 309, 323, 326, 329, 330, 338, 350, 354, 359, 371, 375, 378, 386, 393, 398, 410, 411, 413, 414, 419, 426, 429, 431, 438, 441, 443, 453, 470, 473, 483, 491, 495, 509, 515, 519, 530, 531, 543, 545, 554, 558, 561, 575, 585, 593, 606, 611, 614, 615, 618, 629, 638, 639, 641, 645, 650, 651, 653, 659, 683, 686, 690, 713, 719, 723, 725, 726, 741, 743, 746, 749, 755, 761, 765, 771, 774, 779, 783, 785, 791, 803, 809, 810, 818, 831, 833, 834, 846, 866, 870, 873, 879, 891, 893, 911, 923, 930, 933, 935, 938, 939, 950, 953, 965, 974, 975, 986, 989, 993, 998

3.  The maximum order for n<=1000 is for n=995 and equals
33304588136113480079057828402455235901094586193989341420861500394831872
i.e. 2^52 * 3^115 = 3.3304588... * 10^70

Note that if Einstein had a Spanish style name e.g. with 42 letters, the number of "cycles" in your example would have been 131072.

Here it is:

2015x.mw

The code works in Maple 2015 and 2019. I don't see how it could fail in 2017. Just copy in 1D input at a Maple prompt.

@mmcdara 

I know what the packing problem is; I did myself such packings using global optimization programs. This one is a (simple) one because a given finite set of balls are "packed"/placed in a given cube. Is it so hard to see this?

I have used Maple 2019 new features. Here is a <2019 version.

restart;  # Versiune Maple < 2018
d:=2; R := [seq(rand(0.1 .. 0.15)(),i=1..20)]; timemax:=100;   # R := [0.1 $ 30]
#d:=3; R := [seq(rand(0.1 .. 0.15)(),i=1..60)]; timemax:=100;
S:=Array(1..1,{1 = [<0.2,0.3>, 0.6]});
#S:=Array(1..1,{1 = [<0.,0.,0>, 0.2]});
####################################
NS:=numelems(S):
Xrnd := r -> LinearAlgebra:-RandomVector(d, generator=-1.+r .. 1.-r): 
nor:=LinearAlgebra:-VectorNorm:
check:=proc(v::Vector, r::positive)
  local w;
  for w in S do
    if nor(w[1]-v,2) < r + w[2] then return false fi
  od;
true;
end:

timelimit(timemax,
proc() local r,v; global NS;
for r in R do v:=S[1][1];
  while not check(v,r) do v:=Xrnd(r) od;
  NS:=NS+1; S(NS) := [v, r];
od end ()
):

if   d=2 then plots:-display(seq( plottools:-disk([seq(s[1])],s[2], color=blue),s=S))  
elif d=3 then plots:-display(seq( plottools:-sphere([seq(s[1])],s[2], style=patchnogrid),s=S))
else "cannot plot" fi;

@sand15 

The radii of the balls must be given, as I did in my solution.
It is actually a sphere packing problem.

@acer 

I agree that it is probably a wise decision, but should be documented.

BTW, the result is hfloat only for a Vector/Matrix, but not Array; even if datatype=algebraic.

@acer 

Is it documented?

UseHardwareFloats;
                            deduced

Vector([5.,1]): lprint(%);
Vector[column](2,{1 = 5., 2 = 1},datatype = anything,storage = rectangular,
order = Fortran_order,shape = [])

1/3 * Array([5.,1]): lprint(%);
Array(1 .. 2,{1 = 1.666666667, 2 = 1/3})

1/3 * Vector([5.,1]): lprint(%);
Vector[column](2,{1 = HFloat(1.66666666666666652), 2 = HFloat(.333333333333333315)},datatype = float[8],storage = rectangular,order = Fortran_order,shape = [])

# I would not expect 1/3 * 1 --> hfloat

@Carl Love 

You are right, sorry. I think I've read (partly) the code of your first proc ;-)