

Comes across a fun math puzzle. :)
Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
cin(X) must be smooth function, and at least 10 digits accurate.
In other words, cin(cin(cin(X))) / sin(X) = 1.00000 00000
What is the value of cin(2.019) ?


Op Do. 4 Apr. 2019 om 01:14 het Albert Chan < [hidden email]> geskryf:
>
> Comes across a fun math puzzle. :)
>
> Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
>
> cin(X) must be smooth function, and at least 10 digits accurate.
>
> In other words, cin(cin(cin(X))) / sin(X) = 1.00000 00000
>
> What is the value of cin(2.019) ?
Since you post it on this list, one assumes that the calculation musr
be done in Lua :)
Am I allowed the luxury of requiring luafft (available from LuaRocks)?


> On Apr 4, 2019, at 3:03 AM, Dirk Laurie < [hidden email]> wrote:
>
> Op Do. 4 Apr. 2019 om 01:14 het Albert Chan < [hidden email]> geskryf:
>>
>> Comes across a fun math puzzle. :)
>>
>> Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
>>
>> cin(X) must be smooth function, and at least 10 digits accurate.
>>
>> In other words, cin(cin(cin(X))) / sin(X) = 1.00000 00000
>>
>> What is the value of cin(2.019) ?
>
> Since you post it on this list, one assumes that the calculation musr
> be done in Lua :)
>
> Am I allowed the luxury of requiring luafft (available from LuaRocks)?
>
Yes, all help is allowed.
I use Mathematica to design and test for accuracy, before coding cin(X) in Lua.


On Thu, Apr 4, 2019 at 2:14 AM Albert Chan wrote: Comes across a fun math puzzle. :)
Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
The most straightforward approach is
to build Taylor series of cin(x) one term at a time.
local maclaurin_of_cin
do
local c, d, s, a = {[0] = 1}, {[0] = 1}, 1
local function g(k, m)
if k == 0 or m < 0 then
return m == 0 and 1 or 0
else
local i = k..";"..m
local r = a[i]
if not r then
r = 0
for j = 0, #c do
r = r + c[j] * g(k1, mj)
end
a[i] = r
end
return r
end
end
local function f(o)
local r = 0
for j = 0, #c do
r = r + o[j] * g(2*j+1, #c+1j)
end
return r
end
function maclaurin_of_cin(k)
for n = #c + 1, k do
a = {}
local e, h = f(c), f(d)
s, a = s/(2*n)/(2*n+1)
local t = (she)/3
assert(math.abs(t) < 0.056)
c[n], d[n] = t, e + 2*t
end
return c[k]
end
end
local function cin(x)
local r, p, s, n, R = x, x, x*x, 0
repeat
R, n, p = r, n+1, p*s
r = r + maclaurin_of_cin(n) * p
until r == R
return r
end
 (0.7) < x < 0.7
local x = 0.7
print("x = "..x)
print("cin(x) = "..cin(x))
print("cin(cin(x)) = "..cin(cin(x)))
print("cin(cin(cin(x))) = "..cin(cin(cin(x))))
print("sin(x) = "..math.sin(x))
This implementation works only if x is in the range from (0.7) to (+0.7)
That's because of floating point arithmetic is approximate.
Exact arithmetic of fractions (with arbitrary long numerator and denominator) must be used instead.


> On Apr 4, 2019, at 5:25 PM, Egor Skriptunoff < [hidden email]> wrote:
>
> This implementation works only if x is in the range from (0.7) to (+0.7)
> That's because of floating point arithmetic is approximate.
> Exact arithmetic of fractions (with arbitrary long numerator and denominator) must be used instead.
>
Building Taylor coefficents on the fly ! Amazing !
How do you figure the valid range of 0.7 to 0.7
However, this only solve 1/2 of the puzzle, unable to do cin(2.019)
Here is a hint for the second half of the puzzle:
sin(cin(X)) = cin(sin(X))


> On Apr 4, 2019, at 12:03 AM, Dirk Laurie < [hidden email]> wrote:
>
> Op Do. 4 Apr. 2019 om 01:14 het Albert Chan < [hidden email]> geskryf:
>>
>> Comes across a fun math puzzle. :)
>>
>> Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
>>
>> cin(X) must be smooth function, and at least 10 digits accurate.
>>
>> In other words, cin(cin(cin(X))) / sin(X) = 1.00000 00000
>>
>> What is the value of cin(2.019) ?
>
> Since you post it on this list, one assumes that the calculation musr
> be done in Lua :)
>
> Am I allowed the luxury of requiring luafft (available from LuaRocks)?
>
There is of course the “cheat” answer, which highlights how careful you need to be when setting puzzles:
do
local c = 0
function cin(x) c = c + 1; return ((c % 3) == 0) and sin(x) or x end
end
—Tim


