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

Question

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

目前，我将数据集存储在列表映射中，每个输出值都用作输入值列表的键。这似乎很慢，并且似乎使用了整个堆栈空间。有没有更有效的方法来创建/存储/访问我的数据集？

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

? sparsearray(2^16,x->x\7);
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

在某种程度上，您看到的行为是意料之中的。 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 2

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

这是使用libpari 函数vec_equiv() 的解决方案，它返回向量中相同对象的等价类。

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.

Answer 3

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

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