Ask Your Question
1

canonicalize_radical() and evaluation return different answers

asked 4 years ago

Alex Karenin gravatar image

updated 4 years ago

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.

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
3

answered 4 years ago

FrédéricC gravatar image

updated 4 years ago

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.

Preview: (hide)
link

Comments

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

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 4 years ago )

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: 4 years ago

Seen: 564 times

Last updated: Nov 01 '20