############################################################################
#                                                                          #
#  http://www.research.att.com/~njas/sequences/a080069.py.txt              #
#                                                                          #
#  Coded by Antti Karttunen (His-firstname.His-surname@gmail.com),         #
#  August-September 2006.                                                  #
#                                                                          #
#  This file contains the Python-functions that compute the sequences      #
#  A080069, A080070, A122229, A122232, A122235, A122239, A122242, A122245  #
#                             found in                                     #
#     Neil Sloane's On-Line Encyclopedia of Integer Sequences (OEIS)       #
#                            available at                                  #
#             http://www.research.att.com/~njas/sequences/                 #
#                                                                          #
#  Copy of THIS source file is also available at:                          #
#  http://www.iki.fi/kartturi/matikka/Nekomorphisms/a080069.py             #
#                                                                          #
#  This Python-code is in Public Domain and runs (at least)                #
#  in Python 2.4.2 (#67, Sep 28 2005, 12:41:11)                            #
#  (See www.python.org)                                                    #
#                                                                          #
#  For drawing images, you need also PIL-libary,                           #
#  (tested with version 1.1.5, on Python 2.4.)                             #
#  See: http://www.pythonware.com/library/pil/handbook/index.htm           #
#                                                                          #
#  Last edited  Dec 20 2006 by Antti Karttunen.                            #
#                                                                          #
############################################################################

from math import *



# floor(log(n)/log(2)) does not return correct results after n grows over certain limit!

def A000523(n):
  '''Log_2(n) rounded down. We return -1 for n=0, although -infinity would be correct.'''
  i = -1
  while 0 != n:
    i += 1
    n >>= 1
  return(i)

def A007088(n):
  '''Converts n to binary form. Or equally: nth decimal number using no other digits than 0 and 1.'''
  if 0 == n: return(0)
  else: return((n%2)+10*A007088(n/2))

# Note that A057548(n) = A080300(A057547(n)) = A080300(A079946(A014486(n)))

def A079946(n):
  '''Surround n's binary expansion with 1...0. Result's binary expansion of n has form 11**...*0.'''
  return(2*((1<<(1+A000523(n)))+n))


# Actually, the next function is A036044: write in binary, complement, reverse.
# and A057164(n) = A080300(A036044(A014486(n))) = A080300(A056539(A014486(n)))

# For these routines in C, see http://www.research.att.com/~njas/sequences/a089408.c.txt


def A030101(a):
  '''Reverse a's binary expansion.'''
  b = 0
  while 0 != a:
     b <<= 1
     b |= ((a)&1)
     a >>= 1
  return(b)


def A036044(a):
  '''Reverse and complement a totally balanced binary string.'''
  b = 0
  while 0 != a:
     b <<= 1
     b |= ((~a)&1)
     a >>= 1
  return(b)


def A000265(n):
  '''Largest odd divisor of n; or odd part of n.'''
  if 0 == n: return(n)
  while 0 == (n%2): n >>= 1
  return(n)

def A006519(n):
  '''Highest power of 2 dividing n: 1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,16,...'''
  if 0 == n: return(n)
  p = 1
  while 0 == (n%2):
    n >>= 1
    p <<= 1
  return(p)

def A007814(n):
  '''Exponent of highest power of 2 dividing n (the binary carry sequence).'''
  if 0 == n: return(n)
  p = 0
  while 0 == (n%2):
    n >>= 1
    p += 1
  return(p)

def A036987(n):
  '''Fredholm-Rueppel sequence. a(n)=1 iff n = 2^m - 1.'''
  if 0 == (n & (n+1)): return(1)
  else: return(0)

