Postgis中计算某些特殊点的距离时特别慢的诡异问题

背景

在某个系统中，在Postgis中计算某些特殊点的距离时特别慢，比正常慢了40几倍，刚开始一直百思不得其解。中启乘数科技的工程师通过debug数据库程序，发现是CentOS7.X的glibc中的数学库libm中三角函数cos计算某些值特别慢的问题。

场景重新

由于涉及最终客户的一些信息，我们这里就不给出客户的具体SQL了。我们自己造一个SQL来重现此问题。

建一个存储过程：

CREATE OR REPLACE FUNCTION test_distance(arg_p1 geometry, arg_p2 geometry)
  RETURNS void AS
$BODY$
DECLARE
   v_dist double precision;
BEGIN
FOR i IN 1..10000 LOOP
    v_dist := ST_DistanceSphere (arg_p1, arg_p2);
END LOOP;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE;

上面的存储过程是重复计算两个点之间的距离1万次。

我们随便计算两个点的距离，如下：

postgres=# SELECT test_distance(ST_GeometryFromText('POINT(0 0)', 4326), ST_GeometryFromText('POINT(100 45)', 4326));
 test_distance
---------------
(1 row)
Time: 36.038 ms
postgres=# SELECT test_distance(ST_GeometryFromText('POINT(0 0)', 4326), ST_GeometryFromText('POINT(30 30)', 4326));
 test_distance
---------------
(1 row)
Time: 33.948 ms

可以看出上面的时间只需要30~40ms。

但是当我们换一个特殊的点“POINT(119.297915 42.042785)”时发现非常慢，居然要1.5秒：

postgres=# SELECT test_distance(ST_GeometryFromText('POINT(0 0)', 4326), ST_GeometryFromText('POINT(119.297915 42.042785)', 4326));
 test_distance
---------------
(1 row)
Time: 1548.751 ms (00:01.549)

问题定位

通过gdb去debug PostgreSQL数据库，发现是运行函数sphere_distance慢：

double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
{
    struct  timeval   tv_begin,tv_end;
    int64_t my_usec;
    gettimeofday(&tv_begin,NULL);    
    double d_lon = e->lon - s->lon;
    double cos_d_lon = cos(d_lon);
    double cos_lat_e = cos(e->lat);
    double sin_lat_e = sin(e->lat);
    double cos_lat_s = cos(s->lat);
    double sin_lat_s = sin(s->lat);
    double a1 = POW2(cos_lat_e * sin(d_lon));
    double a2 = POW2(cos_lat_s * sin_lat_e - sin_lat_s * cos_lat_e * cos_d_lon);
    double a = sqrt(a1 + a2);
    double b = sin_lat_s * sin_lat_e + cos_lat_s * cos_lat_e * cos_d_lon;
    gettimeofday(&tv_end,NULL);
    my_usec = (int64_t)tv_end.tv_sec * 1000000LL + (int64_t)tv_end.tv_usec - ((int64_t)tv_begin.tv_sec * 1000000LL + (int64_t)tv_begin.tv_usec);
    g_my_cnt += my_usec;
    return atan2(a, b);
}

进一步debug发现是运行上面的三角函数cos慢。也就是cos函数对于某些特殊的值时非常慢。
为此我们在SQL中直接执行cos函数，发现对于一些特殊值来说会非常慢。

先造测试表：

create table test01(f float8);
insert into test01 select 0.73 from generate_series(1, 10000) as seq;

这时运行cos函数会非常快：

postgres=# explain analyze select sum(cos(f)) from test01;
                                                   QUERY PLAN
----------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=197.55..197.56 rows=1 width=8) (actual time=2.628..2.628 rows=1 loops=1)
   ->  Seq Scan on test01  (cost=0.00..146.70 rows=10170 width=8) (actual time=0.017..1.229 rows=10000 loops=1)
 Planning Time: 0.075 ms
 Execution Time: 2.650 ms
(4 rows)
Time: 3.084 ms

当我们把表中的值改成特殊的值0.73378502495808418时，就会非常慢：