There is of course the “cheat” answer, which highlights how careful you need to be when setting puzzles:
do
local c = 0
function cin(x) c = c + 1; return ((c % 3) == 0) and sin(x) or x end
end
Note that if x is a constant, your cin() function will not return a constant value (except when x==0), because it will return x two times before returning sin(x), and we have sin(x)!=x for any x!=0. The problem should have stated that the function cin(x) must return a value that does not depend on any other external or internal variable, but depending only on x, i.e. the mathematical problem should have been to find an "application" and not a "function" (the Lua "functions" are not restricted to be strictly mathematical applications, and can then return values depending on other variables, but the problem was not given as requiring Lua or any specific programming language).
Also the problem does not clearly state the meaning of "smooth". I think it should have been said to be  an application on a specified input set (e.g. all real numbers representable in the programming language, or a subset of it, but the second part of the problem clearly stated that it must handle input values outside the range where sin(x) has the same sign as the input in radians).  and the application must be defined, continuous and infinitely derivable on all values of the input range.
Note: the range (0.7 .. +0.7) is in fact (pi/4 .. pi/4): on this range all derivate orders of sin(x) have a constant sign and any input of sin(x) in that range will output a value within that range, and with the same sign as the input.
The second half of the problem using x=2.019 is outside the range and the derivates are changing their sign, si it would be difficult to reach the convergence. However, sin(x) has a Taylor development with the almost correct sign if you develop it to the 10th degree to reach the 10digit precision of the problem and it is enough to compute the Taylor coefficients of cin(x) with a precision of 11 digits for degrees 0 and 1, and 10 digits for degrees 2 or higher.
The coefficient of degree 0 is of course 0. But for getting a warranty of 10 digits of precision of the input, the needed precision for the coefficient depends on the width of the input range: this precision of 10 digits is UNREACHABLE if the input set is (INFINITE, +INFINITE) and we need a bounded number of Taylor coefficient: with large values of x, 10 degrees is not enough for the Taylor polynom. Now suppose you don't precompute the Taylor coefficients but just use an binary search loop, it's not warrantied we get a convergence because the signs are constantly changing for high values of x. We can however ensure the convergence if the input is restricted to (pi, +pi) and then can reply to the problem with x=2.019.
Anyway, this is not really a problem for Lua, any programming language could be used (including an Excel worksheet, or Javascript, Java, C, C++, C#, VB, PHP, Python, and even Bash or Powershell...) provided they have a numeric floatting point type capable of representing numbers with 10 digits of precision (the IEEE 32bit "float" is not enough, you need the IEEE 64bit "double" or better). So it's not really a problem in scope with this list.


On Apr 4, 2019, at 5:25 PM, Egor Skriptunoff < [hidden email]> wrote:
> local function cin(x)
> local r, p, s, n, R = x, x, x*x, 0
> repeat
> R, n, p = r, n+1, p*s
> r = r + maclaurin_of_cin(n) * p
> until r == R
> return r
> end
How does your maclaurin_of_cin work ?
It looks like magic to me :)
To handle big cin argument, I use identity cin(x) = asin(cin(sin(x)))
So, adding this line to your cin(x) is all that is needed.
if math.abs(x) >= 0.7 then return math.asin(cin(math.sin(x))) end
With this 1 line patch:
x = 2.019
cin(x) = 1.02692 331869 35764
cin(cin(x)) = 0.95662 892999 61186
cin(cin(cin(x))) = 0.90122 698939 98124
math.sin(x) = 0.90122 698939 98126
Amazing accuracy !


