如何在 Pari/GP 中表示稀疏数组？

Question

我有一个 returns 整数值到整数输入的函数。输出值相对稀疏；对于输入值 1....2^16，该函数仅 returns 大约 2^14 个唯一输出。我想创建一个数据集，让我可以快速找到产生任何给定输出的输入。

目前，我将我的数据集存储在列表映射中，每个输出值都用作输入值列表的键。这看起来很慢并且似乎使用了整个堆栈 space。 create/store/access 我的数据集是否有更有效的方法？

已添加： 事实证明，我的 sparesearray() 函数所花费的时间在输出值（即键）与输入值（存储在列表中的值）的比率上变化很大。下面是一个需要很多列表的函数所花费的时间，每个列表只有几个值：

? sparsearray(2^16,x->x);
time = 126 ms.

这是一个只需要几个列表的函数所花费的时间，每个列表都有很多值：

? sparsearray(2^12,x->x%7);
time = 218 ms.
? sparsearray(2^13,x->x%7);
time = 892 ms.
? sparsearray(2^14,x->x%7);
time = 3,609 ms.

如你所见，时间呈指数增长！

这是我的代码：

\ sparsearray takes two arguments, an integer "n"  and a closure "myfun", 
\ and returns a Map() in which each key a number, and each key is associated 
\ with a List() of the input numbers for which the closure produces that output. 
\ E.g.:
\ ? sparsearray(10,x->x%3)
\ %1 = Map([0, List([3, 6, 9]); 1, List([1, 4, 7, 10]); 2, List([2, 5, 8])])
sparsearray(n,myfun=(x)->x)=
{
    my(m=Map(),output,oldvalue=List());
    for(loop=1,n,
        output=myfun(loop);                      
        if(!mapisdefined(m,output), 
        /* then */
            oldvalue=List(),
        /* else */    
            oldvalue=mapget(m,output));
        listput(oldvalue,loop);
        mapput(m,output,oldvalue));
    m
}

Answer 1

在使用 Map() 创建的地图上使用 mapput、mapget 和 mapisdefined 方法。如果需要多个维度，则使用多项式或向量键。

我想这就是您已经在做的事情，我不确定是否有更好的方法。你有一些代码吗？根据个人经验，2^16 个值和 2^14 个键在速度或内存方面应该不是问题 - 在您的实现中可能会进行一些不必要的复制。

Answer 2

在某种程度上，您所看到的行为是可以预料的。 PARI 似乎按值而不是引用传递列表和映射，除了用于操作它们的特殊内置函数。这可以通过创建像 mylistput(list,item)=listput(list,item); 这样的包装函数来看出。当您尝试使用此功能时，您会发现它不起作用，因为它是在列表的副本上运行。可以说，这是 PARI 中的一个错误，但也许他们有他们的理由。这种行为的结果是每次您将元素添加到存储在地图中的列表之一时，整个列表都会被复制，可能复制两次。

以下是避免此问题的解决方案。

sparsearray(n,myfun=(x)->x)=
{
   my(vi=vector(n, i, i)); \ input values
   my(vo=vector(n, i, myfun(vi[i]))); \ output values
   my(perm=vecsort(vo,,1)); \ obtain order of output values as a permutation
   my(list=List(), bucket=List(), key);
   for(loop=1, #perm, 
      if(loop==1||vo[perm[loop]]<>key, 
          if(#bucket, listput(list,[key,Vec(bucket)]);bucket=List()); key=vo[perm[loop]]);
      listput(bucket,vi[perm[loop]])
   );

   if(#bucket, listput(list,[key,Vec(bucket)])); 
   Mat(Col(list))
}

输出是一个与地图格式相同的矩阵 - 如果您更喜欢地图，则可以使用 Map(...) 进行转换，但您可能需要一个矩阵进行处理，因为没有内置的在地图上获取键列表的函数。

我对上面的内容做了一些修改，试图在 C# 中制作更类似于 GroupBy 的内容。（一个可以用于很多事情的函数）

VecGroupBy(v, f)={
   my(g=vector(#v, i, f(v[i]))); \ groups
   my(perm=vecsort(g,,1)); 
   my(list=List(), bucket=List(), key);
   for(loop=1, #perm, 
      if(loop==1||g[perm[loop]]<>key, 
          if(#bucket, listput(list,[key,Vec(bucket)]);bucket=List()); key=g[perm[loop]]);
      listput(bucket, v[perm[loop]])
   );
   if(#bucket, listput(list,[key,Vec(bucket)])); 
   Mat(Col(list))
}

您可以像 VecGroupBy([1..300],i->i%7) 那样使用它。

Answer 3

由于垃圾收集的方式，没有好的本地 GP 解决方案，因为在 GP 的内存模型中必须限制通过引用传递参数（从 2.13 版本开始，使用 ~ 修饰符，但不适用于地图组件）。

这是一个使用 libpari 函数 vec_equiv() 的解决方案，其中 returns 向量中相同对象的等价性类。

install(vec_equiv,G);
sparsearray(n, f=x->x)=
{
  my(v = vector(n, x, f(x)), e  = vec_equiv(v));
  [vector(#e, i, v[e[i][1]]), e];
}

? sparsearray(10, x->x%3)
%1 = [[0, 1, 2], [Vecsmall([3, 6, 9]), Vecsmall([1, 4, 7, 10]), Vecsmall([2, 5, 8])]]

（您有 3 个值对应于给定的 3 组索引）

行为如预期的那样是线性的

 ? sparsearray(2^20,x->x%7);
 time = 307 ms.
 ? sparsearray(2^21,x->x%7);
 time = 670 ms.
 ? sparsearray(2^22,x->x%7);
 time = 1,353 ms.

如何在 Pari/GP 中表示稀疏数组？

How do I represent sparse arrays in Pari/GP?

pari-gp