1

I try to compute real 1-Brjuno function for real numbers in x in [0, 1] range

  G(x):= block(
            [g],
        if (x=0) 
            then g : 0
            else g : float(1/x - floor(1/x)),
        return( g)
        )$
    
   bet(j,x) := block(
        [r],
        if (j=-1) 
            then r:1
            else r : product (G(x)^i, i, 0, j),
        return( r)
        )$
    
    A(i,x) := bet(i-1,x)*log(1/G(x)^i)$
    B(x):= sum(A(i,x), i, 0,10)$

There are numerical errors :

                        0
    expt: undefined: 0.0

caused by A function for some x values : 1/2,1/5,1/10,1/3,1/9,1/4,0

How can I solve this problem ?

Adam
  • 1,254
  • 12
  • 25
  • 1
    Try `simp:false; tellsimp(0^0, 1); tellsimp(0.0^0, 1); simp:true;`. – Robert Dodier Sep 22 '21 at 21:15
  • @RobertDodier expt: undefined: 0 to a negative exponent. – Adam Sep 23 '21 at 15:33
  • Okay, that means you have something/0 at some point. Is there a limit that applies in that case? or some other kind of special case? You might need to find out whether "something" is itself zero or nonzero. By the way, I tried a few examples, and it seems like the "0 to a negative exponent" occurs for "nice" values like 1/2, but for random values I don't get the error. – Robert Dodier Sep 23 '21 at 16:42
  • @RobertDodier I use simply plot2d(B(x),[x,0,1]); then there are no errors but the diagram looks different then in the paper. So I changed to discrete set of points ( rational umbers) in [0,1] range, – Adam Sep 23 '21 at 16:45

2 Answers2

1

The Brjuno function has " logarithmic singularities at rational points" so using simple "

G(x):= block(

    [g],
    if (x=0) 
        then g : 0
        else g : float(1/x - floor(1/x)),
    return( g)
    )$
    

beta(j,x) := block(

    [r],
    if (j=-1) 
        then r:1
        else r : product (G(x)^i, i, 0, j),
    return( r)
    )$
    
    
A(i,x) := beta(i-1,x)*log(1/G(x)^i)$
B(x):= sum(A(i,x), i, 0,10)$
plot2d(B(x),[x,0,1], [y, -1, 10]);

gives plot : enter image description here

But the diagram is not the same as Figure 1 ( left ) in the paper.

Adam
  • 1,254
  • 12
  • 25
1

This produces a plot similar to the paper. I used an equation from wikipedia and a continued fraction approximation:

cflength: 10 $
B(x):= if length(x)=1 then 0 else (local(x), x[1]: 0, -log(cf2num(x))+cf2num(x)*B(rest(x))) $
cf2num(e):=ev(cfdisrep(e), numer, infeval) $
x0: makelist(i*sqrt(2)/1421, i, 1, 1000), numer $
x: map(cf, x0) $
y: map(B, x) $
draw2d(points(x0, y)) $

enter image description here

slitvinov
  • 5,693
  • 20
  • 31