postgres=# update test01 set f=0.73378502495808418;
UPDATE 10000
Time: 23.088 ms
postgres=# explain analyze select sum(cos(f)) from test01;
                                                   QUERY PLAN
----------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=385.67..385.68 rows=1 width=8) (actual time=369.834..369.835 rows=1 loops=1)
   ->  Seq Scan on test01  (cost=0.00..286.78 rows=19778 width=8) (actual time=0.934..2.461 rows=10000 loops=1)
 Planning Time: 0.072 ms
 Execution Time: 369.859 ms
(4 rows)
Time: 370.771 ms

进一步定位发现就是三角函数cos对于一些特殊的值非常慢。写一个C的测试程序testcos.c：

#include <stdio.h>
#include <stdint.h>
double  test_cos(double lat)
{
    double cos_lat_e = cos(lat);
    return cos_lat_e;
}
int main(int argc, char * argv[])
{
    int i;
    double dis;
    double lat;
    int flag = atoi(argv[1]);
    int cnt = atoi(argv[2]);
    uint64_t u64lat;
    if (flag) 
    {
        /*慢*/
        //lat = 0.73378502495808418; 
        u64lat=0x3FE77B2ABB8FAA1C;
        //临时
        //u64lat=0x3FE77B2ABB8FBA1C;
    }
    else
    {   /*快*/
        //lat = 0.73378502670341339;
        u64lat=0x3FE77B2ABC7F8A6C;
    }
    lat = *(double *) &u64lat;
   //printf("size of double=%d\n", sizeof(lat));
   //u64lat = *(uint64_t *)&lat;
   printf("u64lat=%lX\n", u64lat);
   printf("lat=%.15f\n", lat);
   for (i=0; i<cnt; i++)
   {
       dis = test_cos(lat);
   }
   return 0;
}

注意上面把cos计算放在一个单独的函数中，避免编译器把多次循环给优化掉。

编译：

gcc testcos.c -o testcos -lm

[codetest@pgdev gistest]$ time ./testcos 1 50000
lat=0.733785024958084
real    0m1.853s
user    0m1.852s
sys     0m0.000s
[codetest@pgdev gistest]$ time ./testcos 0 50000
lat=0.733785026703413
real    0m0.002s
user    0m0.000s
sys     0m0.001s

testcos中第一个参数为1是慢的情况，第一个参数为0是快的情况。
从上面的测试可以看出不同的值，cos函数执行时间相差上千倍。

上面的测试时是在CentOS7.6上测试的。然后我们又测试了操作系统Rocky Linux8.6，发现没有这个慢的问题。
所以基本断定是glibc的数学库cos的问题。通过查询网上的信息，说这是一个bug:

• https://sourceware.org/bugzilla/show_bug.cgi?id=5781
• https://sourceware.org/bugzilla/show_bug.cgi?id=14412
• https://sourceware.org/bugzilla/show_bug.cgi?id=15936

同时查看glibc中也就是一些特别的值做sin或cos计算时，发现误差过大，然后做了更精确的计算,导致慢了很多。

在CentOS7.X下的glibc版本为glibc-2.17，可以编译升级glibc到2.28的版本，就没有此问题了。注意替换glibc需要小心操作，不小心操作很容易导致操作系统无法启动的问题。

如果不想替换整个操作系统的glibc，可以只替换PostgreSQL程序，让其链接到新版本glibc的数学动态库：

编译glibc-2.28，假设编译完了glib-2.28安装在/home/codetest/glibc-2.28，把postgres依赖的动态库libm.so.6指向glibc-2.28中的libm.so

patchelf --replace-needed libm.so.6 /home/codetest/glibc-2.28/lib/libm-2.28.so /mypg/pgsql/bin/postgres

问题结论

实际上此问题并不是PostgreSQL数据库的问题，也不是PostGIS的问题，而是glibc-2.17版本的问题，在低版本的glibc中，计算一些特殊值的sin或cos时会非常慢，需要特别注意。

glibc是非常底层的库，当这些底层的库有一些不易发现的缺陷时，会导致一些很奇怪的问题。从此例可以看出，随着国产化的深入需要更多的关心一些底层的技术，才能有效保证计算机系统的安全稳定运行。