## Note: A125974(10) should return 12, not 4!
def A125974(n):
  '''Function whose restriction to A014486 induces Kreweras bijection A125976.'''
  if 0 == n: return(n)

  chosen = A000265(n)         # Initially ones, get rid of lsb-0's.
  others = n >> A007814(n+1)  # Initially zeros, get rid of lsb-1's.
  s = 0     # the resulting sum
  b = n%2   # n's parity.
  p = 1     # powers of two.

  while (chosen != 0) or (others != 0):
    if (1 == chosen) or (1 == A036987(chosen+1)):  # Last one or zero at hand.
       chosen = others
       others = 0
       nb = 1 - b
    elif (0 == (chosen%4)) or (3 == (chosen%4)):   # Source run continues, dest changes.
       tmp = chosen
       chosen = others
       others = tmp >> 1
       nb = 1 - b
    elif (1 == (chosen%4)): # Source run changes, from ones to zeros, skip past zeros.
       chosen = A000265(chosen-1)
       nb = b
    else: # i.e. if (2 == (chosen%4)) source run changes, from zeros to ones, skip past ones.
       chosen = chosen >> A007814(chosen+2)
       nb = b

    s += b*p
    p <<= 1
    b = nb

  return(s)



def tb_A057164(a):
  return(A036044(a))

def tb_A057163(a):
  '''Reflect a binary tree represented as a totally balanced binary string.
     Keep two "parallel" stacks, the other for the reconstructed
     totally balanced binary strings, and the other for their
     corresponding sizes.
     Start scanning the argument "a" from the end, pushing
     zeros (leaves) to the stack, and joining to topmost
     subtrees (-binary strings) popped from the stack
     when 1 is encountered.'''


  tree_stack = [0] # The last leaf is implicit, not marked in a.
  size_stack = [1] # And its size is 1.

  while (0 != a):
     if(0 == (a&1)): # Push zeros (leaves) to stack.
         tree_stack.append(0)
         size_stack.append(1)
     else: # It's 1, join two branches in swapped order.
         leftchild  = tree_stack.pop()
         rightchild = tree_stack.pop()
         leftsize   = size_stack.pop()
         rightsize  = size_stack.pop()

         size_stack.append(1+leftsize+rightsize)
         tree_stack.append((1 << (leftsize+rightsize)) + (rightchild << leftsize) + leftchild)
     a >>= 1

  return(tree_stack.pop() >> 1) # Discard the last leaf by halving.


def tb_A057117_aux(n,i,r):
  '''Apply the gatomorphism A057117 to A014486-codes. Aux. function that does most of the work.'''

  if(0 == ((n)&1)): return(0)

  c = i
  j = 1
  while j <= r:
     c += (n & 1)
     n >>= 1
     j += 1

  w = c << 1  # w = 2*c
  i <<= 1     # i = 2*i

# Now w = twice the count of ones on preceding row, the width of the next one.
# n points to the beginning of the next row.

  c = 0
  j = 1
  while j <= i:
     c += (n & 1)
     n >>= 1
     j += 1

# Now the 1-bit at the beginning of n is "c":th 1 in whole n (zero-based).

  x = tb_A057117_aux(n,c,(w-(j-1)))
  y = tb_A057117_aux(n>>1,c+(n&1),(w-j))

  if(0 == y): ywidth = 1
  else: ywidth = A000523(y)+1

  if(0 == x): xwidth = 1
  else: xwidth = A000523(x)+1

  return((1 << (ywidth+xwidth)) + (x << ywidth) + y)


def tb_A057117(a):
  '''Apply the gatomorphism A057117 to A014486-codes.'''
  return(tb_A057117_aux(A030101(a),0,1) >> 1)




def tb_A082356(a):
  '''Implement the gatomorphism A082356 on A014486-codes.
     Keep two "parallel" stacks, the other for the reconstructed
     totally balanced binary strings, and the other for their
     corresponding sizes.
     Start scanning the argument "a" as base-4 number, from the
     least-significant end, ...
     See the function quatexpA->parenthesization in gatorank.scm
     and http://www.research.att.com/~njas/sequences/A085184'''

  if(0 == a): return(a)

  a >>= 1 # Discard the least-significant bit (0) before starting.

  tree_stack = [4] # The last two leaves are implicit, 100 in binary.
  size_stack = [3] # And their total size in bits is 3.


  while (a > 2):
     if(0 == (a&3)): # Double-zeros stand for double-leaves (\/), just push them to stack.
         leftchild  = 0
         leftsize   = 1
         rightchild = 0
         rightsize  = 1

     elif(1 == (a&3)): # (case 01)
         rightchild = tree_stack.pop()
         rightsize  = size_stack.pop()
         leftchild  = 0
         leftsize   = 1

     elif(2 == (a&3)): # (case 10)
         rightchild = 0
         rightsize  = 1
         leftchild  = tree_stack.pop()
         leftsize   = size_stack.pop()

     else: # It's 3 (case 11), join two branches in normal order.
         leftchild  = tree_stack.pop()
         leftsize   = size_stack.pop()
         rightchild = tree_stack.pop()
         rightsize  = size_stack.pop()

     size_stack.append(1+leftsize+rightsize)
     tree_stack.append((1 << (leftsize+rightsize)) + (leftchild << rightsize) + rightchild)

     a >>= 2 # Next digit in base-4.

  return(tree_stack.pop() >> 1) # Discard the last leaf by halving.



