Ask Your Question
1

canonicalize_radical() and evaluation return different answers

asked 2020-10-30 00:10:19 +0200

Alex Karenin gravatar image

updated 2020-10-30 09:09:47 +0200

FrédéricC gravatar image

This code return proof that f^8(x)=f(f(f(f(f(f(f(f(x))))))))=x if f(x)=(x(sqrt(2)-1)-1)/(x-1):

def yx(x): 
    cst=sqrt(2)-1
    if(x==1):
        return -Infinity
    return (cst*x-1)/(x-1)
tmp=x 
for t in range(1,9): 
    tmp=(yx(tmp)).canonicalize_radical() 
    print(t, tmp.canonicalize_radical())

But this code looks directly at orbit of a point and shows that the statement above is false:

tmp=1.5
lst=[]
for t in range(1,9): 
    tmp=(yx(tmp)).canonicalize_radical() 
    print(t, tmp.canonicalize_radical().n(80))
    lst.append((t,tmp.n()))

By giving the next output:

1 -0.75735931288071485359493
2 0.74754689570642838213503
3 2.7345908033901357840028
4 0.076504843704676803069591
5 1.0485281374238889523706
6 0.13254046199725214558851
7 -0.090964954340477885594991
8 0.35391568118499587774957

Why the output is so different? Bit precision don't seem to be a problem.

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
3

answered 2020-10-30 08:17:05 +0200

FrédéricC gravatar image

updated 2020-11-01 14:04:42 +0200

Use exact fields and polynomial rings, avoid the symbolic ring and numerical reals:

sage: K.<a> = QuadraticField(2)                                                   
sage: x = polygen(K, 'x')                                                          
sage: def yx(x):  
....:     cst = a - 1 
....:     if x == 1: 
....:         return -Infinity 
....:     return (cst * x - 1) / (x - 1) 
....: tmp = x  
....: for t in range(1,9):  
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 ((a - 1)*x - 1)/(x - 1)
2 (a*x - 1)/x
3 (-x - 1)/(x - a - 1)
4 (1/2*a*x - 1)/(x - 1/2*a)
5 ((a + 1)*x - 1)/(x + 1)
6 (-1)/(x - a)
7 (x - 1)/(x - a + 1)
8 x

EDIT: and when using this with numerical input, the answers are more reasonable

sage: tmp = 1.5                                                                 
sage: for t in range(1,9): 
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 -0.757359312880715
2 0.747546895706428
3 2.73459080339014
4 0.0765048437046765
5 1.04852813742386
6 -11.6568542494924
7 0.460495713220364
8 1.50000000000000

The correct answer differs at the 6th value from the wrong computation. Not clear how and why.

EDIT:

Comparing the correct answer and the procedure using the symbolic sqrt(2), one can see precision loss happening, with large exponents for the floating point coefficients of sqrt(2)

(-0.757359312880715, 3.0*sqrt(2) - 5.0)
(0.747546895706428, (-8.0*sqrt(2) + 10.0)/(3.0*sqrt(2) - 6.0))
(2.73459080339014, (-150.0*sqrt(2) + 210.0)/(114.0*sqrt(2) - 162.0))
(0.0765048437046765, (-79524.0*sqrt(2) + 112464.0)/(85176.0*sqrt(2) - 120456.0))
(1.04852813742386, (-25732492128.0*sqrt(2) + 36391239360.0)/(39678297120.0*sqrt(2) - 56113585920.0))
(-11.6568542494924, (-(2.51898761979489e+21)*sqrt(2) + 3.56238645536385e+21)/((7.34086788505619e+21)*sqrt(2) - 1.03815549226356e+22))
(0.460495713220364, ((2.61510083243401e+43)*sqrt(2) - 3.69831106420135e+43)/((2.04721262905725e+44)*sqrt(2) - 2.89519586507424e+44))
(1.50000000000000, ((1.55098758806785e+89)*sqrt(2) - 2.19342768211789e+89)/((1.03399172537856e+89)*sqrt(2) - 1.46228512141192e+89))

So this problem is essentially caused by using a mixture of symbolic sqrt(2) and floating points.

edit flag offensive delete link more

Comments

The reason(s) to avoid SR is(are) distinct of the reason(s) to avoid (approximate) numericals.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-10-30 09:32:19 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2020-10-30 00:10:19 +0200

Seen: 466 times

Last updated: Nov 01 '20