a = false
cin = function(x)
if not a then a = (math.sin(x)/x)^(1.0/3.0) end
return a*x
end
x = 2.019
print(cin(x),'check:',cin(cin(cin(x)))/math.sin(x))
> 1.543010716654 check: 1.0
:)
Em qua, 3 de abr de 2019 às 19:14, Albert Chan
< [hidden email]> escreveu:
>
> Comes across a fun math puzzle. :)
>
> Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
>
> cin(X) must be smooth function, and at least 10 digits accurate.
>
> In other words, cin(cin(cin(X))) / sin(X) = 1.00000 00000
>
> What is the value of cin(2.019) ?
>
>

Rodrigo Azevedo Moreira da Silva


This is my cin(x) version, by reducing x to below 0.1 radian.
If x >= 0.1, it's reduced value = (sin(0.1) to 0.1) ~ (0.099833, 0.1)
Taylor series last term is adjusted to reduce errors for this tight domain.
local abs, sin, asin = math.abs, math.sin, math.asin
function cin(x)
local n = 0  count nested sin's
while abs(x) >= 0.1 do x=sin(x); n=n+1 end
local x2 = x*x
x = x  x*x2*(1/18 + x2*(7/1080 + x2*(643/408240 + x2*0.0004635)))
for i=1,n do x=asin(x) end
return x
end
Result has about 13 digits accuracy:
x = 2.019
cin(x) = 1.02692 331869 35365
cin(cin(x)) = 0.95662 892999 61478
cin(cin(cin(x))) = 0.90122 698939 98782
math.sin(x) = 0.90122 698939 98126