def tb_A074684(a):
  '''Implement the gatomorphism A074684 on A014486-codes. A variant of above.'''

  if(0 == a): return(a)

  a >>= 1 # Discard the least-significant bit (0) before starting.

  tree_stack = [4] # The last two leaves are implicit, 100 in binary.
  size_stack = [3] # And their total size in bits is 3.


  while (a > 2):
     if(0 == (a&3)): # Double-zeros stand for double-leaves (\/), just push them to stack.
         leftchild  = 0
         leftsize   = 1
         rightchild = 0
         rightsize  = 1

     elif(1 == (a&3)): # (case 01)
         rightchild = 0
         rightsize  = 1
         leftchild  = tree_stack.pop()
         leftsize   = size_stack.pop()

     elif(2 == (a&3)): # (case 10)
         rightchild = tree_stack.pop()
         rightsize  = size_stack.pop()
         leftchild  = 0
         leftsize   = 1

     else: # It's 3 (case 11), join two branches in normal order.
         leftchild  = tree_stack.pop()
         leftsize   = size_stack.pop()
         rightchild = tree_stack.pop()
         rightsize  = size_stack.pop()

     size_stack.append(1+leftsize+rightsize)
     tree_stack.append((1 << (leftsize+rightsize)) + (leftchild << rightsize) + rightchild)

     a >>= 2 # Next digit in base-4.

  return(tree_stack.pop() >> 1) # Discard the last leaf by halving.



def tb_A082358(a):
  '''Implement the gatomorphism A082358 on A014486-codes. A variant of above.'''
  return(tb_A057163(tb_A082356(a)))

def tb_A082360(a):
  '''Implement the gatomorphism A082360 on A014486-codes. A variant of above ones.'''
  return(tb_A057163(tb_A074684(a)))


########################################################################

# For testing:

seqA014486 = [0,2,10,12,42,44,50,52,56,170,172,178,180,184,202,
              204,210,212,216,226,228,232,240,682,684,690,692,
              696,714,716,722,724,728,738,740,744,752,810,812,
              818,820,824,842,844,850,852,856,866,868,872,880,
              906,908,914,916,920,930,932,936,944,962,964,968,976,992]

# For example:
# l = [tb_A057117(n) for n in seqA014486]
# [seqA014486.index(l[i]) for i in range(len(l))]
# Gives the signature-permutation of A057117:
# [0, 1, 2, 3, 4, 5, 7, 8, 6, 9, 10, 12, 13, 11, 17, 18, 21, 22, 20, 14, 15, 16, 19, 23, 24, 26, 27, 25, 31, 32, 35, 36, 34, 28, 29, 30, 33, 45, 46, 49, 50, 48, 58, 59, 63, 64, 62, 54, 55, 57, 61, 37, 38, 40, 41, 39, 44, 47, 42, 43, 56, 60, 51, 52, 53]

# l = [tb_A082356(n) for n in seqA014486]
# [seqA014486.index(l[i]) for i in range(len(l))]

# l = [tb_A082358(n) for n in seqA014486]
# [seqA014486.index(l[i]) for i in range(len(l))]

# l = [tb_A082360(n) for n in seqA014486]
# [seqA014486.index(l[i]) for i in range(len(l))]

# l = [tb_A074684(n) for n in seqA014486]
# [seqA014486.index(l[i]) for i in range(len(l))]

########################################################################


