Bill Allombert on Sat, 22 Nov 2014 20:58:23 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Growing ordered set


On Fri, Nov 21, 2014 at 03:36:40PM +0100, Bill Allombert wrote:
> On Fri, Nov 21, 2014 at 12:18:29PM +0100, Bill Allombert wrote:
> > On Fri, Nov 21, 2014 at 01:04:26AM +0100, Karim Belabas wrote:
> > Well, but writing everything in C is cheating...
> > 
> > Running my code through gp2c leads
> > 
> > ? setrand(1);#test(10^11)
> > %1 = 537776
> > ? ##
> >   ***   last result computed in 1,484 ms.
> > ? setrand(1);#test2(10^11)
> > %2 = 537776
> > ? ##
> >   ***   last result computed in 41,908 ms.
> > 
> > Thus we could provide a C interface to t_LIST trees that would be not much slower
> > than your code.
> 
> To give an example, the attached C file implement a function treeadd()
> that can be used in GP as follow:

The attached file implement self-balancing AVL trees.
The interface is the same, but this is not subject to the worst-cases issues
of the previous (not self-balancing) version.

You can compile the file with
gp2c-run treeavl.c

Cheers,
Bill.
/*-*- compile-command: "/usr/bin/gcc -c -o ../tree.gp.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I"/home/bill/pari/amd64/include" ../tree.gp.c && /usr/bin/gcc -o ../tree.gp.so -shared -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -Wl,-shared ../tree.gp.o -lc -lm -L/home/bill/pari/amd64/lib -lpari"; -*-*/
#include <pari/pari.h>
#include <pari/paripriv.h>
/*
GP;install("treeadd","lWG","treeadd","./treeavl.so");
*/

#define value(i)  gmael(list_data(T),(i),1)
#define left(i)   mael3(list_data(T),(i),2,1)
#define right(i)  mael3(list_data(T),(i),2,2)
#define height(i) mael3(list_data(T),(i),2,3)

static long
treeheight(GEN T, long i)
{ return i? height(i): 0; }

static long
new_leaf(GEN T, GEN x)
{ 
  pari_sp av = avma;
  listput(T, mkvec2(x, mkvecsmall3(0,0,1)), 0);
  avma = av;
  return lg(list_data(T))-1;
}

static void
fix_height(GEN T, long x)
{
  height(x) = maxss(treeheight(T,left(x)), treeheight(T,right(x)))+1;
}
static long
treebalance(GEN T, long i)
{
  return i ? treeheight(T,left(i)) - treeheight(T,right(i)): 0;
}

static long
rotright(GEN T, long y)
{
  long x = left(y);
  long t = right(x);
  right(x) = y;
  left(y)  = t;
  fix_height(T, y);
  fix_height(T, x);
  return x;
}

static long
rotleft(GEN T, long x)
{
  long y = right(x);
  long t = left(y);
  left(y)  = x;
  right(x) = t;
  fix_height(T, x);
  fix_height(T, y);
  return y;
}

static long
treeinsert(GEN T, GEN x, long i, long *d)
{
  long b, c;
  if (i==0 || !list_data(T))
    return new_leaf(T, x);
  c = cmp_universal(x, value(i));
  if (c < 0)
  {
    long s = treeinsert(T, x, left(i), d);
    if (s < 0) return s;
    left(i) = s;
  }
  else if (c > 0)
  {
    long s = treeinsert(T, x, right(i), d);
    if (s < 0) return s;
    right(i) = s;
  }
  else return -i;
  fix_height(T, i);
  b = treebalance(T, i);
  if (b > 1)
  {
    if (*d > 0)
      left(i) = rotleft(T, left(i));
    return rotright(T, i);
  }
  if (b < -1) 
  {
    if (*d < 0)
      right(i) = rotright(T, right(i));
    return rotleft(T, i);
  }
  *d = c;
  return i; 
}

long
treeadd(GEN T, GEN x)
{
  GEN d;
  long c = 0;
  long r = treeinsert(T, x, 1, &c);
  if (r < 0) return -r;
  if (r == 1) return 0;
  d = list_data(T);
  /* By convention we want the root to be 1 */
  swap(gel(d,1), gel(d,r)); 
  if (left(1) == 1) left(1)=r;
  else if (right(1) == 1) right(1)=r;
  else pari_err_BUG("treeadd");
  return 0;
}