Op Do. 4 Apr. 2019 om 23:25 het Egor Skriptunoff
< [hidden email]> geskryf:
>
> On Thu, Apr 4, 2019 at 2:14 AM Albert Chan wrote:
>>
>> Comes across a fun math puzzle. :)
>> Design a function cin(X), such that cin(cin(cin(X))) = sin(X)
>>
>
> The most straightforward approach is
> to build Taylor series of cin(x) one term at a time.
[code snipped]
> This implementation works only if x is in the range from (0.7) to (+0.7)
> That's because of floating point arithmetic is approximate.
> Exact arithmetic of fractions (with arbitrary long numerator and denominator) must be used instead.
No, the reason is that the convergence radius of the Maclaurin series
is only about 0.7. There are at least two ways around that:
1. Use a nonlinear transformation (epsilon algorithm, Levin's U etc)
to sum the series.
2. Apply analytic continuation. E.g. use this method to calculate cin
and enough derivatives at say pi/6, then calculate a Taylor series
round pi/6, etc.
I can't guarantee that either approach will work. Unfortunately I am
travelling and am typing this in an idle halfhour, so I can't try it
out now.
Problems of this kind are subtle and often have no unique solution.
Additional hypotheses are sometimes needed. The convergence
acceleration method in effect assumes that cin has no singularities
worse than poles near the region of interest. The analytic
continuation approach works if for example cin is analytic in an
ellipse with foci at (0,0) and (2.019,0). There is no reason to
suppose that the two answers will be the same. The very neat trick
which the designer of the puzzle certainly wanted solvers to see, does
not require that cin is analytic anywhere except at 0. This answer may
well be different from the other two.
Personally, I am very fond of the U transform because it produces an
asymptotic series, i.e. the answers initially seem to converge but
then diverge. The user is thus warned that there is some nonpolar
singularity or excessive roundoff propagation, and the point where
this happens is merely a "best practical" answer for the given data.
Sorry, we are straying from Lua, but since for example discussions on
finite state machines have been conducted on the list recently, a
little computational complex analysis can't hurt :)
 Dirk


On Apr 6, 2019, at 2:45 AM, Dirk Laurie < [hidden email]> wrote:
>
>> This implementation works only if x is in the range from (0.7) to (+0.7)
>> That's because of floating point arithmetic is approximate.
>> Exact arithmetic of fractions (with arbitrary long numerator and denominator) must be used instead.
>
> No, the reason is that the convergence radius of the Maclaurin series
> is only about 0.7. There are at least two ways around that:
>
> 1. Use a nonlinear transformation (epsilon algorithm, Levin's U etc)
> to sum the series.
> 2. Apply analytic continuation. E.g. use this method to calculate cin
> and enough derivatives at say pi/6, then calculate a Taylor series
> round pi/6, etc.
>
> I can't guarantee that either approach will work. Unfortunately I am
> travelling and am typing this in an idle halfhour, so I can't try it
> out now.
Good catch !
According to the author of this cin(x) puzzle, Valentin Albiolo,
radius of convergence is actually 0 !
In other words, you cannot keep building more Taylor coefficients,
and expect better accuracy ...
Here is link to official solution:
[VA] Short and Sweet Math Challenge #24, post #28
http://www.hpmuseum.org/forum/thread12656post114578.html#pid114578


> On Apr 4, 2019, at 5:25 PM, Egor Skriptunoff < [hidden email]> wrote:
>
> The most straightforward approach is
> to build Taylor series of cin(x) one term at a time.
>
> function maclaurin_of_cin(k)
> for n = #c + 1, k do
> a = {}
> local e, h = f(c), f(d)
> s, a = s/(2*n)/(2*n+1)
> local t = (she)/3
> assert(math.abs(t) < 0.056)
> c[n], d[n] = t, e + 2*t
> end
> return c[k]
> end
Say, I wanted a function, din(x), such that din(din(din(din(x)))) = sin(x)
What is required to modify above to do maclaurin_of_din(k) ?


Op Di. 9 Apr. 2019 om 00:09 het Albert Chan < [hidden email]> geskryf:
>
>
> > On Apr 4, 2019, at 5:25 PM, Egor Skriptunoff < [hidden email]> wrote:
> >
> > The most straightforward approach is
> > to build Taylor series of cin(x) one term at a time.
> >
> > function maclaurin_of_cin(k)
> > for n = #c + 1, k do
> > a = {}
> > local e, h = f(c), f(d)
> > s, a = s/(2*n)/(2*n+1)
> > local t = (she)/3
> > assert(math.abs(t) < 0.056)
> > c[n], d[n] = t, e + 2*t
> > end
> > return c[k]
> > end
>
> Say, I wanted a function, din(x), such that din(din(din(din(x)))) = sin(x)
>
> What is required to modify above to do maclaurin_of_din(k) ?
I can't speak for the code, but the usual method is equivalent to
Newton's method for roots, applied to functions. It will work for any
function s whose Maclaurin series starts at x.
The wellknown algorithm for square roots of real numbers is
x ← x + (a/xx)/2.
For nth roots, it becomes
x ← x + (a/x^(n1)x)/n
For nth root of a function, it becomes:
f ← f + (g(g(...(s)))f)/n  n1 applications of g
where g is the inverse function of f.
You can do this analytically, but usually for one step only; the
inverse function soon starts to become intractable.
Examples for din:
f(x)=x, g(x)=x, next approximation is x + (sin(x)x)/3
f(x) = 2*sin(x/2), g(x)=2*asin(x/2), next approximation is
2*sin(x/2) + (2*asin(2*asin(2*asin(sin(x)/2)/2)/2)2*sin(x/2))/4
This second approximation is not bad at all, applied four times
starting at 1 gives 0.84465, should be 0.84147. You can draw a
reasonable graph with that. Applied four times at pi/2, gives 1.058,
still not outrageously bad.
What we do instead is to work with a power series for f; then the
reversal of the series is a truncated power series for g. We add just
one term to both series at a time. You already know up to x^k the
identity holds, so you compare coefficients of x^(k+1). I think that
is the basis for Egor's code, but my letters f and g probably mean
different things to his f and g.


On Tue, Apr 9, 2019 at 1:09 AM Albert Chan wrote:
> The most straightforward approach is
> to build Taylor series of cin(x) one term at a time.
>
> function maclaurin_of_cin(k)
> for n = #c + 1, k do
> a = {}
> local e, h = f(c), f(d)
> s, a = s/(2*n)/(2*n+1)
> local t = (she)/3
> assert(math.abs(t) < 0.056)
> c[n], d[n] = t, e + 2*t
> end
> return c[k]
> end
Say, I wanted a function, din(x), such that din(din(din(din(x)))) = sin(x)
I don't use smart methods like Newton's (mentioned by Dirk).
My implementation is dumb and straightforward.
I solve trivially constructed system of equations to find Tailor coefficients of din().
(I simply substitute one Tailor series into another.) If n first coefficients are already known, the formula for the (n+1)th one is not very complex.
local maclaurin_of_din
do
local c, d, s, a = {[0] = 1}, {[0] = 1}, 1
local function g(k, m, o)
if k == 0 or m < 0 then
return m == 0 and 1 or 0
else
local i = k..";"..m
local r = a[i]
if not r then
r = 0
for j = 0, #o do
r = r + o[j] * g(k1, mj, o)
end
a[i] = r
end
return r
end
end
local function f(o)
local r = 0
a = {}
for j = 0, #o do
r = r + o[j] * g(2*j+1, #o+1j, o)
end
a = nil
return r
end
function maclaurin_of_din(k)
for n = #c + 1, k do
local e, h = f(c), f(d)
s = s/(2*n)/(2*n+1)
local t = (sh2*e)/4
assert(t*t < 1)
c[n], d[n] = t, e + 2*t
end
return c[k]
end
end
local function din(x)
local r, p, s, n, R = x, x, x*x, 0
repeat
R, n, p = r, n+1, p*s
r = r + maclaurin_of_din(n) * p
until r == R
return r
end
local x = 0.7
print("x = "..x)
print("din(x) = "..din(x))
print("din(din(x)) = "..din(din(x)))
print("din(din(din(x))) = "..din(din(din(x))))
print("din(din(din(din(x)))) = "..din(din(din(din(x)))))
print("sin(x) = "..math.sin(x))


On Tuesday, April 9, 2019, 5:55:08 PM EDT, Egor Skriptunoff <[hidden email]> wrote:
> I don't use smart methods like Newton's (mentioned by Dirk). > My implementation is dumb and straightforward. > I solve trivially constructed system of equations to find Tailor coefficients of din(). > (I simply substitute one Tailor series into another.) > If n first coefficients are already known, the formula for the (n+1)th one is not very > complex.
Thank you for the code. Does not look dumb / straightforward to me ! Still don't understand what d table is used for ...
Anyway, I think I also did it ... uh, by guessing :)
First 20 din(x) coefficients matched your version (ignoring rounding errors).
A bonus is that coefficients generation is almost twice as fast.
Using your code as template, just replace f() and maclaurin_of_din():
local function f(coef) local r = 0 for j = 1, #c do r = r + coef[j] * g(2*j+1, #c+1j, c) end return r end local d2 = {[0] = 1}
function maclaurin_of_din(k) for n = #c + 1, k do a = {} local r, r2, r3 = f(c), f(d), f(d2) s, a = (2*n)*(2*n+1)*s local t = (1/srr2r3)/4 assert(t*t < 1) c[n], d[n], d2[n] = t, r + 3*t, r2 + 2*t end return c[k] end