def take(n,g):
  '''Returns a list composed of n next elements returned by generator g. Inspired by Haskell'''
  z = []
  if 0 == n: return(z)
  for x in g:
    z.append(x)
    if n > 1:
      n = n-1
    else:
      return(z)


# Note that A080068 can be also obtained as iteration of A072795 o A057506.
# This works as A057506(n) = A057163(A057164(n), and instead of adding a right stick ./
# to the binary tree at the middle, we can add a left stick \. into it in the beginning
# or at the end of each iteration.

# E.g. take(10,genA080069()) gives: [2, 10, 44, 178, 740, 2868, 11852, 47522, 190104, 735842]

def genA080069():
    '''Yield successive terms of A080069, starting from A080069(1)=2.'''
    i = 2
    while True:
       yield i
       i = tb_A057163(A079946(tb_A057164(i)))

def genA080070():
    '''Yield successive terms of A080070, starting from A080070(1)=10.'''
    i = 2
    while True:
       yield A007088(i)
       i = tb_A057163(A079946(tb_A057164(i)))


def genA122229():
    '''Yield successive terms of A122229, starting from A122229(1)=2.'''
    i = 2 # Boring to look at, but included for completeness!
    while True:
       yield i
       i = A079946(tb_A057117(i))


def genA122232():
    '''Yield successive terms of A122232, starting from A122232(1)=42.'''
    i = 42 # Chaotic...
    while True:
       yield i
       i = A079946(tb_A057117(i))


def genA122235():
    '''Yield successive terms of A122235, starting from A122235(1)=44.'''
    i = 44 # Similar looking. Should compute for starting values 50 and 52 also.
    while True:
       yield i
       i = A079946(tb_A057117(i))


def genA122239():
    '''Yield successive terms of A122239, starting from A122239(1)=52.'''
    i = 52 # Similar looking. Should compute for starting value 50 also.
    while True:
       yield i
       i = A079946(tb_A057117(i))

# Neither A082356 nor A074684 dissipate the change.
# But A082358 is surely interesting!

def genA122242():
    '''Yield successive terms of A122242, starting from A122242(1)=42.'''
    i = 42
    while True:
       yield i
       i = A079946(tb_A082358(i))


def genA122245():
    '''Yield successive terms of A122245, starting from A122245(1)=44.'''
    i = 44 # Similar looking. Should compute for starting values 50 and 52 also.
    while True:
       yield i
       i = A079946(tb_A082358(i))


def genA0new1():
    '''Yield successive terms of A0new1, starting from A0new1(1)=2.'''
    i = 2 # Very regular, like some class-1 1D-CA.
    while True:
       yield i
       i = A079946(tb_A082360(i))



########################################################################
#
#
# Part for drawing the above sequence in "Wolframesque 1D-CA manner".
#
# C.f. also Symbolic Systems in chapter 3 "The World of Simple Programs"
# in Stephen Wolfram, A New Kind of Science, Wolfram Media, 2002;
# pp. 102-104, 896-898
#
#
########################################################################

# Uses the PIL-libary, version 1.1.5, on Python 2.4.
# See: http://www.pythonware.com/library/pil/handbook/index.htm
# The graphics interface uses the same coordinate system as PIL itself,
# with (0, 0) in the upper left corner.

import Image, ImageDraw, ImageFont

def draw_bin_string(draw,row,scale,width,height,ymargin,binstr):
  x = (width-1) - ( (width-scale*(A000523(binstr)+1)) / 2 ) # Get it into center.
# x = (width/2) + (scale*row) - 1 # Simpler!
  y = scale*(row - 1) + ymargin
  black = (0,0,0)
  white = (256,256,256)
  pixrange = range(scale)
  while 0 != binstr:
    if (1 == (binstr%2)): color = black # 1's are black.
    else:                 color = white # 0's are white.
    for x_off in pixrange:
     for y_off in pixrange: draw.point([(x-x_off,y+y_off)], fill=color)
    binstr >>= 1
    x -= scale


def draw_up_to_n(gen,upto_n,scale,filebase,captext):
  '''Draw binary strings produced by generator gen, up to upto_n:th row, saving
     the image to the file, with additional caption "captext", if present.'''

  bfileout = open("b"+filebase+".txt",'w')

  row = 1
  xmargin = 0
  ymargin = 1


# Take the first integer returned by the generator gen:
  for binstr in gen:
    bfileout.write(str(row)+" "+str(binstr)+"\n")
    break

  firstwid = (A000523(binstr)+1)

  width  = 2*(scale*upto_n) + (scale*firstwid) + 2*xmargin
  height = (scale*upto_n) + 2*ymargin
  image = Image.new("RGB",(width,height),(128,000,000)) # Nice red background
  draw = ImageDraw.Draw(image)
  draw_bin_string(draw,row,scale,width,height,ymargin,binstr)

  row += 1

# x = (width-1) - (width-(A000523(binstr)+1))/2 # Get it into center.

# And then the rest:
  for binstr in gen:
    bfileout.write(str(row)+" "+str(binstr)+"\n")
    draw_bin_string(draw,row,scale,width,height,ymargin,binstr)
    row += 1
    if row > upto_n: break

  bfileout.close()

  if(captext):
    # font = ImageFont.load("some_larger_font.pil") # But we don't have it!
    font = ImageFont.load_default()
    draw.text((20,10), captext, fill=(0,0,0), font=font) # Text in black.
    draw.text((20,25), "First "+str(upto_n)+" iterations, 1 bit = "
                       + str(scale) + "x" + str(scale) + " pixels.",
                       fill=(0,0,0), font=font)

  del draw
  image.save("a" + filebase + "_" + str(upto_n) + ".png","png")



# Very fine texture, like a pyramid of sand polished by the wind!
def do_it_for_A080069(upto_n,scale):
  draw_up_to_n(genA080069(),upto_n,scale,"080069","See: http://www.research.att.com/~njas/sequences/A080069")

# Boringly regular:
def do_it_for_A122229(upto_n,scale):
  draw_up_to_n(genA122229(),upto_n,scale,"122229","See: http://www.research.att.com/~njas/sequences/A122229")

# Chaotic:
def do_it_for_A122232(upto_n,scale):
  draw_up_to_n(genA122232(),upto_n,scale,"122232","See: http://www.research.att.com/~njas/sequences/A122232")

def do_it_for_A122235(upto_n,scale):
  draw_up_to_n(genA122235(),upto_n,scale,"122235","See: http://www.research.att.com/~njas/sequences/A122235")

# Are the "circles" real?
def do_it_for_A122239(upto_n,scale):
  draw_up_to_n(genA122239(),upto_n,scale,"122239","See: http://www.research.att.com/~njas/sequences/A122239")

def do_it_for_A122242(upto_n,scale):
  draw_up_to_n(genA122242(),upto_n,scale,"122242","See: http://www.research.att.com/~njas/sequences/A122242")

def do_it_for_A122245(upto_n,scale):
  draw_up_to_n(genA122245(),upto_n,scale,"122245","See: http://www.research.att.com/~njas/sequences/A122245")



########################################################################
#                                                                      #
#  For the comparison, here are Scheme-implementations of some of the  #
#  Catalan automorphisms defined above.                                #
#                                                                      #
########################################################################


# (define (*a082356! s) ;; *a082352! applied to recursion case 0.
#   (cond ((pair? s) (*a082352! s) (*a082356! (car s)) (*a082356! (cdr s))))
#   s
# )
# 
# 
# 
# ;; D   C
# ;;  \ /
# ;;   P   B         D  C B  A  []  B        B   A
# ;;    \ /           \ / \ /    \ /          \ /         and by default:
# ;;     M   A   -->   Q   N      M   A   -->  N  []     []  A       []  A
# ;;      \ /           \ /        \ /          \ /       \ /   -->   \ /
# ;;       X             Y          X            Y         X           Y
# 
# 
# (define (*a082352! s)
#   (cond ((not (pair? s)) s)
#         ((not (pair? (car s))) s)
#         ((not (pair? (caar s))) (swap! (robr! s)))
#         (else (robr! s))
#   )
# )
# 

# And relevant code from gatorank.scm

# (define (binexp->parenthesization n)
#    (let loop ((n n) (stack (list (list))))
#        (cond ((zero? n) (car stack))
#              ((zero? (modulo n 2))
#                 (loop (floor->exact (/ n 2)) (cons (list) stack))
#              )
#              (else
#                 (loop (floor->exact (/ n 2)) (cons2top! stack))
#              )
#        )
#    )
# )
# 
# 
# ;; Experimental, quaternary zigzag-tree code, variant A:
# ;; Differs from A057117 first time at position 56,
# ;; whence this has 42 while it has 44.
# ;; This is now stored in OEIS as A082356.
# ;; (map CatalanRankGlobal (map parenthesization->binexp (map quatexpA->parenthesization (map A014486 (iota0 64)))))
# ;; (0 1 2 3 4 5 7 8 6 9 10 12 13 11 17 18 21 22 20 14 15 16 19 23 24 26 27 25 31 32 35 36 34 28 29 30 33 45 46 49 50 48 58 59 63 64 62 54 55 57 61 37 38 40 41 39 42 43 44 47 51 52 56 60 53)
# 
# (define (quatexpA->parenthesization n)
#    (if (zero? n) (list)
#        (let loop ((n (floor->exact (/ n 2)))
#                   (stack (list (cons (list) (list))))
#                  )
#            (cond ((< n 2) (car stack))
#                  (else
#                    (case (modulo n 4)
#                      ((0) ;; 00
#                         (loop (floor->exact (/ n 4))
#                               (cons (cons (list) (list)) stack)
#                         )
#                      )
#                      ((1) ;; 01
#                         (loop (floor->exact (/ n 4))
#                               (cons2top! (cons (list) stack))
#                         )
#                      )
#                      ((2) ;; 10
#                         (loop (floor->exact (/ n 4))
#                               (flip!topmost (cons2top! (cons (list) stack)))
#                         )
#                      )
#                      ((3) ;; 11
#                         (loop (floor->exact (/ n 4))
#                               (cons2top! stack)
#                         )
#                      )
#                    )
#                  )
#            )
#        )
#    )
# )
# 
# ;; (map CatalanRankGlobal (map parenthesization->binexp (map quatexpB->parenthesization (map A014486 (iota0 64)))))
# ;; (0 1 3 2 8 7 5 4 6 22 21 18 17 20 13 12 10 9 11 15 14 19 16 64 63 59 58 62 50 49 46 45 48 55 54 61 57 36 35 32 31 34 27 26 24 23 25 29 28 33 30 41 40 38 37 39 52 51 60 56 43 42 47 44 53)
# ;; Is equal to A074684.
# ;; Variant B, with the roles of 01 and 10 swapped: (but 11 is not flipped!)
# (define (quatexpB->parenthesization n)
#    (if (zero? n) (list)
#        (let loop ((n (floor->exact (/ n 2)))
#                   (stack (list (cons (list) (list))))
#                  )
#            (cond ((< n 2) (car stack))
#                  (else
#                    (case (modulo n 4)
#                      ((0) ;; 00
#                         (loop (floor->exact (/ n 4))
#                               (cons (cons (list) (list)) stack)
#                         )
#                      )
#                      ((1) ;; 01
#                         (loop (floor->exact (/ n 4))
#                               (flip!topmost (cons2top! (cons (list) stack)))
#                         )
#                      )
#                      ((2) ;; 10
#                         (loop (floor->exact (/ n 4))
#                               (cons2top! (cons (list) stack))
#                         )
#                      )
#                      ((3) ;; 11
#                         (loop (floor->exact (/ n 4))
#                               (cons2top! stack)
#                         )
#                      )
#                    )
#                  )
#            )
#        )
#    )
# )
# 

# 
# ;; Convert (a . (b . rest)) --> ((a . b) . rest)
# ;; with no cons cells wasted.
# 
# (define (cons2top! stack)
#   (let ((ex-cdr (cdr stack)))
#       (set-cdr! stack (car ex-cdr))
#       (set-car! ex-cdr stack)
#       ex-cdr
#   )
# )
# 
# (define (flip!topmost stack)
#   (let* ((topmost (car stack))
#          (ex-cdr (cdr topmost))
#         )
#       (set-cdr! topmost (car topmost))
#       (set-car! topmost ex-cdr)
#       stack
#   )
# )
#