From ebdc2d2115b850384f3910f33fbaefc3d57687f3 Mon Sep 17 00:00:00 2001
From: Tomas Vondra <tomas@2ndquadrant.com>
Date: Sat, 1 Sep 2018 22:46:45 +0200
Subject: [PATCH 2/2] multivariate histograms

---
 doc/src/sgml/catalogs.sgml                       |   14 +-
 doc/src/sgml/func.sgml                           |   85 +
 doc/src/sgml/planstats.sgml                      |  106 +
 doc/src/sgml/ref/create_statistics.sgml          |   29 +-
 src/backend/catalog/system_views.sql             |   11 +
 src/backend/commands/statscmds.c                 |   33 +-
 src/backend/nodes/outfuncs.c                     |    2 +-
 src/backend/optimizer/util/plancat.c             |   44 +-
 src/backend/parser/parse_utilcmd.c               |    2 +
 src/backend/statistics/Makefile                  |    2 +-
 src/backend/statistics/README                    |    4 +
 src/backend/statistics/README.histogram          |  305 +++
 src/backend/statistics/dependencies.c            |    2 +-
 src/backend/statistics/extended_stats.c          |  252 +-
 src/backend/statistics/histogram.c               | 3019 ++++++++++++++++++++++
 src/backend/statistics/mcv.c                     |   94 +-
 src/backend/utils/adt/ruleutils.c                |   10 +
 src/backend/utils/adt/selfuncs.c                 |    7 +-
 src/bin/psql/describe.c                          |    9 +-
 src/include/catalog/pg_cast.dat                  |    4 +
 src/include/catalog/pg_proc.dat                  |   24 +
 src/include/catalog/pg_statistic_ext.h           |    2 +
 src/include/catalog/pg_type.dat                  |    7 +
 src/include/nodes/relation.h                     |    7 +-
 src/include/statistics/extended_stats_internal.h |   15 +
 src/include/statistics/statistics.h              |   63 +-
 src/include/utils/selfuncs.h                     |    4 +
 src/test/regress/expected/create_table_like.out  |    2 +-
 src/test/regress/expected/opr_sanity.out         |    3 +-
 src/test/regress/expected/stats_ext.out          |  209 +-
 src/test/regress/expected/type_sanity.out        |    3 +-
 src/test/regress/sql/stats_ext.sql               |  133 +-
 32 files changed, 4403 insertions(+), 103 deletions(-)
 create mode 100644 src/backend/statistics/README.histogram
 create mode 100644 src/backend/statistics/histogram.c

diff --git a/doc/src/sgml/catalogs.sgml b/doc/src/sgml/catalogs.sgml
index dc7bbe5173..0edd28ad0e 100644
--- a/doc/src/sgml/catalogs.sgml
+++ b/doc/src/sgml/catalogs.sgml
@@ -6571,8 +6571,9 @@ SCRAM-SHA-256$<replaceable>&lt;iteration count&gt;</replaceable>:<replaceable>&l
         An array containing codes for the enabled statistic kinds;
         valid values are:
         <literal>d</literal> for n-distinct statistics,
-        <literal>f</literal> for functional dependency statistics, and
-        <literal>m</literal> for most common values (MCV) list statistics
+        <literal>f</literal> for functional dependency statistics,
+        <literal>m</literal> for most common values (MCV) list statistics, and
+        <literal>h</literal> for histogram statistics
       </entry>
      </row>
 
@@ -6605,6 +6606,15 @@ SCRAM-SHA-256$<replaceable>&lt;iteration count&gt;</replaceable>:<replaceable>&l
       </entry>
      </row>
 
+     <row>
+      <entry><structfield>stxhistogram</structfield></entry>
+      <entry><type>pg_histogram</type></entry>
+      <entry></entry>
+      <entry>
+       Histogram, serialized as <structname>pg_histogram</structname> type.
+      </entry>
+     </row>
+
     </tbody>
    </tgroup>
   </table>
diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
index 69cfe7bbe9..0a5a62685a 100644
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -20992,6 +20992,91 @@ SELECT m.* FROM pg_statistic_ext,
    </para>
   </sect2>
 
+  <sect2 id="functions-statistics-histogram">
+   <title>Inspecting histograms</title>
+
+   <indexterm>
+     <primary>pg_histogram</primary>
+     <secondary>pg_histogram_buckets</secondary>
+   </indexterm>
+
+   <para>
+    <function>pg_histogram_buckets</function> returns a list of all buckets
+    stored in a multi-column histogram, and returns the following columns:
+
+    <informaltable>
+     <tgroup cols="3">
+      <thead>
+       <row>
+        <entry>Name</entry>
+        <entry>Type</entry>
+        <entry>Description</entry>
+       </row>
+      </thead>
+      <tbody>
+       <row>
+        <entry><literal>index</literal></entry>
+        <entry><type>int</type></entry>
+        <entry>index of the item in the histogram buckets</entry>
+       </row>
+       <row>
+        <entry><literal>minvals</literal></entry>
+        <entry><type>text[]</type></entry>
+        <entry>lower boundaries of the histogram bucket</entry>
+       </row>
+       <row>
+        <entry><literal>maxvals</literal></entry>
+        <entry><type>text[]</type></entry>
+        <entry>upper boundaries of the histogram bucket</entry>
+       </row>
+       <row>
+        <entry><literal>nullsonly</literal></entry>
+        <entry><type>boolean[]</type></entry>
+        <entry>flags identifying <literal>NULL</literal>-only dimensions</entry>
+       </row>
+       <row>
+        <entry><literal>mininclusive</literal></entry>
+        <entry><type>boolean[]</type></entry>
+        <entry>flags identifying which lower boundaries are inclusive</entry>
+       </row>
+       <row>
+        <entry><literal>maxinclusive</literal></entry>
+        <entry><type>boolean[]</type></entry>
+        <entry>flags identifying which upper boundaries are inclusive</entry>
+       </row>
+       <row>
+        <entry><literal>frequency</literal></entry>
+        <entry><type>double precision</type></entry>
+        <entry>frequency of this histogram bucket</entry>
+       </row>
+       <row>
+        <entry><literal>density</literal></entry>
+        <entry><type>double precision</type></entry>
+        <entry>density of this histogram bucket (frequency / volume)</entry>
+       </row>
+       <row>
+        <entry><literal>bucket_volume</literal></entry>
+        <entry><type>double precision</type></entry>
+        <entry>volume of this histogram bucket (a measure of size)</entry>
+       </row>
+      </tbody>
+     </tgroup>
+    </informaltable>
+   </para>
+
+   <para>
+    The <function>pg_histogram_buckets</function> function can be used like this:
+
+<programlisting>
+SELECT m.* FROM pg_statistic_ext,
+                pg_histogram_buckets(stxhistogram) m WHERE stxname = 'stts3';
+</programlisting>
+
+     Values of the <type>pg_histogram</type> can be obtained only from the
+     <literal>pg_statistic.stxhistogram</literal> column.
+   </para>
+  </sect2>
+
   </sect1>
 
 </chapter>
diff --git a/doc/src/sgml/planstats.sgml b/doc/src/sgml/planstats.sgml
index de8ef165c9..67a4f7219c 100644
--- a/doc/src/sgml/planstats.sgml
+++ b/doc/src/sgml/planstats.sgml
@@ -695,6 +695,112 @@ EXPLAIN (ANALYZE, TIMING OFF) SELECT * FROM t WHERE a &lt;= 49 AND b &gt; 49;
 
   </sect2>
 
+  <sect2 id="mv-histograms">
+   <title>Histograms</title>
+
+   <para>
+    <acronym>MCV</acronym> lists, introduced in <xref linkend="mcv-lists"/>,
+    work very well for columns with only a few distinct values, and for
+    columns with only few common values. In those cases, <acronym>MCV</acronym>
+    lists are a very accurate approximation of the real distribution.
+    Histograms, briefly described in <xref linkend="row-estimation-examples"/>,
+    are meant to address the high-cardinality case.
+   </para>
+
+   <para>
+    Although the example data we've used in <xref linkend="mcv-lists"/> does
+    not quality as a high-cardinality case, we can try creating a histogram
+    instead of the <acronym>MCV</acronym> list. With the histogram in place,
+    you may get a plan like this:
+
+<programlisting>
+DROP STATISTICS stts2;
+CREATE STATISTICS stts3 (histogram) ON a, b FROM t;
+ANALYZE t;
+EXPLAIN (ANALYZE, TIMING OFF) SELECT * FROM t WHERE a = 1 AND b = 1;
+                                  QUERY PLAN
+-------------------------------------------------------------------------------
+ Seq Scan on t  (cost=0.00..195.00 rows=100 width=8) (actual rows=100 loops=1)
+   Filter: ((a = 1) AND (b = 1))
+   Rows Removed by Filter: 9900
+</programlisting>
+
+    Which seems quite accurate. For other combinations of values the
+    estimates may be worse, as illustrated by the following query:
+
+<programlisting>
+EXPLAIN (ANALYZE, TIMING OFF) SELECT * FROM t WHERE a = 1 AND b = 10;
+                                 QUERY PLAN                                 
+-----------------------------------------------------------------------------
+ Seq Scan on t  (cost=0.00..195.00 rows=100 width=8) (actual rows=0 loops=1)
+   Filter: ((a = 1) AND (b = 10))
+   Rows Removed by Filter: 10000
+</programlisting>
+
+    This happens due to histograms tracking ranges of values, which makes it
+    impossible to decide how tuples with the exact combination of values there
+    are in the bucket.
+   </para>
+
+   <para>
+    It's also possible to build <acronym>MCV</acronym> lists and a histogram, in
+    which case <command>ANALYZE</command> will build a <acronym>MCV</acronym> lists
+    with the most frequent values, and a histogram on the remaining part of the
+    sample.
+
+<programlisting>
+DROP STATISTICS stts3;
+CREATE STATISTICS stts4 (mcv, histogram) ON a, b FROM t;
+ANALYZE t;
+EXPLAIN (ANALYZE, TIMING OFF) SELECT * FROM t WHERE a = 1 AND b = 10;
+                                QUERY PLAN                                 
+---------------------------------------------------------------------------
+ Seq Scan on t  (cost=0.00..195.00 rows=1 width=8) (actual rows=0 loops=1)
+   Filter: ((a = 1) AND (b = 10))
+   Rows Removed by Filter: 10000
+</programlisting>
+
+    In this case the <acronym>MCV</acronym> list and histogram are treated as a
+    single composed statistics.
+   </para>
+
+   <para>
+    Similarly to <acronym>MCV</acronym> lists, it is possible to inspect
+    histograms using a function called <function>pg_histogram_buckets</function>,
+    which simply lists buckets of a histogram, along with information about
+    boundaries, frequencies, volume etc. When applied to the histogram from
+    <varname>stts3</varname>, you should get something like this:
+
+<programlisting>
+SELECT m.* FROM pg_statistic_ext,
+                pg_histogram_buckets(stxhistogram) m WHERE stxname = 'stts3';
+ index | minvals | maxvals | nullsonly | mininclusive | maxinclusive | frequency | density  | bucket_volume 
+-------+---------+---------+-----------+--------------+--------------+-----------+----------+---------------
+     0 | {0,0}   | {3,1}   | {f,f}     | {t,t}        | {f,f}        |      0.01 |    2.635 |      0.003795
+     1 | {50,0}  | {99,51} | {f,f}     | {t,t}        | {t,f}        |      0.01 | 0.034444 |      0.290323
+     2 | {0,25}  | {26,37} | {f,f}     | {t,t}        | {f,f}        |      0.01 | 0.292778 |      0.034156
+...
+    61 | {60,56} | {99,62} | {f,f}     | {t,t}        | {t,f}        |      0.02 | 0.752857 |      0.026565
+    62 | {35,25} | {50,37} | {f,f}     | {t,t}        | {f,f}        |      0.02 | 0.878333 |       0.02277
+    63 | {81,85} | {87,99} | {f,f}     | {t,t}        | {f,t}        |      0.02 | 1.756667 |      0.011385
+(64 rows)
+</programlisting>
+
+    Which shows that there are 64 buckets, with frequencies ranging between 1%
+    and 2%. The <structfield>minvals</structfield> and <structfield>maxvals</structfield>
+    show the bucket boundaries, <structfield>nullsonly</structfield> shows which
+    columns contain only null values (in the given bucket).
+   </para>
+
+   <para>
+    Similarly to <acronym>MCV</acronym> lists, the planner applies all conditions
+    to the buckets, and sums the frequencies of the matching ones. For details,
+    see <function>histogram_clauselist_selectivity</function> function in
+    <filename>src/backend/statistics/histogram.c</filename>.
+   </para>
+
+  </sect2>
+
  </sect1>
 
  <sect1 id="planner-stats-security">
diff --git a/doc/src/sgml/ref/create_statistics.sgml b/doc/src/sgml/ref/create_statistics.sgml
index fcbfa569d0..ef84341551 100644
--- a/doc/src/sgml/ref/create_statistics.sgml
+++ b/doc/src/sgml/ref/create_statistics.sgml
@@ -84,7 +84,8 @@ CREATE STATISTICS [ IF NOT EXISTS ] <replaceable class="parameter">statistics_na
       <literal>ndistinct</literal>, which enables n-distinct statistics, and
       <literal>dependencies</literal>, which enables functional
       dependency statistics, and <literal>mcv</literal> which enables
-      most-common values lists.
+      most-common values lists, and <literal>histogram</literal> which
+      enables histograms.
       If this clause is omitted, all supported statistics kinds are
       included in the statistics object.
       For more information, see <xref linkend="planner-stats-extended"/>
@@ -190,6 +191,32 @@ EXPLAIN ANALYZE SELECT * FROM t2 WHERE (a = 1) AND (b = 2);
 </programlisting>
   </para>
 
+  <para>
+   Create table <structname>t3</structname> with two strongly correlated
+   columns, and a histogram on those two columns:
+
+<programlisting>
+CREATE TABLE t3 (
+    a   float,
+    b   float
+);
+
+INSERT INTO t3 SELECT mod(i,1000), mod(i,1000) + 50 * (r - 0.5) FROM (
+                   SELECT i, random() r FROM generate_series(1,1000000) s(i)
+                 ) foo;
+
+CREATE STATISTICS s3 WITH (histogram) ON (a, b) FROM t3;
+
+ANALYZE t3;
+
+-- small overlap
+EXPLAIN ANALYZE SELECT * FROM t3 WHERE (a &lt; 500) AND (b &gt; 500);
+
+-- no overlap
+EXPLAIN ANALYZE SELECT * FROM t3 WHERE (a &lt; 400) AND (b &gt; 600);
+</programlisting>
+  </para>
+
  </refsect1>
 
  <refsect1>
diff --git a/src/backend/catalog/system_views.sql b/src/backend/catalog/system_views.sql
index 7251552419..d823e42125 100644
--- a/src/backend/catalog/system_views.sql
+++ b/src/backend/catalog/system_views.sql
@@ -1119,6 +1119,17 @@ LANGUAGE INTERNAL
 STRICT IMMUTABLE PARALLEL SAFE
 AS 'jsonb_insert';
 
+CREATE OR REPLACE FUNCTION
+  pg_histogram_buckets(histogram pg_histogram, otype integer DEFAULT 0,
+            OUT index integer, OUT minvals text[], OUT maxvals text[],
+            OUT nullsonly boolean[], OUT mininclusive boolean[],
+            OUT maxinclusive boolean[], OUT frequency double precision,
+            OUT density double precision, OUT bucket_volume double precision)
+RETURNS SETOF record
+LANGUAGE INTERNAL
+STRICT IMMUTABLE PARALLEL SAFE
+AS 'pg_histogram_buckets';
+
 --
 -- The default permissions for functions mean that anyone can execute them.
 -- A number of functions shouldn't be executable by just anyone, but rather
diff --git a/src/backend/commands/statscmds.c b/src/backend/commands/statscmds.c
index 903d8155e0..d7d504360d 100644
--- a/src/backend/commands/statscmds.c
+++ b/src/backend/commands/statscmds.c
@@ -70,12 +70,13 @@ CreateStatistics(CreateStatsStmt *stmt)
 	Oid			relid;
 	ObjectAddress parentobject,
 				myself;
-	Datum		types[3];		/* one for each possible type of statistic */
+	Datum		types[4];		/* one for each possible type of statistic */
 	int			ntypes;
 	ArrayType  *stxkind;
 	bool		build_ndistinct;
 	bool		build_dependencies;
 	bool		build_mcv;
+	bool		build_histogram;
 	bool		requested_type = false;
 	int			i;
 	ListCell   *cell;
@@ -271,6 +272,7 @@ CreateStatistics(CreateStatsStmt *stmt)
 	build_ndistinct = false;
 	build_dependencies = false;
 	build_mcv = false;
+	build_histogram = false;
 	foreach(cell, stmt->stat_types)
 	{
 		char	   *type = strVal((Value *) lfirst(cell));
@@ -290,6 +292,11 @@ CreateStatistics(CreateStatsStmt *stmt)
 			build_mcv = true;
 			requested_type = true;
 		}
+		else if (strcmp(type, "histogram") == 0)
+		{
+			build_histogram = true;
+			requested_type = true;
+		}
 		else
 			ereport(ERROR,
 					(errcode(ERRCODE_SYNTAX_ERROR),
@@ -302,6 +309,7 @@ CreateStatistics(CreateStatsStmt *stmt)
 		build_ndistinct = true;
 		build_dependencies = true;
 		build_mcv = true;
+		build_histogram = true;
 	}
 
 	/* construct the char array of enabled statistic types */
@@ -312,6 +320,8 @@ CreateStatistics(CreateStatsStmt *stmt)
 		types[ntypes++] = CharGetDatum(STATS_EXT_DEPENDENCIES);
 	if (build_mcv)
 		types[ntypes++] = CharGetDatum(STATS_EXT_MCV);
+	if (build_histogram)
+		types[ntypes++] = CharGetDatum(STATS_EXT_HISTOGRAM);
 	Assert(ntypes > 0 && ntypes <= lengthof(types));
 	stxkind = construct_array(types, ntypes, CHAROID, 1, true, 'c');
 
@@ -331,6 +341,7 @@ CreateStatistics(CreateStatsStmt *stmt)
 	nulls[Anum_pg_statistic_ext_stxndistinct - 1] = true;
 	nulls[Anum_pg_statistic_ext_stxdependencies - 1] = true;
 	nulls[Anum_pg_statistic_ext_stxmcv - 1] = true;
+	nulls[Anum_pg_statistic_ext_stxhistogram - 1] = true;
 
 	/* insert it into pg_statistic_ext */
 	statrel = heap_open(StatisticExtRelationId, RowExclusiveLock);
@@ -435,8 +446,9 @@ RemoveStatisticsById(Oid statsOid)
  * values, this assumption could fail.  But that seems like a corner case
  * that doesn't justify zapping the stats in common cases.)
  *
- * For MCV lists that's not the case, as those statistics store the datums
- * internally. In this case we simply reset the statistics value to NULL.
+ * For MCV lists and histograms that's not the case, as those statistics
+ * store the datums internally. In those cases we simply reset those
+ * statistics to NULL.
  */
 void
 UpdateStatisticsForTypeChange(Oid statsOid, Oid relationOid, int attnum,
@@ -473,9 +485,10 @@ UpdateStatisticsForTypeChange(Oid statsOid, Oid relationOid, int attnum,
 
 	/*
 	 * We can also leave the record as it is if there are no statistics
-	 * including the datum values, like for example MCV lists.
+	 * including the datum values, like for example MCV and histograms.
 	 */
-	if (statext_is_kind_built(oldtup, STATS_EXT_MCV))
+	if (statext_is_kind_built(oldtup, STATS_EXT_MCV) ||
+		statext_is_kind_built(oldtup, STATS_EXT_HISTOGRAM))
 		reset_stats = true;
 
 	/*
@@ -496,11 +509,11 @@ UpdateStatisticsForTypeChange(Oid statsOid, Oid relationOid, int attnum,
 	memset(replaces, 0, Natts_pg_statistic_ext * sizeof(bool));
 	memset(values, 0, Natts_pg_statistic_ext * sizeof(Datum));
 
-	if (statext_is_kind_built(oldtup, STATS_EXT_MCV))
-	{
-		replaces[Anum_pg_statistic_ext_stxmcv - 1] = true;
-		nulls[Anum_pg_statistic_ext_stxmcv - 1] = true;
-	}
+	replaces[Anum_pg_statistic_ext_stxmcv - 1] = true;
+	replaces[Anum_pg_statistic_ext_stxhistogram - 1] = true;
+
+	nulls[Anum_pg_statistic_ext_stxmcv - 1] = true;
+	nulls[Anum_pg_statistic_ext_stxhistogram - 1] = true;
 
 	rel = heap_open(StatisticExtRelationId, RowExclusiveLock);
 
diff --git a/src/backend/nodes/outfuncs.c b/src/backend/nodes/outfuncs.c
index b5af904c18..053cbc498f 100644
--- a/src/backend/nodes/outfuncs.c
+++ b/src/backend/nodes/outfuncs.c
@@ -2437,7 +2437,7 @@ _outStatisticExtInfo(StringInfo str, const StatisticExtInfo *node)
 	/* NB: this isn't a complete set of fields */
 	WRITE_OID_FIELD(statOid);
 	/* don't write rel, leads to infinite recursion in plan tree dump */
-	WRITE_CHAR_FIELD(kind);
+	WRITE_INT_FIELD(kinds);
 	WRITE_BITMAPSET_FIELD(keys);
 }
 
diff --git a/src/backend/optimizer/util/plancat.c b/src/backend/optimizer/util/plancat.c
index 0112450419..1079183ccc 100644
--- a/src/backend/optimizer/util/plancat.c
+++ b/src/backend/optimizer/util/plancat.c
@@ -1324,6 +1324,9 @@ get_relation_statistics(RelOptInfo *rel, Relation relation)
 		HeapTuple	htup;
 		Bitmapset  *keys = NULL;
 		int			i;
+		int			kind = 0;
+
+		StatisticExtInfo *info = makeNode(StatisticExtInfo);
 
 		htup = SearchSysCache1(STATEXTOID, ObjectIdGetDatum(statOid));
 		if (!htup)
@@ -1338,42 +1341,25 @@ get_relation_statistics(RelOptInfo *rel, Relation relation)
 		for (i = 0; i < staForm->stxkeys.dim1; i++)
 			keys = bms_add_member(keys, staForm->stxkeys.values[i]);
 
-		/* add one StatisticExtInfo for each kind built */
+		/* now build the bitmask of statistics kinds */
 		if (statext_is_kind_built(htup, STATS_EXT_NDISTINCT))
-		{
-			StatisticExtInfo *info = makeNode(StatisticExtInfo);
-
-			info->statOid = statOid;
-			info->rel = rel;
-			info->kind = STATS_EXT_NDISTINCT;
-			info->keys = bms_copy(keys);
-
-			stainfos = lcons(info, stainfos);
-		}
+			kind |= STATS_EXT_INFO_NDISTINCT;
 
 		if (statext_is_kind_built(htup, STATS_EXT_DEPENDENCIES))
-		{
-			StatisticExtInfo *info = makeNode(StatisticExtInfo);
-
-			info->statOid = statOid;
-			info->rel = rel;
-			info->kind = STATS_EXT_DEPENDENCIES;
-			info->keys = bms_copy(keys);
-
-			stainfos = lcons(info, stainfos);
-		}
+			kind |= STATS_EXT_INFO_DEPENDENCIES;
 
 		if (statext_is_kind_built(htup, STATS_EXT_MCV))
-		{
-			StatisticExtInfo *info = makeNode(StatisticExtInfo);
+			kind |= STATS_EXT_INFO_MCV;
 
-			info->statOid = statOid;
-			info->rel = rel;
-			info->kind = STATS_EXT_MCV;
-			info->keys = bms_copy(keys);
+		if (statext_is_kind_built(htup, STATS_EXT_HISTOGRAM))
+			kind |= STATS_EXT_INFO_HISTOGRAM;
 
-			stainfos = lcons(info, stainfos);
-		}
+		info->statOid = statOid;
+		info->rel = rel;
+		info->kinds = kind;
+		info->keys = bms_copy(keys);
+
+		stainfos = lcons(info, stainfos);
 
 		ReleaseSysCache(htup);
 		bms_free(keys);
diff --git a/src/backend/parser/parse_utilcmd.c b/src/backend/parser/parse_utilcmd.c
index 29877126d7..58b06fca92 100644
--- a/src/backend/parser/parse_utilcmd.c
+++ b/src/backend/parser/parse_utilcmd.c
@@ -1683,6 +1683,8 @@ generateClonedExtStatsStmt(RangeVar *heapRel, Oid heapRelid,
 			stat_types = lappend(stat_types, makeString("dependencies"));
 		else if (enabled[i] == STATS_EXT_MCV)
 			stat_types = lappend(stat_types, makeString("mcv"));
+		else if (enabled[i] == STATS_EXT_HISTOGRAM)
+			stat_types = lappend(stat_types, makeString("histogram"));
 		else
 			elog(ERROR, "unrecognized statistics kind %c", enabled[i]);
 	}
diff --git a/src/backend/statistics/Makefile b/src/backend/statistics/Makefile
index d2815265fb..3e5ad454cd 100644
--- a/src/backend/statistics/Makefile
+++ b/src/backend/statistics/Makefile
@@ -12,6 +12,6 @@ subdir = src/backend/statistics
 top_builddir = ../../..
 include $(top_builddir)/src/Makefile.global
 
-OBJS = extended_stats.o dependencies.o mcv.o mvdistinct.o
+OBJS = extended_stats.o dependencies.o histogram.o mcv.o mvdistinct.o
 
 include $(top_srcdir)/src/backend/common.mk
diff --git a/src/backend/statistics/README b/src/backend/statistics/README
index 8f153a9e85..9de750614f 100644
--- a/src/backend/statistics/README
+++ b/src/backend/statistics/README
@@ -20,6 +20,8 @@ There are currently two kinds of extended statistics:
 
     (c) MCV lists (README.mcv)
 
+    (d) histograms (README.histogram)
+
 
 Compatible clause types
 -----------------------
@@ -30,6 +32,8 @@ Each type of statistics may be used to estimate some subset of clause types.
 
     (b) MCV lists - equality and inequality clauses (AND, OR, NOT), IS NULL
 
+    (c) histogram - equality and inequality clauses (AND, OR, NOT), IS NULL
+
 Currently, only OpExprs in the form Var op Const, or Const op Var are
 supported, however it's feasible to expand the code later to also estimate the
 selectivities on clauses such as Var op Var.
diff --git a/src/backend/statistics/README.histogram b/src/backend/statistics/README.histogram
new file mode 100644
index 0000000000..e1a4504502
--- /dev/null
+++ b/src/backend/statistics/README.histogram
@@ -0,0 +1,305 @@
+Multivariate histograms
+=======================
+
+Histograms on individual attributes consist of buckets represented by ranges,
+covering the domain of the attribute. That is, each bucket is a [min,max]
+interval, and contains all values in this range. The histogram is built in such
+a way that all buckets have about the same frequency.
+
+Multivariate histograms are an extension into n-dimensional space - the buckets
+are n-dimensional intervals (i.e. n-dimensional rectagles), covering the domain
+of the combination of attributes. That is, each bucket has a vector of lower
+and upper boundaries, denoted min[i] and max[i] (where i = 1..n).
+
+In addition to the boundaries, each bucket tracks additional info:
+
+    * frequency (fraction of tuples in the bucket)
+    * whether the boundaries are inclusive or exclusive
+    * whether the dimension contains only NULL values
+    * number of distinct values in each dimension (for building only)
+
+It's possible that in the future we'll multiple histogram types, with different
+features. We do however expect all the types to share the same representation
+(buckets as ranges) and only differ in how we build them.
+
+The current implementation builds non-overlapping buckets, that may not be true
+for some histogram types and the code should not rely on this assumption. There
+are interesting types of histograms (or algorithms) with overlapping buckets.
+
+When used on low-cardinality data, histograms usually perform considerably worse
+than MCV lists (which are a good fit for this kind of data). This is especially
+true on label-like values, where ordering of the values is mostly unrelated to
+meaning of the data, as proper ordering is crucial for histograms.
+
+On high-cardinality data the histograms are usually a better choice, because MCV
+lists can't represent the distribution accurately enough.
+
+
+Selectivity estimation
+----------------------
+
+The estimation is implemented in clauselist_mv_selectivity_histogram(), and
+works very similarly to clauselist_mv_selectivity_mcvlist.
+
+The main difference is that while MCV lists support exact matches, histograms
+often result in approximate matches - e.g. with equality we can only say if
+the constant would be part of the bucket, but not whether it really is there
+or what fraction of the bucket it corresponds to. In this case we rely on
+some defaults just like in the per-column histograms.
+
+The current implementation uses histograms to estimates those types of clauses
+(think of WHERE conditions):
+
+    (a) equality clauses    WHERE (a = 1) AND (b = 2)
+    (b) inequality clauses  WHERE (a < 1) AND (b >= 2)
+    (c) NULL clauses        WHERE (a IS NULL) AND (b IS NOT NULL)
+    (d) OR-clauses          WHERE (a = 1)  OR (b = 2)
+
+Similarly to MCV lists, it's possible to add support for additional types of
+clauses, for example:
+
+    (e) multi-var clauses   WHERE (a > b)
+
+and so on. These are tasks for the future, not yet implemented.
+
+
+When evaluating a clause on a bucket, we may get one of three results:
+
+    (a) FULL_MATCH - The bucket definitely matches the clause.
+
+    (b) PARTIAL_MATCH - The bucket matches the clause, but not necessarily all
+                        the tuples it represents.
+
+    (c) NO_MATCH - The bucket definitely does not match the clause.
+
+This may be illustrated using a range [1, 5], which is essentially a 1-D bucket.
+With clause
+
+    WHERE (a < 10) => FULL_MATCH (all range values are below
+                      10, so the whole bucket matches)
+
+    WHERE (a < 3)  => PARTIAL_MATCH (there may be values matching
+                      the clause, but we don't know how many)
+
+    WHERE (a < 0)  => NO_MATCH (the whole range is above 1, so
+                      no values from the bucket can match)
+
+Some clauses may produce only some of those results - for example equality
+clauses may never produce FULL_MATCH as we always hit only part of the bucket
+(we can't match both boundaries at the same time). This results in less accurate
+estimates compared to MCV lists, where we can hit a MCV items exactly (there's
+no PARTIAL match in MCV).
+
+There are also clauses that may not produce any PARTIAL_MATCH results. A nice
+example of that is 'IS [NOT] NULL' clause, which either matches the bucket
+completely (FULL_MATCH) or not at all (NO_MATCH), thanks to how the NULL-buckets
+are constructed.
+
+Computing the total selectivity estimate is trivial - simply sum selectivities
+from all the FULL_MATCH and PARTIAL_MATCH buckets (but for buckets marked with
+PARTIAL_MATCH, multiply the frequency by 0.5 to minimize the average error).
+
+
+Building a histogram
+---------------------
+
+The algorithm of building a histogram in general is quite simple:
+
+    (a) create an initial bucket (containing all sample rows)
+
+    (b) create NULL buckets (by splitting the initial bucket)
+
+    (c) repeat
+
+        (1) choose bucket to split next
+
+        (2) terminate if no bucket that might be split found, or if we've
+            reached the maximum number of buckets (16384)
+
+        (3) choose dimension to partition the bucket by
+
+        (4) partition the bucket by the selected dimension
+
+The main complexity is hidden in steps (c.1) and (c.3), i.e. how we choose the
+bucket and dimension for the split, as discussed in the next section.
+
+
+Partitioning criteria
+---------------------
+
+Similarly to one-dimensional histograms, we want to produce buckets with roughly
+the same frequency.
+
+We also need to produce "regular" buckets, because buckets with one dimension
+much longer than the others are very likely to match a lot of conditions (which
+increases error, even if the bucket frequency is very low).
+
+This is especially important when handling OR-clauses, because in that case each
+clause may add buckets independently. With AND-clauses all the clauses have to
+match each bucket, which makes this issue somewhat less concenrning.
+
+To achieve this, we choose the largest bucket (containing the most sample rows),
+but we only choose buckets that can actually be split (have at least 3 different
+combinations of values).
+
+Then we choose the "longest" dimension of the bucket, which is computed by using
+the distinct values in the sample as a measure.
+
+For details see functions select_bucket_to_partition() and partition_bucket(),
+which also includes further discussion.
+
+
+The current limit on number of buckets (16384) is mostly arbitrary, but chosen
+so that it guarantees we don't exceed the number of distinct values indexable by
+uint16 in any of the dimensions. In practice we could handle more buckets as we
+index each dimension separately and the splits should use the dimensions evenly.
+
+Also, histograms this large (with 16k values in multiple dimensions) would be
+quite expensive to build and process, so the 16k limit is rather reasonable.
+
+The actual number of buckets is also related to statistics target, because we
+require MIN_BUCKET_ROWS (10) tuples per bucket before a split, so we can't have
+more than (2 * 300 * target / 10) buckets. For the default target (100) this
+evaluates to ~6k.
+
+
+NULL handling (create_null_buckets)
+-----------------------------------
+
+When building histograms on a single attribute, we first filter out NULL values.
+In the multivariate case, we can't really do that because the rows may contain
+a mix of NULL and non-NULL values in different columns (so we can't simply
+filter all of them out).
+
+For this reason, the histograms are built in a way so that for each bucket, each
+dimension only contains only NULL or non-NULL values. Building the NULL-buckets
+happens as the first step in the build, by the create_null_buckets() function.
+The number of NULL buckets, as produced by this function, has a clear upper
+boundary (2^N) where N is the number of dimensions (attributes the histogram is
+built on). Or rather 2^K where K is the number of attributes that are not marked
+as not-NULL.
+
+The buckets with NULL dimensions are then subject to the same build algorithm
+(i.e. may be split into smaller buckets) just like any other bucket, but may
+only be split by non-NULL dimension.
+
+
+Serialization
+-------------
+
+To store the histogram in pg_statistic_ext table, it is serialized into a more
+efficient form. We also use the representation for estimation, i.e. we don't
+fully deserialize the histogram.
+
+For example the boundary values are deduplicated to minimize the required space.
+How much redundancy is there, actually? Let's assume there are no NULL values,
+so we start with a single bucket - in that case we have 2*N boundaries. Each
+time we split a bucket we introduce one new value (in the "middle" of one of
+the dimensions), and keep boundries for all the other dimensions. So after K
+splits, we have up to
+
+    2*N + K
+
+unique boundary values (we may have fewe values, if the same value is used for
+several splits). But after K splits we do have (K+1) buckets, so
+
+    (K+1) * 2 * N
+
+boundary values. Using e.g. N=4 and K=999, we arrive to those numbers:
+
+    2*N + K       = 1007
+    (K+1) * 2 * N = 8000
+
+wich means a lot of redundancy. It's somewhat counter-intuitive that the number
+of distinct values does not really depend on the number of dimensions (except
+for the initial bucket, but that's negligible compared to the total).
+
+By deduplicating the values and replacing them with 16-bit indexes (uint16), we
+reduce the required space to
+
+    1007 * 8 + 8000 * 2 ~= 24kB
+
+which is significantly less than 64kB required for the 'raw' histogram (assuming
+the values are 8B).
+
+While the bytea compression (pglz) might achieve the same reduction of space,
+the deduplicated representation is used to optimize the estimation by caching
+results of function calls for already visited values. This significantly
+reduces the number of calls to (often quite expensive) operators.
+
+Note: Of course, this reasoning only holds for histograms built by the algorithm
+that simply splits the buckets in half. Other histograms types (e.g. containing
+overlapping buckets) may behave differently and require different serialization.
+
+Serialized histograms are marked with 'magic' constant, to make it easier to
+check the bytea value really is a serialized histogram.
+
+
+varlena compression
+-------------------
+
+This serialization may however disable automatic varlena compression, the array
+of unique values is placed at the beginning of the serialized form. Which is
+exactly the chunk used by pglz to check if the data is compressible, and it
+will probably decide it's not very compressible. This is similar to the issue
+we had with JSONB initially.
+
+Maybe storing buckets first would make it work, as the buckets may be better
+compressible.
+
+On the other hand the serialization is actually a context-aware compression,
+usually compressing to ~30% (or even less, with large data types). So the lack
+of additional pglz compression may be acceptable.
+
+
+Deserialization
+---------------
+
+The deserialization is not a perfect inverse of the serialization, as we keep
+the deduplicated arrays. This reduces the amount of memory and also allows
+optimizations during estimation (e.g. we can cache results for the distinct
+values, saving expensive function calls).
+
+
+Inspecting the histogram
+------------------------
+
+Inspecting the regular (per-attribute) histograms is trivial, as it's enough
+to select the columns from pg_stats. The data is encoded as anyarrays, and
+all the items have the same data type, so anyarray provides a simple way to
+get a text representation.
+
+With multivariate histograms the columns may use different data types, making
+it impossible to use anyarrays. It might be possible to produce similar
+array-like representation, but that would complicate further processing and
+analysis of the histogram.
+
+So instead the histograms are stored in a custom data type (pg_histogram),
+which however makes it more difficult to inspect the contents. To make that
+easier, there's a SRF returning detailed information about the histogram.
+
+    SELECT * FROM pg_histogram_buckets();
+
+It has two input parameters:
+
+    histogram - OID of the histogram (pg_statistic_ext.stxhistogram)
+    otype     - type of output
+
+and produces a table with these columns:
+
+    - bucket ID                (0...nbuckets-1)
+    - lower bucket boundaries  (string array)
+    - upper bucket boundaries  (string array)
+    - nulls only dimensions    (boolean array)
+    - lower boundary inclusive (boolean array)
+    - upper boundary includive (boolean array)
+    - frequency                (double precision)
+    - density                  (double precision)
+    - volume                   (double precision)
+
+The 'otype' accepts three values, determining what will be returned in the
+lower/upper boundary arrays:
+
+    - 0 - values stored in the histogram, encoded as text
+    - 1 - indexes into the deduplicated arrays
+    - 2 - idnexes into the deduplicated arrays, scaled to [0,1]
diff --git a/src/backend/statistics/dependencies.c b/src/backend/statistics/dependencies.c
index 29e816c4f7..5505feb913 100644
--- a/src/backend/statistics/dependencies.c
+++ b/src/backend/statistics/dependencies.c
@@ -932,7 +932,7 @@ dependencies_clauselist_selectivity(PlannerInfo *root,
 	int			listidx;
 
 	/* check if there's any stats that might be useful for us. */
-	if (!has_stats_of_kind(rel->statlist, STATS_EXT_DEPENDENCIES))
+	if (!has_stats_of_kind(rel->statlist, STATS_EXT_INFO_DEPENDENCIES))
 		return 1.0;
 
 	list_attnums = (AttrNumber *) palloc(sizeof(AttrNumber) *
diff --git a/src/backend/statistics/extended_stats.c b/src/backend/statistics/extended_stats.c
index 0b66000705..25237cd53c 100644
--- a/src/backend/statistics/extended_stats.c
+++ b/src/backend/statistics/extended_stats.c
@@ -38,7 +38,6 @@
 #include "utils/selfuncs.h"
 #include "utils/syscache.h"
 
-
 /*
  * Used internally to refer to an individual statistics object, i.e.,
  * a pg_statistic_ext entry.
@@ -58,7 +57,7 @@ static VacAttrStats **lookup_var_attr_stats(Relation rel, Bitmapset *attrs,
 					  int nvacatts, VacAttrStats **vacatts);
 static void statext_store(Relation pg_stext, Oid relid,
 			  MVNDistinct *ndistinct, MVDependencies *dependencies,
-			  MCVList * mcvlist, VacAttrStats **stats);
+			  MCVList * mcvlist, MVHistogram * histogram, VacAttrStats **stats);
 
 
 /*
@@ -92,10 +91,14 @@ BuildRelationExtStatistics(Relation onerel, double totalrows,
 		StatExtEntry *stat = (StatExtEntry *) lfirst(lc);
 		MVNDistinct *ndistinct = NULL;
 		MVDependencies *dependencies = NULL;
+		MVHistogram *histogram = NULL;
 		MCVList    *mcv = NULL;
 		VacAttrStats **stats;
 		ListCell   *lc2;
 
+		bool		build_mcv = false;
+		bool		build_histogram = false;
+
 		/*
 		 * Check if we can build these stats based on the column analyzed. If
 		 * not, report this fact (except in autovacuum) and move on.
@@ -131,12 +134,49 @@ BuildRelationExtStatistics(Relation onerel, double totalrows,
 				dependencies = statext_dependencies_build(numrows, rows,
 														  stat->columns, stats);
 			else if (t == STATS_EXT_MCV)
-				mcv = statext_mcv_build(numrows, rows, stat->columns, stats,
-										totalrows);
+				build_mcv = true;
+			else if (t == STATS_EXT_HISTOGRAM)
+				build_histogram = true;
+		}
+
+		/*
+		 * If asked to build both MCV and histogram, first build the MCV part
+		 * and then histogram on the remaining rows.
+		 */
+		if (build_mcv && build_histogram)
+		{
+			HeapTuple  *rows_filtered = NULL;
+			int			numrows_filtered;
+
+			mcv = statext_mcv_build(numrows, rows, stat->columns, stats,
+									&rows_filtered, &numrows_filtered,
+									totalrows);
+
+			/*
+			 * Only build the histogram when there are rows not covered by
+			 * MCV.
+			 */
+			if (rows_filtered)
+			{
+				Assert(numrows_filtered > 0);
+
+				histogram = statext_histogram_build(numrows_filtered, rows_filtered,
+													stat->columns, stats, numrows);
+
+				/* free this immediately, as we may be building many stats */
+				pfree(rows_filtered);
+			}
 		}
+		else if (build_mcv)
+			mcv = statext_mcv_build(numrows, rows, stat->columns, stats,
+									NULL, NULL, totalrows);
+		else if (build_histogram)
+			histogram = statext_histogram_build(numrows, rows, stat->columns,
+												stats, numrows);
 
 		/* store the statistics in the catalog */
-		statext_store(pg_stext, stat->statOid, ndistinct, dependencies, mcv, stats);
+		statext_store(pg_stext, stat->statOid, ndistinct, dependencies, mcv,
+					  histogram, stats);
 	}
 
 	heap_close(pg_stext, RowExclusiveLock);
@@ -168,6 +208,10 @@ statext_is_kind_built(HeapTuple htup, char type)
 			attnum = Anum_pg_statistic_ext_stxmcv;
 			break;
 
+		case STATS_EXT_HISTOGRAM:
+			attnum = Anum_pg_statistic_ext_stxhistogram;
+			break;
+
 		default:
 			elog(ERROR, "unexpected statistics type requested: %d", type);
 	}
@@ -233,7 +277,8 @@ fetch_statentries_for_relation(Relation pg_statext, Oid relid)
 		{
 			Assert((enabled[i] == STATS_EXT_NDISTINCT) ||
 				   (enabled[i] == STATS_EXT_DEPENDENCIES) ||
-				   (enabled[i] == STATS_EXT_MCV));
+				   (enabled[i] == STATS_EXT_MCV) ||
+				   (enabled[i] == STATS_EXT_HISTOGRAM));
 			entry->types = lappend_int(entry->types, (int) enabled[i]);
 		}
 
@@ -308,7 +353,7 @@ lookup_var_attr_stats(Relation rel, Bitmapset *attrs,
 static void
 statext_store(Relation pg_stext, Oid statOid,
 			  MVNDistinct *ndistinct, MVDependencies *dependencies,
-			  MCVList * mcv, VacAttrStats **stats)
+			  MCVList * mcv, MVHistogram * histogram, VacAttrStats **stats)
 {
 	HeapTuple	stup,
 				oldtup;
@@ -347,10 +392,18 @@ statext_store(Relation pg_stext, Oid statOid,
 		values[Anum_pg_statistic_ext_stxmcv - 1] = PointerGetDatum(data);
 	}
 
+	if (histogram != NULL)
+	{
+		/* histogram already is a bytea value, not need to serialize */
+		nulls[Anum_pg_statistic_ext_stxhistogram - 1] = (histogram == NULL);
+		values[Anum_pg_statistic_ext_stxhistogram - 1] = PointerGetDatum(histogram);
+	}
+
 	/* always replace the value (either by bytea or NULL) */
 	replaces[Anum_pg_statistic_ext_stxndistinct - 1] = true;
 	replaces[Anum_pg_statistic_ext_stxdependencies - 1] = true;
 	replaces[Anum_pg_statistic_ext_stxmcv - 1] = true;
+	replaces[Anum_pg_statistic_ext_stxhistogram - 1] = true;
 
 	/* there should already be a pg_statistic_ext tuple */
 	oldtup = SearchSysCache1(STATEXTOID, ObjectIdGetDatum(statOid));
@@ -465,6 +518,19 @@ compare_scalars_simple(const void *a, const void *b, void *arg)
 								 (SortSupport) arg);
 }
 
+/*
+ * qsort_arg comparator for sorting data when partitioning a MV bucket
+ */
+int
+compare_scalars_partition(const void *a, const void *b, void *arg)
+{
+	Datum		da = ((ScalarItem *) a)->value;
+	Datum		db = ((ScalarItem *) b)->value;
+	SortSupport ssup = (SortSupport) arg;
+
+	return ApplySortComparator(da, false, db, false, ssup);
+}
+
 int
 compare_datums_simple(Datum a, Datum b, SortSupport ssup)
 {
@@ -590,10 +656,11 @@ build_sorted_items(int numrows, HeapTuple *rows, TupleDesc tdesc,
 
 /*
  * has_stats_of_kind
- *		Check whether the list contains statistic of a given kind
+ *		Check whether the list contains statistic of a given kind (at least
+ * one of those specified statistics types).
  */
 bool
-has_stats_of_kind(List *stats, char requiredkind)
+has_stats_of_kind(List *stats, int requiredkinds)
 {
 	ListCell   *l;
 
@@ -601,7 +668,7 @@ has_stats_of_kind(List *stats, char requiredkind)
 	{
 		StatisticExtInfo *stat = (StatisticExtInfo *) lfirst(l);
 
-		if (stat->kind == requiredkind)
+		if (stat->kinds & requiredkinds)
 			return true;
 	}
 
@@ -623,7 +690,7 @@ has_stats_of_kind(List *stats, char requiredkind)
  * further tiebreakers are needed.
  */
 StatisticExtInfo *
-choose_best_statistics(List *stats, Bitmapset *attnums, char requiredkind)
+choose_best_statistics(List *stats, Bitmapset *attnums, int requiredkinds)
 {
 	ListCell   *lc;
 	StatisticExtInfo *best_match = NULL;
@@ -637,8 +704,8 @@ choose_best_statistics(List *stats, Bitmapset *attnums, char requiredkind)
 		int			numkeys;
 		Bitmapset  *matched;
 
-		/* skip statistics that are not of the correct type */
-		if (info->kind != requiredkind)
+		/* skip statistics that do not match any of the requested types */
+		if ((info->kinds & requiredkinds) == 0)
 			continue;
 
 		/* determine how many attributes of these stats can be matched to */
@@ -843,7 +910,7 @@ statext_is_compatible_clause_internal(Node *clause, Index relid, Bitmapset **att
 
 /*
  * statext_is_compatible_clause
- *		Determines if the clause is compatible with MCV lists.
+ *		Determines if the clause is compatible with MCV lists and histograms
  *
  * Only OpExprs with two arguments using an equality operator are supported.
  * When returning True attnum is set to the attribute number of the Var within
@@ -873,6 +940,89 @@ statext_is_compatible_clause(Node *clause, Index relid, Bitmapset **attnums)
 }
 
 /*
+ * examine_equality_clause
+ *		Extract variable from a simple top-level equality clause.
+ *
+ * For simple equality clause (Var = Const) or (Const = Var) extracts
+ * the Var. For other clauses returns NULL.
+ */
+static Var *
+examine_equality_clause(PlannerInfo *root, RestrictInfo *rinfo)
+{
+	OpExpr	   *expr;
+	Var		   *var;
+	bool		ok;
+	bool		varonleft = true;
+
+	if (!IsA(rinfo->clause, OpExpr))
+		return NULL;
+
+	expr = (OpExpr *) rinfo->clause;
+
+	if (list_length(expr->args) != 2)
+		return NULL;
+
+	/* see if it actually has the right */
+	ok = (NumRelids((Node *) expr) == 1) &&
+		(is_pseudo_constant_clause(lsecond(expr->args)) ||
+		 (varonleft = false,
+		  is_pseudo_constant_clause(linitial(expr->args))));
+
+	/* unsupported structure (two variables or so) */
+	if (!ok)
+		return NULL;
+
+	if (get_oprrest(expr->opno) != F_EQSEL)
+		return NULL;
+
+	var = (varonleft) ? linitial(expr->args) : lsecond(expr->args);
+
+	return var;
+}
+
+/*
+ * estimate_equality_groups
+ *		Estimates number of groups for attributes in equality clauses.
+ *
+ * Extracts simple top-level equality clauses, and estimates ndistinct
+ * for that combination (using simplified estimate_num_groups). Then
+ * returns number of attributes with an equality clause, and a lists
+ * of equality clauses (to use as conditions for histograms) and also
+ * remaining non-equality clauses.
+ */
+static double
+estimate_equality_groups(PlannerInfo *root, List *clauses,
+						 List **eqclauses, List **neqclauses)
+{
+	List   *vars = NIL;
+	ListCell *lc;
+
+	*eqclauses = NIL;
+	*neqclauses = NIL;
+
+	foreach(lc, clauses)
+	{
+		Var	   *var;
+		RestrictInfo *rinfo = (RestrictInfo *) lfirst(lc);
+
+		Assert(IsA(rinfo, RestrictInfo));
+
+		var = examine_equality_clause(root, rinfo);
+
+		/* is it a simple equality clause */
+		if (var)
+		{
+			vars = lappend(vars, var);
+			*eqclauses = lappend(*eqclauses, rinfo);
+		}
+		else
+			*neqclauses = lappend(*neqclauses, rinfo);
+	}
+
+	return estimate_num_groups_simple(root, vars);
+}
+
+/*
  * statext_clauselist_selectivity
  *		Estimate clauses using the best multi-column statistics.
  *
@@ -937,13 +1087,14 @@ statext_clauselist_selectivity(PlannerInfo *root, List *clauses, int varRelid,
 				mcv_sel,
 				mcv_basesel,
 				mcv_totalsel,
+				histogram_sel,
 				other_sel,
 				sel;
 
-	/* we're interested in MCV lists */
-	int			types = STATS_EXT_MCV;
+	/* we're interested in MCV lists and histograms */
+	int			types = (STATS_EXT_INFO_MCV | STATS_EXT_INFO_HISTOGRAM);
 
-	/* check if there's any stats that might be useful for us. */
+	/* Check if there's any stats that might be useful for us. */
 	if (!has_stats_of_kind(rel->statlist, types))
 		return (Selectivity) 1.0;
 
@@ -994,8 +1145,8 @@ statext_clauselist_selectivity(PlannerInfo *root, List *clauses, int varRelid,
 	if (!stat)
 		return (Selectivity) 1.0;
 
-	/* We only understand MCV lists for now. */
-	Assert(stat->kind == STATS_EXT_MCV);
+	/* We only understand MCV lists and histograms for now. */
+	Assert(stat->kinds & (STATS_EXT_INFO_MCV | STATS_EXT_INFO_HISTOGRAM));
 
 	/* now filter the clauses to be estimated using the selected MCV */
 	stat_clauses = NIL;
@@ -1018,28 +1169,59 @@ statext_clauselist_selectivity(PlannerInfo *root, List *clauses, int varRelid,
 	}
 
 	/*
-	 * First compute "simple" selectivity, i.e. without the extended statistics,
-	 * and essentially assuming independence of the columns/clauses. We'll then
-	 * use the various selectivities computed from MCV list to improve it.
+	 * For statistics with MCV list, we'll estimate the MCV and non-MCV parts.
 	 */
-	simple_sel = clauselist_selectivity_simple(root, stat_clauses, varRelid,
-											   jointype, sjinfo, NULL);
+	if (stat->kinds & STATS_EXT_INFO_MCV)
+	{
+		/*
+		 * First compute "simple" selectivity, i.e. without the extended statistics,
+		 * and essentially assuming independence of the columns/clauses. We'll then
+		 * use the various selectivities computed from MCV list to improve it.
+		 */
+		simple_sel = clauselist_selectivity_simple(root, stat_clauses, varRelid,
+												   jointype, sjinfo, NULL);
+
+		/*
+		 * Now compute the multi-column estimate from the MCV list, along with the
+		 * other selectivities (base & total selectivity).
+		 */
+		mcv_sel = mcv_clauselist_selectivity(root, stat, stat_clauses, varRelid,
+											 jointype, sjinfo, rel,
+											 &mcv_basesel, &mcv_totalsel);
+
+		/* Estimated selectivity of values not covered by MCV matches */
+		other_sel = simple_sel - mcv_basesel;
+		CLAMP_PROBABILITY(other_sel);
+
+		/* The non-MCV selectivity can't exceed the 1 - mcv_totalsel. */
+		if (other_sel > 1.0 - mcv_totalsel)
+			other_sel = 1.0 - mcv_totalsel;
+	}
+	else
+	{
+		/* Otherwise just remember there was no MCV list. */
+		mcv_totalsel = 0.0;
+	}
 
 	/*
-	 * Now compute the multi-column estimate from the MCV list, along with the
-	 * other selectivities (base & total selectivity).
+	 * If we have a histogram, we'll use it to improve the non-MCV estimate.
 	 */
-	mcv_sel = mcv_clauselist_selectivity(root, stat, stat_clauses, varRelid,
-										 jointype, sjinfo, rel,
-										 &mcv_basesel, &mcv_totalsel);
+	if (stat->kinds & STATS_EXT_INFO_HISTOGRAM)
+	{
+		List   *eqclauses,
+			   *neqclauses;
+		double	ngroups;
 
-	/* Estimated selectivity of values not covered by MCV matches */
-	other_sel = simple_sel - mcv_basesel;
-	CLAMP_PROBABILITY(other_sel);
+		ngroups = estimate_equality_groups(root, stat_clauses,
+										   &eqclauses, &neqclauses);
 
-	/* The non-MCV selectivity can't exceed the 1 - mcv_totalsel. */
-	if (other_sel > 1.0 - mcv_totalsel)
-		other_sel = 1.0 - mcv_totalsel;
+		histogram_sel = histogram_clauselist_selectivity(root, stat,
+														 neqclauses, eqclauses,
+														 varRelid, jointype,
+														 sjinfo, rel);
+
+		other_sel = (1 / ngroups) * histogram_sel;
+	}
 
 	/* Overall selectivity is the combination of MCV and non-MCV estimates. */
 	sel = mcv_sel + other_sel;
diff --git a/src/backend/statistics/histogram.c b/src/backend/statistics/histogram.c
new file mode 100644
index 0000000000..1ff34a53c0
--- /dev/null
+++ b/src/backend/statistics/histogram.c
@@ -0,0 +1,3019 @@
+/*-------------------------------------------------------------------------
+ *
+ * histogram.c
+ *	  POSTGRES multivariate histograms
+ *
+ * Portions Copyright (c) 1996-2017, PostgreSQL Global Development Group
+ * Portions Copyright (c) 1994, Regents of the University of California
+ *
+ * IDENTIFICATION
+ *	  src/backend/statistics/histogram.c
+ *-------------------------------------------------------------------------
+ */
+#include "postgres.h"
+
+#include <math.h>
+
+#include "access/htup_details.h"
+#include "catalog/pg_collation.h"
+#include "catalog/pg_statistic_ext.h"
+#include "fmgr.h"
+#include "funcapi.h"
+#include "optimizer/clauses.h"
+#include "statistics/extended_stats_internal.h"
+#include "statistics/statistics.h"
+#include "utils/builtins.h"
+#include "utils/bytea.h"
+#include "utils/fmgroids.h"
+#include "utils/lsyscache.h"
+#include "utils/selfuncs.h"
+#include "utils/syscache.h"
+#include "utils/typcache.h"
+
+
+/*
+ * Multivariate histograms
+ */
+typedef struct MVBucketBuild
+{
+	/* Frequencies of this bucket. */
+	float		frequency;
+
+	/*
+	 * Information about dimensions being NULL-only. Not yet used.
+	 */
+	bool	   *nullsonly;
+
+	/* lower boundaries - values and information about the inequalities */
+	Datum	   *min;
+	bool	   *min_inclusive;
+
+	/* upper boundaries - values and information about the inequalities */
+	Datum	   *max;
+	bool	   *max_inclusive;
+
+	/* number of distinct values in each dimension */
+	uint32	   *ndistincts;
+
+	/* number of distinct combination of values */
+	uint32		ndistinct;
+
+	/* aray of sample rows (for this bucket) */
+	HeapTuple  *rows;
+	uint32		numrows;
+
+}			MVBucketBuild;
+
+typedef struct MVHistogramBuild
+{
+	int32		vl_len_;		/* unused: ensure same alignment as
+								 * MVHistogram for serialization */
+	uint32		magic;			/* magic constant marker */
+	uint32		type;			/* type of histogram (BASIC) */
+	uint32		nbuckets;		/* number of buckets (buckets array) */
+	uint32		ndimensions;	/* number of dimensions */
+	Oid			types[STATS_MAX_DIMENSIONS];	/* OIDs of data types */
+	MVBucketBuild **buckets;	/* array of buckets */
+}			MVHistogramBuild;
+
+static MVBucketBuild * create_initial_ext_bucket(int numrows, HeapTuple *rows,
+												 Bitmapset *attrs, VacAttrStats **stats);
+
+static MVBucketBuild * select_bucket_to_partition(int nbuckets, MVBucketBuild * *buckets);
+
+static MVBucketBuild * partition_bucket(MVBucketBuild * bucket, Bitmapset *attrs,
+										VacAttrStats **stats,
+										int *ndistvalues, Datum **distvalues);
+
+static MVBucketBuild * copy_ext_bucket(MVBucketBuild * bucket, uint32 ndimensions);
+
+static void update_bucket_ndistinct(MVBucketBuild * bucket, Bitmapset *attrs,
+						VacAttrStats **stats);
+
+static void update_dimension_ndistinct(MVBucketBuild * bucket, int dimension,
+						   Bitmapset *attrs, VacAttrStats **stats,
+						   bool update_boundaries);
+
+static void create_null_buckets(MVHistogramBuild * histogram, int bucket_idx,
+					Bitmapset *attrs, VacAttrStats **stats);
+
+static Datum *build_ndistinct(int numrows, HeapTuple *rows, Bitmapset *attrs,
+				VacAttrStats **stats, int i, int *nvals);
+
+static MVHistogram * serialize_histogram(MVHistogramBuild * histogram,
+										 VacAttrStats **stats);
+
+/*
+ * Computes size of a serialized histogram bucket, depending on the number
+ * of dimentions (columns) the statistic is defined on. The datum values
+ * are stored in a separate array (deduplicated, to minimize the size), and
+ * so the serialized buckets only store uint16 indexes into that array.
+ *
+ * Each serialized bucket needs to store (in this order):
+ *
+ * - number of tuples     (float)
+ * - number of distinct   (float)
+ * - min inclusive flags  (ndim * sizeof(bool))
+ * - max inclusive flags  (ndim * sizeof(bool))
+ * - null dimension flags (ndim * sizeof(bool))
+ * - min boundary indexes (2 * ndim * sizeof(uint16))
+ * - max boundary indexes (2 * ndim * sizeof(uint16))
+ *
+ * So in total:
+ *
+ *	 ndim * (4 * sizeof(uint16) + 3 * sizeof(bool)) + (2 * sizeof(float))
+ *
+ * XXX We might save a bit more space by using proper bitmaps instead of
+ * boolean arrays.
+ */
+#define BUCKET_SIZE(ndims)	\
+	(ndims * (4 * sizeof(uint16) + 3 * sizeof(bool)) + sizeof(float))
+
+/*
+ * Macros for convenient access to parts of a serialized bucket.
+ */
+#define BUCKET_FREQUENCY(b)		(*(float*)b)
+#define BUCKET_MIN_INCL(b,n)	((bool*)(b + sizeof(float)))
+#define BUCKET_MAX_INCL(b,n)	(BUCKET_MIN_INCL(b,n) + n)
+#define BUCKET_NULLS_ONLY(b,n)	(BUCKET_MAX_INCL(b,n) + n)
+#define BUCKET_MIN_INDEXES(b,n) ((uint16*)(BUCKET_NULLS_ONLY(b,n) + n))
+#define BUCKET_MAX_INDEXES(b,n) ((BUCKET_MIN_INDEXES(b,n) + n))
+
+/*
+ * Minimal number of rows per bucket (can't split smaller buckets).
+ *
+ * XXX The single-column statistics (std_typanalyze) pretty much says we
+ * need 300 rows per bucket. Should we use the same value here?
+ */
+#define MIN_BUCKET_ROWS			10
+
+/*
+ * Represents match info for a histogram bucket.
+ */
+typedef struct bucket_match
+{
+	bool		match;		/* true/false */
+	double		fraction;	/* fraction of bucket */
+} bucket_match;
+
+/*
+ * Builds a multivariate histogram from the set of sampled rows.
+ *
+ * The build algorithm is iterative - initially a single bucket containing all
+ * sample rows is formed, and then repeatedly split into smaller buckets. In
+ * each round the largest bucket is split into two smaller ones.
+ *
+ * The criteria for selecting the largest bucket (and the dimension for the
+ * split) needs to be elaborate enough to produce buckets of roughly the same
+ * size, and also regular shape (not very narrow in just one dimension).
+ *
+ * The current algorithm works like this:
+ *
+ *   a) build NULL-buckets (create_null_buckets)
+ *
+ *   b) while [maximum number of buckets not reached]
+ *
+ *   c) choose bucket to partition (largest bucket)
+ *
+ *       c.1) if no bucket eligible to split, terminate the build
+ *
+ *       c.2) choose bucket dimension to partition (largest dimension)
+ *
+ *       c.3) split the bucket into two buckets
+ *
+ * See the discussion at select_bucket_to_partition and partition_bucket for
+ * more details about the algorithm.
+ *
+ * The function does not update the interan pointers, hence the histogram
+ * is suitable only for storing. Before using it for estimation, it needs
+ * to go through statext_histogram_deserialize() first.
+ */
+MVHistogram *
+statext_histogram_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
+						VacAttrStats **stats, int numrows_total)
+{
+	int			i;
+	int			numattrs = bms_num_members(attrs);
+
+	int		   *ndistvalues;
+	Datum	  **distvalues;
+
+	MVHistogramBuild *histogram;
+	HeapTuple  *rows_copy;
+
+	/* not supposed to build of too few or too many columns */
+	Assert((numattrs >= 2) && (numattrs <= STATS_MAX_DIMENSIONS));
+
+	/* we need to make a copy of the row array, as we'll modify it */
+	rows_copy = (HeapTuple *) palloc0(numrows * sizeof(HeapTuple));
+	memcpy(rows_copy, rows, sizeof(HeapTuple) * numrows);
+
+	/* build the histogram header */
+
+	histogram = (MVHistogramBuild *) palloc0(sizeof(MVHistogramBuild));
+
+	histogram->magic = STATS_HIST_MAGIC;
+	histogram->type = STATS_HIST_TYPE_BASIC;
+	histogram->ndimensions = numattrs;
+	histogram->nbuckets = 1;	/* initially just a single bucket */
+
+	/*
+	 * Allocate space for maximum number of buckets (better than repeatedly
+	 * doing repalloc for short-lived objects).
+	 */
+	histogram->buckets
+		= (MVBucketBuild * *) palloc0(STATS_HIST_MAX_BUCKETS * sizeof(MVBucketBuild));
+
+	/* Create the initial bucket, covering all sampled rows */
+	histogram->buckets[0]
+		= create_initial_ext_bucket(numrows, rows_copy, attrs, stats);
+
+	/*
+	 * Collect info on distinct values in each dimension (used later to pick
+	 * dimension to partition).
+	 */
+	ndistvalues = (int *) palloc0(sizeof(int) * numattrs);
+	distvalues = (Datum **) palloc0(sizeof(Datum *) * numattrs);
+
+	for (i = 0; i < numattrs; i++)
+		distvalues[i] = build_ndistinct(numrows, rows, attrs, stats, i,
+										&ndistvalues[i]);
+
+	/*
+	 * Split the initial bucket into buckets that don't mix NULL and non-NULL
+	 * values in a single dimension.
+	 *
+	 * XXX Maybe this should be happening before the build_ndistinct()?
+	 */
+	create_null_buckets(histogram, 0, attrs, stats);
+
+	/*
+	 * Split the buckets into smaller and smaller buckets. The loop will end
+	 * when either all buckets are too small (MIN_BUCKET_ROWS), or there are
+	 * too many buckets in total (STATS_HIST_MAX_BUCKETS).
+	 */
+	while (histogram->nbuckets < STATS_HIST_MAX_BUCKETS)
+	{
+		MVBucketBuild *bucket = select_bucket_to_partition(histogram->nbuckets,
+														   histogram->buckets);
+
+		/* no bucket eligible for partitioning */
+		if (bucket == NULL)
+			break;
+
+		/* we modify the bucket in-place and add one new bucket */
+		histogram->buckets[histogram->nbuckets++]
+			= partition_bucket(bucket, attrs, stats, ndistvalues, distvalues);
+	}
+
+	/* Finalize the histogram build - compute bucket frequencies etc. */
+	for (i = 0; i < histogram->nbuckets; i++)
+	{
+		/*
+		 * The frequency has to be computed from the whole sample, in case
+		 * some of the rows were filtered out in the MCV build.
+		 */
+		histogram->buckets[i]->frequency
+			= (histogram->buckets[i]->numrows * 1.0) / numrows_total;
+	}
+
+	return serialize_histogram(histogram, stats);
+}
+
+/*
+ * build_ndistinct
+ *		build array of ndistinct values in a particular column, count them
+ *
+ */
+static Datum *
+build_ndistinct(int numrows, HeapTuple *rows, Bitmapset *attrs,
+				VacAttrStats **stats, int i, int *nvals)
+{
+	int			j;
+	int			nvalues,
+				ndistinct;
+	Datum	   *values,
+			   *distvalues;
+	int		   *attnums;
+
+	TypeCacheEntry *type;
+	SortSupportData ssup;
+
+	type = lookup_type_cache(stats[i]->attrtypid, TYPECACHE_LT_OPR);
+
+	/* initialize sort support, etc. */
+	memset(&ssup, 0, sizeof(ssup));
+	ssup.ssup_cxt = CurrentMemoryContext;
+
+	/* We always use the default collation for statistics */
+	ssup.ssup_collation = DEFAULT_COLLATION_OID;
+	ssup.ssup_nulls_first = false;
+
+	PrepareSortSupportFromOrderingOp(type->lt_opr, &ssup);
+
+	nvalues = 0;
+	values = (Datum *) palloc0(sizeof(Datum) * numrows);
+
+	attnums = build_attnums(attrs);
+
+	/* collect values from the sample rows, ignore NULLs */
+	for (j = 0; j < numrows; j++)
+	{
+		Datum		value;
+		bool		isnull;
+
+		/*
+		 * remember the index of the sample row, to make the partitioning
+		 * simpler
+		 */
+		value = heap_getattr(rows[j], attnums[i],
+							 stats[i]->tupDesc, &isnull);
+
+		if (isnull)
+			continue;
+
+		values[nvalues++] = value;
+	}
+
+	/* if no non-NULL values were found, free the memory and terminate */
+	if (nvalues == 0)
+	{
+		pfree(values);
+		return NULL;
+	}
+
+	/* sort the array of values using the SortSupport */
+	qsort_arg((void *) values, nvalues, sizeof(Datum),
+			  compare_scalars_simple, (void *) &ssup);
+
+	/* count the distinct values first, and allocate just enough memory */
+	ndistinct = 1;
+	for (j = 1; j < nvalues; j++)
+		if (compare_scalars_simple(&values[j], &values[j - 1], &ssup) != 0)
+			ndistinct += 1;
+
+	distvalues = (Datum *) palloc0(sizeof(Datum) * ndistinct);
+
+	/* now collect distinct values into the array */
+	distvalues[0] = values[0];
+	ndistinct = 1;
+
+	for (j = 1; j < nvalues; j++)
+	{
+		if (compare_scalars_simple(&values[j], &values[j - 1], &ssup) != 0)
+		{
+			distvalues[ndistinct] = values[j];
+			ndistinct += 1;
+		}
+	}
+
+	pfree(values);
+
+	*nvals = ndistinct;
+	return distvalues;
+}
+
+/*
+ * statext_histogram_load
+ *		Load the histogram list for the indicated pg_statistic_ext tuple
+*/
+MVHistogram *
+statext_histogram_load(Oid mvoid)
+{
+	bool		isnull = false;
+	Datum		histogram;
+	HeapTuple	htup = SearchSysCache1(STATEXTOID, ObjectIdGetDatum(mvoid));
+
+	if (!HeapTupleIsValid(htup))
+		elog(ERROR, "cache lookup failed for statistics object %u", mvoid);
+
+	histogram = SysCacheGetAttr(STATEXTOID, htup,
+								Anum_pg_statistic_ext_stxhistogram, &isnull);
+
+	ReleaseSysCache(htup);
+
+	if (isnull)
+		return NULL;
+
+	return statext_histogram_deserialize(DatumGetByteaP(histogram));
+}
+
+/*
+ * Serialize the MV histogram into a bytea value. The basic algorithm is quite
+ * simple, and mostly mimincs the MCV serialization:
+ *
+ * (1) perform deduplication for each attribute (separately)
+ *
+ *   (a) collect all (non-NULL) attribute values from all buckets
+ *   (b) sort the data (using 'lt' from VacAttrStats)
+ *   (c) remove duplicate values from the array
+ *
+ * (2) serialize the arrays into a bytea value
+ *
+ * (3) process all buckets
+ *
+ *   (a) replace min/max values with indexes into the arrays
+ *
+ * Each attribute has to be processed separately, as we're mixing different
+ * datatypes, and we we need to use the right operators to compare/sort them.
+ * We're also mixing pass-by-value and pass-by-ref types, and so on.
+ *
+ * TODO Consider packing boolean flags (NULL) for each item into 'char' or
+ * a longer type (instead of using an array of bool items).
+ */
+static MVHistogram *
+serialize_histogram(MVHistogramBuild * histogram, VacAttrStats **stats)
+{
+	int			dim,
+				i;
+	Size		total_length = 0;
+
+	bytea	   *output = NULL;	/* serialized histogram as bytea */
+	char	   *data = NULL;
+
+	DimensionInfo *info;
+	SortSupport ssup;
+
+	int			nbuckets = histogram->nbuckets;
+	int			ndims = histogram->ndimensions;
+
+	/* allocated for serialized bucket data */
+	int			bucketsize = BUCKET_SIZE(ndims);
+	char	   *bucket = palloc0(bucketsize);
+
+	/* values per dimension (and number of non-NULL values) */
+	Datum	  **values = (Datum **) palloc0(sizeof(Datum *) * ndims);
+	int		   *counts = (int *) palloc0(sizeof(int) * ndims);
+
+	/* info about dimensions (for deserialize) */
+	info = (DimensionInfo *) palloc0(sizeof(DimensionInfo) * ndims);
+
+	/* sort support data */
+	ssup = (SortSupport) palloc0(sizeof(SortSupportData) * ndims);
+
+	/* collect and deduplicate values for each dimension separately */
+	for (dim = 0; dim < ndims; dim++)
+	{
+		int			b;
+		int			count;
+		TypeCacheEntry *type;
+
+		type = lookup_type_cache(stats[dim]->attrtypid, TYPECACHE_LT_OPR);
+
+		/* OID of the data types */
+		histogram->types[dim] = stats[dim]->attrtypid;
+
+		/* keep important info about the data type */
+		info[dim].typlen = stats[dim]->attrtype->typlen;
+		info[dim].typbyval = stats[dim]->attrtype->typbyval;
+
+		/*
+		 * Allocate space for all min/max values, including NULLs (we won't
+		 * use them, but we don't know how many are there), and then collect
+		 * all non-NULL values.
+		 */
+		values[dim] = (Datum *) palloc0(sizeof(Datum) * nbuckets * 2);
+
+		for (b = 0; b < histogram->nbuckets; b++)
+		{
+			/* skip buckets where this dimension is NULL-only */
+			if (!histogram->buckets[b]->nullsonly[dim])
+			{
+				values[dim][counts[dim]] = histogram->buckets[b]->min[dim];
+				counts[dim] += 1;
+
+				values[dim][counts[dim]] = histogram->buckets[b]->max[dim];
+				counts[dim] += 1;
+			}
+		}
+
+		/* there are just NULL values in this dimension */
+		if (counts[dim] == 0)
+			continue;
+
+		/* sort and deduplicate */
+		ssup[dim].ssup_cxt = CurrentMemoryContext;
+		ssup[dim].ssup_collation = DEFAULT_COLLATION_OID;
+		ssup[dim].ssup_nulls_first = false;
+
+		PrepareSortSupportFromOrderingOp(type->lt_opr, &ssup[dim]);
+
+		qsort_arg(values[dim], counts[dim], sizeof(Datum),
+				  compare_scalars_simple, &ssup[dim]);
+
+		/*
+		 * Walk through the array and eliminate duplicitate values, but keep
+		 * the ordering (so that we can do bsearch later). We know there's at
+		 * least 1 item, so we can skip the first element.
+		 */
+		count = 1;				/* number of deduplicated items */
+		for (i = 1; i < counts[dim]; i++)
+		{
+			/* if it's different from the previous value, we need to keep it */
+			if (compare_datums_simple(values[dim][i - 1], values[dim][i], &ssup[dim]) != 0)
+			{
+				/* XXX: not needed if (count == j) */
+				values[dim][count] = values[dim][i];
+				count += 1;
+			}
+		}
+
+		/* make sure we fit into uint16 */
+		Assert(count <= UINT16_MAX);
+
+		/* keep info about the deduplicated count */
+		info[dim].nvalues = count;
+
+		/* compute size of the serialized data */
+		if (info[dim].typlen > 0)
+			/* byval or byref, but with fixed length (name, tid, ...) */
+			info[dim].nbytes = info[dim].nvalues * info[dim].typlen;
+		else if (info[dim].typlen == -1)
+			/* varlena, so just use VARSIZE_ANY */
+			for (i = 0; i < info[dim].nvalues; i++)
+				info[dim].nbytes += VARSIZE_ANY(values[dim][i]);
+		else if (info[dim].typlen == -2)
+			/* cstring, so simply strlen */
+			for (i = 0; i < info[dim].nvalues; i++)
+				info[dim].nbytes += strlen(DatumGetPointer(values[dim][i]));
+		else
+			elog(ERROR, "unknown data type typbyval=%d typlen=%d",
+				 info[dim].typbyval, info[dim].typlen);
+	}
+
+	/*
+	 * Now we finally know how much space we'll need for the serialized
+	 * histogram, as it contains these fields:
+	 *
+	 * - length (4B) for varlena
+	 * - magic (4B)
+	 * - type (4B)
+	 * - ndimensions (4B)
+	 * - nbuckets (4B)
+	 * - info (ndim * sizeof(DimensionInfo)
+	 * - arrays of values for each dimension
+	 * - serialized buckets (nbuckets * bucketsize)
+	 *
+	 * So the 'header' size is 20B + ndim * sizeof(DimensionInfo) and then
+	 * we'll place the data (and buckets).
+	 */
+	total_length = (offsetof(MVHistogram, buckets)
+					+ ndims * sizeof(DimensionInfo)
+					+ nbuckets * bucketsize);
+
+	/* account for the deduplicated data */
+	for (dim = 0; dim < ndims; dim++)
+		total_length += info[dim].nbytes;
+
+	/*
+	 * Enforce arbitrary limit of 1MB on the size of the serialized MCV list.
+	 * This is meant as a protection against someone building MCV list on long
+	 * values (e.g. text documents).
+	 *
+	 * XXX Should we enforce arbitrary limits like this one? Maybe it's not
+	 * even necessary, as long values are usually unique and so won't make it
+	 * into the MCV list in the first place. In the end, we have a 1GB limit
+	 * on bytea values.
+	 */
+	if (total_length > (1024 * 1024))
+		elog(ERROR, "serialized histogram exceeds 1MB (%ld > %d)",
+			 total_length, (1024 * 1024));
+
+	/* allocate space for the serialized histogram list, set header */
+	output = (bytea *) palloc0(total_length);
+
+	/*
+	 * we'll use 'data' to keep track of the place to write data
+	 *
+	 * XXX No VARDATA() here, as MVHistogramBuild includes the length.
+	 */
+	data = (char *) output;
+
+	memcpy(data, histogram, offsetof(MVHistogramBuild, buckets));
+	data += offsetof(MVHistogramBuild, buckets);
+
+	memcpy(data, info, sizeof(DimensionInfo) * ndims);
+	data += sizeof(DimensionInfo) * ndims;
+
+	/* serialize the deduplicated values for all attributes */
+	for (dim = 0; dim < ndims; dim++)
+	{
+#ifdef USE_ASSERT_CHECKING
+		char	   *tmp = data;
+#endif
+		for (i = 0; i < info[dim].nvalues; i++)
+		{
+			Datum		v = values[dim][i];
+
+			if (info[dim].typbyval) /* passed by value */
+			{
+				memcpy(data, &v, info[dim].typlen);
+				data += info[dim].typlen;
+			}
+			else if (info[dim].typlen > 0)	/* pased by reference */
+			{
+				memcpy(data, DatumGetPointer(v), info[dim].typlen);
+				data += info[dim].typlen;
+			}
+			else if (info[dim].typlen == -1)	/* varlena */
+			{
+				memcpy(data, DatumGetPointer(v), VARSIZE_ANY(v));
+				data += VARSIZE_ANY(values[dim][i]);
+			}
+			else if (info[dim].typlen == -2)	/* cstring */
+			{
+				memcpy(data, DatumGetPointer(v), strlen(DatumGetPointer(v)) + 1);
+				data += strlen(DatumGetPointer(v)) + 1;
+			}
+		}
+
+		/* make sure we got exactly the amount of data we expected */
+		Assert((data - tmp) == info[dim].nbytes);
+	}
+
+	/* finally serialize the items, with uint16 indexes instead of the values */
+	for (i = 0; i < nbuckets; i++)
+	{
+		/* don't write beyond the allocated space */
+		Assert(data <= (char *) output + total_length - bucketsize);
+
+		/* reset the values for each item */
+		memset(bucket, 0, bucketsize);
+
+		BUCKET_FREQUENCY(bucket) = histogram->buckets[i]->frequency;
+
+		for (dim = 0; dim < ndims; dim++)
+		{
+			/* do the lookup only for non-NULL values */
+			if (!histogram->buckets[i]->nullsonly[dim])
+			{
+				uint16		idx;
+				Datum	   *v = NULL;
+
+				/* min boundary */
+				v = (Datum *) bsearch_arg(&histogram->buckets[i]->min[dim],
+										  values[dim], info[dim].nvalues, sizeof(Datum),
+										  compare_scalars_simple, &ssup[dim]);
+
+				Assert(v != NULL);	/* serialization or deduplication error */
+
+				/* compute index within the array */
+				idx = (v - values[dim]);
+
+				Assert((idx >= 0) && (idx < info[dim].nvalues));
+
+				BUCKET_MIN_INDEXES(bucket, ndims)[dim] = idx;
+
+				/* max boundary */
+				v = (Datum *) bsearch_arg(&histogram->buckets[i]->max[dim],
+										  values[dim], info[dim].nvalues, sizeof(Datum),
+										  compare_scalars_simple, &ssup[dim]);
+
+				Assert(v != NULL);	/* serialization or deduplication error */
+
+				/* compute index within the array */
+				idx = (v - values[dim]);
+
+				Assert((idx >= 0) && (idx < info[dim].nvalues));
+
+				BUCKET_MAX_INDEXES(bucket, ndims)[dim] = idx;
+			}
+		}
+
+		/* copy flags (nulls, min/max inclusive) */
+		memcpy(BUCKET_NULLS_ONLY(bucket, ndims),
+			   histogram->buckets[i]->nullsonly, sizeof(bool) * ndims);
+
+		memcpy(BUCKET_MIN_INCL(bucket, ndims),
+			   histogram->buckets[i]->min_inclusive, sizeof(bool) * ndims);
+
+		memcpy(BUCKET_MAX_INCL(bucket, ndims),
+			   histogram->buckets[i]->max_inclusive, sizeof(bool) * ndims);
+
+		/* copy the item into the array */
+		memcpy(data, bucket, bucketsize);
+
+		data += bucketsize;
+	}
+
+	/* at this point we expect to match the total_length exactly */
+	Assert((data - (char *) output) == total_length);
+
+	/* free the values/counts arrays here */
+	pfree(counts);
+	pfree(info);
+	pfree(ssup);
+
+	for (dim = 0; dim < ndims; dim++)
+		pfree(values[dim]);
+
+	pfree(values);
+
+	/* make sure the length is correct */
+	SET_VARSIZE(output, total_length);
+
+	return (MVHistogram *)output;
+}
+
+/*
+* Reads serialized histogram into MVHistogram structure.
+ 
+ * Returns histogram in a partially-serialized form (keeps the boundary values
+ * deduplicated, so that it's possible to optimize the estimation part by
+ * caching function call results across buckets etc.).
+ */
+MVHistogram *
+statext_histogram_deserialize(bytea *data)
+{
+	int			dim,
+				i;
+
+	Size		expected_size;
+	char	   *tmp = NULL;
+
+	MVHistogram *histogram;
+	DimensionInfo *info;
+
+	int			nbuckets;
+	int			ndims;
+	int			bucketsize;
+
+	/* temporary deserialization buffer */
+	int			bufflen;
+	char	   *buff;
+	char	   *ptr;
+
+	if (data == NULL)
+		return NULL;
+
+	/*
+	 * We can't possibly deserialize a histogram if there's not even a
+	 * complete header.
+	 */
+	if (VARSIZE_ANY_EXHDR(data) < offsetof(MVHistogram, buckets))
+		elog(ERROR, "invalid histogram size %ld (expected at least %ld)",
+			 VARSIZE_ANY_EXHDR(data), offsetof(MVHistogram, buckets));
+
+	/* read the histogram header */
+	histogram
+		= (MVHistogram *) palloc(sizeof(MVHistogram));
+
+	/* initialize pointer to data (varlena header is included) */
+	tmp = (char *) data;
+
+	/* get the header and perform basic sanity checks */
+	memcpy(histogram, tmp, offsetof(MVHistogram, buckets));
+	tmp += offsetof(MVHistogram, buckets);
+
+	if (histogram->magic != STATS_HIST_MAGIC)
+		elog(ERROR, "invalid histogram magic %d (expected %dd)",
+			 histogram->magic, STATS_HIST_MAGIC);
+
+	if (histogram->type != STATS_HIST_TYPE_BASIC)
+		elog(ERROR, "invalid histogram type %d (expected %dd)",
+			 histogram->type, STATS_HIST_TYPE_BASIC);
+
+	if (histogram->ndimensions == 0)
+		ereport(ERROR,
+				(errcode(ERRCODE_DATA_CORRUPTED),
+				 errmsg("invalid zero-length dimension array in histogram")));
+	else if (histogram->ndimensions > STATS_MAX_DIMENSIONS)
+		ereport(ERROR,
+				(errcode(ERRCODE_DATA_CORRUPTED),
+				 errmsg("invalid length (%d) dimension array in histogram",
+						histogram->ndimensions)));
+
+	if (histogram->nbuckets == 0)
+		ereport(ERROR,
+				(errcode(ERRCODE_DATA_CORRUPTED),
+				 errmsg("invalid zero-length bucket array in histogram")));
+	else if (histogram->nbuckets > STATS_HIST_MAX_BUCKETS)
+		ereport(ERROR,
+				(errcode(ERRCODE_DATA_CORRUPTED),
+				 errmsg("invalid length (%d) bucket array in histogram",
+						histogram->nbuckets)));
+
+	nbuckets = histogram->nbuckets;
+	ndims = histogram->ndimensions;
+	bucketsize = BUCKET_SIZE(ndims);
+
+	/*
+	 * What size do we expect with those parameters (it's incomplete, as we
+	 * yet have to count the array sizes (from DimensionInfo records).
+	 */
+	expected_size = offsetof(MVHistogram, buckets) +
+		ndims * sizeof(DimensionInfo) +
+		(nbuckets * bucketsize);
+
+	/* check that we have at least the DimensionInfo records */
+	if (VARSIZE_ANY(data) < expected_size)
+		elog(ERROR, "invalid histogram size %ld (expected %ld)",
+			 VARSIZE_ANY_EXHDR(data), expected_size);
+
+	/* Now it's safe to access the dimention info. */
+	info = (DimensionInfo *) (tmp);
+	tmp += ndims * sizeof(DimensionInfo);
+
+	/* account for the value arrays */
+	for (dim = 0; dim < ndims; dim++)
+		expected_size += info[dim].nbytes;
+
+	if (VARSIZE_ANY(data) != expected_size)
+		elog(ERROR, "invalid histogram size %ld (expected %ld)",
+			 VARSIZE_ANY_EXHDR(data), expected_size);
+
+	/* looks OK - not corrupted or something */
+
+	/* a single buffer for all the values and counts */
+	bufflen = (sizeof(int) + sizeof(Datum *)) * ndims;
+
+	for (dim = 0; dim < ndims; dim++)
+		/* don't allocate space for byval types, matching Datum */
+		if (!(info[dim].typbyval && (info[dim].typlen == sizeof(Datum))))
+			bufflen += (sizeof(Datum) * info[dim].nvalues);
+
+	/* also, include space for the result, tracking the buckets */
+	bufflen += nbuckets * (sizeof(MVBucket *) + /* bucket pointer */
+						   sizeof(MVBucket));	/* bucket data */
+
+	buff = palloc0(bufflen);
+	ptr = buff;
+
+	histogram->nvalues = (int *) ptr;
+	ptr += (sizeof(int) * ndims);
+
+	histogram->values = (Datum **) ptr;
+	ptr += (sizeof(Datum *) * ndims);
+
+	/*
+	 * XXX This uses pointers to the original data array (the types not passed
+	 * by value), so when someone frees the memory, e.g. by doing something
+	 * like this:
+	 *
+	 *	bytea * data = ... fetch the data from catalog ...
+	 *	MVHistogramBuild histogram = deserialize_histogram(data);
+	 *	pfree(data);
+	 *
+	 * then 'histogram' references the freed memory. Should copy the pieces.
+	 */
+	for (dim = 0; dim < ndims; dim++)
+	{
+#ifdef USE_ASSERT_CHECKING
+		/* remember where data for this dimension starts */
+		char	   *start = tmp;
+#endif
+
+		histogram->nvalues[dim] = info[dim].nvalues;
+
+		if (info[dim].typbyval)
+		{
+			/* passed by value / Datum - simply reuse the array */
+			if (info[dim].typlen == sizeof(Datum))
+			{
+				histogram->values[dim] = (Datum *) tmp;
+				tmp += info[dim].nbytes;
+
+				/* no overflow of input array */
+				Assert(tmp <= start + info[dim].nbytes);
+			}
+			else
+			{
+				histogram->values[dim] = (Datum *) ptr;
+				ptr += (sizeof(Datum) * info[dim].nvalues);
+
+				for (i = 0; i < info[dim].nvalues; i++)
+				{
+					/* just point into the array */
+					memcpy(&histogram->values[dim][i], tmp, info[dim].typlen);
+					tmp += info[dim].typlen;
+
+					/* no overflow of input array */
+					Assert(tmp <= start + info[dim].nbytes);
+				}
+			}
+		}
+		else
+		{
+			/* all the other types need a chunk of the buffer */
+			histogram->values[dim] = (Datum *) ptr;
+			ptr += (sizeof(Datum) * info[dim].nvalues);
+
+			if (info[dim].typlen > 0)
+			{
+				/* pased by reference, but fixed length (name, tid, ...) */
+				for (i = 0; i < info[dim].nvalues; i++)
+				{
+					/* just point into the array */
+					histogram->values[dim][i] = PointerGetDatum(tmp);
+					tmp += info[dim].typlen;
+
+					/* no overflow of input array */
+					Assert(tmp <= start + info[dim].nbytes);
+				}
+			}
+			else if (info[dim].typlen == -1)
+			{
+				/* varlena */
+				for (i = 0; i < info[dim].nvalues; i++)
+				{
+					/* just point into the array */
+					histogram->values[dim][i] = PointerGetDatum(tmp);
+					tmp += VARSIZE_ANY(tmp);
+
+					/* no overflow of input array */
+					Assert(tmp <= start + info[dim].nbytes);
+				}
+			}
+			else if (info[dim].typlen == -2)
+			{
+				/* cstring */
+				for (i = 0; i < info[dim].nvalues; i++)
+				{
+					/* just point into the array */
+					histogram->values[dim][i] = PointerGetDatum(tmp);
+					tmp += (strlen(tmp) + 1);	/* don't forget the \0 */
+
+					/* no overflow of input array */
+					Assert(tmp <= start + info[dim].nbytes);
+				}
+			}
+		}
+
+		/* check we consumed the serialized data for this dimension exactly */
+		Assert((tmp - start) == info[dim].nbytes);
+	}
+
+	/* now deserialize the buckets and point them into the varlena values */
+	histogram->buckets = (MVBucket * *) ptr;
+	ptr += (sizeof(MVBucket *) * nbuckets);
+
+	for (i = 0; i < nbuckets; i++)
+	{
+		MVBucket   *bucket = (MVBucket *) ptr;
+
+		ptr += sizeof(MVBucket);
+
+		bucket->frequency = BUCKET_FREQUENCY(tmp);
+		bucket->nullsonly = BUCKET_NULLS_ONLY(tmp, ndims);
+		bucket->min_inclusive = BUCKET_MIN_INCL(tmp, ndims);
+		bucket->max_inclusive = BUCKET_MAX_INCL(tmp, ndims);
+
+		bucket->min = BUCKET_MIN_INDEXES(tmp, ndims);
+		bucket->max = BUCKET_MAX_INDEXES(tmp, ndims);
+
+		histogram->buckets[i] = bucket;
+
+		Assert(tmp <= (char *) data + VARSIZE_ANY(data));
+
+		tmp += bucketsize;
+	}
+
+	/* at this point we expect to match the total_length exactly */
+	Assert((tmp - (char *) data) == expected_size);
+
+	/* we should exhaust the output buffer exactly */
+	Assert((ptr - buff) == bufflen);
+
+	return histogram;
+}
+
+/*
+ * create_initial_ext_bucket
+ *		Create an initial bucket, covering all the sampled rows.
+ */
+static MVBucketBuild *
+create_initial_ext_bucket(int numrows, HeapTuple *rows, Bitmapset *attrs,
+						  VacAttrStats **stats)
+{
+	int			i;
+	int			numattrs = bms_num_members(attrs);
+
+	/* TODO allocate bucket as a single piece, including all the fields. */
+	MVBucketBuild *bucket = (MVBucketBuild *) palloc0(sizeof(MVBucketBuild));
+
+	Assert(numrows > 0);
+	Assert(rows != NULL);
+	Assert((numattrs >= 2) && (numattrs <= STATS_MAX_DIMENSIONS));
+
+	/* allocate the per-dimension arrays */
+
+	/* flags for null-only dimensions */
+	bucket->nullsonly = (bool *) palloc0(numattrs * sizeof(bool));
+
+	/* inclusiveness boundaries - lower/upper bounds */
+	bucket->min_inclusive = (bool *) palloc0(numattrs * sizeof(bool));
+	bucket->max_inclusive = (bool *) palloc0(numattrs * sizeof(bool));
+
+	/* lower/upper boundaries */
+	bucket->min = (Datum *) palloc0(numattrs * sizeof(Datum));
+	bucket->max = (Datum *) palloc0(numattrs * sizeof(Datum));
+
+	/* number of distinct values (per dimension) */
+	bucket->ndistincts = (uint32 *) palloc0(numattrs * sizeof(uint32));
+
+	/* all the sample rows fall into the initial bucket */
+	bucket->numrows = numrows;
+	bucket->rows = rows;
+
+	/*
+	 * Update the number of ndistinct combinations in the bucket (which we use
+	 * when selecting bucket to partition), and then number of distinct values
+	 * for each partition (which we use when choosing which dimension to
+	 * split).
+	 */
+	update_bucket_ndistinct(bucket, attrs, stats);
+
+	/* Update ndistinct (and also set min/max) for all dimensions. */
+	for (i = 0; i < numattrs; i++)
+		update_dimension_ndistinct(bucket, i, attrs, stats, true);
+
+	return bucket;
+}
+
+/*
+ * Choose the bucket to partition next.
+ *
+ * The current criteria is rather simple, chosen so that the algorithm produces
+ * buckets with about equal frequency and regular size. We select the bucket
+ * with the highest number of distinct values, and then split it by the longest
+ * dimension.
+ *
+ * The distinct values are uniformly mapped to [0,1] interval, and this is used
+ * to compute length of the value range.
+ *
+ * NOTE: This is not the same array used for deduplication, as this contains
+ *		 values for all the tuples from the sample, not just the boundary values.
+ *
+ * Returns either pointer to the bucket selected to be partitioned, or NULL if
+ * there are no buckets that may be split (e.g. if all buckets are too small
+ * or contain too few distinct values).
+ *
+ *
+ * Tricky example
+ * --------------
+ *
+ * Consider this table:
+ *
+ *	   CREATE TABLE t AS SELECT i AS a, i AS b
+ *						   FROM generate_series(1,1000000) s(i);
+ *
+ *	   CREATE STATISTICS s1 ON t (a,b) WITH (histogram);
+ *
+ *	   ANALYZE t;
+ *
+ * It's a very specific (and perhaps artificial) example, because every bucket
+ * always has exactly the same number of distinct values in all dimensions,
+ * which makes the partitioning tricky.
+ *
+ * Then:
+ *
+ *	   SELECT * FROM t WHERE (a < 100) AND (b < 100);
+ *
+ * is estimated to return ~120 rows, while in reality it returns only 99.
+ *
+ *							 QUERY PLAN
+ *	   -------------------------------------------------------------
+ *		Seq Scan on t  (cost=0.00..19425.00 rows=117 width=8)
+ *					   (actual time=0.129..82.776 rows=99 loops=1)
+ *		  Filter: ((a < 100) AND (b < 100))
+ *		  Rows Removed by Filter: 999901
+ *		Planning time: 1.286 ms
+ *		Execution time: 82.984 ms
+ *	   (5 rows)
+ *
+ * So this estimate is reasonably close. Let's change the query to OR clause:
+ *
+ *	   SELECT * FROM t WHERE (a < 100) OR (b < 100);
+ *
+ *							 QUERY PLAN
+ *	   -------------------------------------------------------------
+ *		Seq Scan on t  (cost=0.00..19425.00 rows=8100 width=8)
+ *					   (actual time=0.145..99.910 rows=99 loops=1)
+ *		  Filter: ((a < 100) OR (b < 100))
+ *		  Rows Removed by Filter: 999901
+ *		Planning time: 1.578 ms
+ *		Execution time: 100.132 ms
+ *	   (5 rows)
+ *
+ * That's clearly a much worse estimate. This happens because the histogram
+ * contains buckets like this:
+ *
+ *	   bucket 592  [3 30310] [30134 30593] => [0.000233]
+ *
+ * i.e. the length of "a" dimension is (30310-3)=30307, while the length of "b"
+ * is (30593-30134)=459. So the "b" dimension is much narrower than "a".
+ * Of course, there are also buckets where "b" is the wider dimension.
+ *
+ * This is partially mitigated by selecting the "longest" dimension but that
+ * only happens after we already selected the bucket. So if we never select the
+ * bucket, this optimization does not apply.
+ *
+ * The other reason why this particular example behaves so poorly is due to the
+ * way we actually split the selected bucket. We do attempt to divide the bucket
+ * into two parts containing about the same number of tuples, but that does not
+ * too well when most of the tuples is squashed on one side of the bucket.
+ *
+ * For example for columns with data on the diagonal (i.e. when a=b), we end up
+ * with a narrow bucket on the diagonal and a huge bucket overing the remaining
+ * part (with much lower density).
+ *
+ * So perhaps we need two partitioning strategies - one aiming to split buckets
+ * with high frequency (number of sampled rows), the other aiming to split
+ * "large" buckets. And alternating between them, somehow.
+ *
+ * TODO Consider using similar lower boundary for row count as for simple
+ * histograms, i.e. 300 tuples per bucket.
+ */
+static MVBucketBuild *
+select_bucket_to_partition(int nbuckets, MVBucketBuild * *buckets)
+{
+	int			i;
+	int			numrows = 0;
+	MVBucketBuild *bucket = NULL;
+
+	for (i = 0; i < nbuckets; i++)
+	{
+		/* if the number of rows is higher, use this bucket */
+		if ((buckets[i]->ndistinct > 2) &&
+			(buckets[i]->numrows > numrows) &&
+			(buckets[i]->numrows >= MIN_BUCKET_ROWS))
+		{
+			bucket = buckets[i];
+			numrows = buckets[i]->numrows;
+		}
+	}
+
+	/* may be NULL if there are not buckets with (ndistinct>1) */
+	return bucket;
+}
+
+/*
+ * A simple bucket partitioning implementation - we choose the longest bucket
+ * dimension, measured using the array of distinct values built at the very
+ * beginning of the build.
+ *
+ * We map all the distinct values to a [0,1] interval, uniformly distributed,
+ * and then use this to measure length. It's essentially a number of distinct
+ * values within the range, normalized to [0,1].
+ *
+ * Then we choose a 'middle' value splitting the bucket into two parts with
+ * roughly the same frequency.
+ *
+ * This splits the bucket by tweaking the existing one, and returning the new
+ * bucket (essentially shrinking the existing one in-place and returning the
+ * other "half" as a new bucket). The caller is responsible for adding the new
+ * bucket into the list of buckets.
+ *
+ * There are multiple histogram options, centered around the partitioning
+ * criteria, specifying both how to choose a bucket and the dimension most in
+ * need of a split. For a nice summary and general overview, see "rK-Hist : an
+ * R-Tree based histogram for multi-dimensional selectivity estimation" thesis
+ * by J. A. Lopez, Concordia University, p.34-37 (and possibly p. 32-34 for
+ * explanation of the terms).
+ *
+ * It requires care to prevent splitting only one dimension and not splitting
+ * another one at all (which might happen easily in case of strongly dependent
+ * columns - e.g. y=x). The current algorithm minimizes this, but may still
+ * happen for perfectly dependent examples (when all the dimensions have equal
+ * length, the first one will be selected).
+ *
+ * TODO Should probably consider statistics target for the columns (e.g.
+ * to split dimensions with higher statistics target more frequently).
+ */
+static MVBucketBuild *
+partition_bucket(MVBucketBuild * bucket, Bitmapset *attrs,
+				 VacAttrStats **stats,
+				 int *ndistvalues, Datum **distvalues)
+{
+	int			i;
+	int			dimension;
+	int			numattrs = bms_num_members(attrs);
+
+	Datum		split_value;
+	MVBucketBuild *new_bucket;
+
+	/* needed for sort, when looking for the split value */
+	bool		isNull;
+	int			nvalues = 0;
+	TypeCacheEntry *type;
+	ScalarItem *values;
+	SortSupportData ssup;
+	int		   *attnums;
+
+	int			nrows = 1;		/* number of rows below current value */
+	double		delta;
+
+	/* needed when splitting the values */
+	HeapTuple  *oldrows = bucket->rows;
+	int			oldnrows = bucket->numrows;
+
+	values = (ScalarItem *) palloc0(bucket->numrows * sizeof(ScalarItem));
+
+	/*
+	 * We can't split buckets with a single distinct value (this also
+	 * disqualifies NULL-only dimensions). Also, there has to be multiple
+	 * sample rows (otherwise, how could there be more distinct values).
+	 */
+	Assert(bucket->ndistinct > 1);
+	Assert(bucket->numrows > 1);
+	Assert((numattrs >= 2) && (numattrs <= STATS_MAX_DIMENSIONS));
+
+	/* Look for the next dimension to split. */
+	delta = 0.0;
+	dimension = -1;
+
+	for (i = 0; i < numattrs; i++)
+	{
+		Datum	   *a,
+				   *b;
+
+		type = lookup_type_cache(stats[i]->attrtypid, TYPECACHE_LT_OPR);
+
+		/* initialize sort support, etc. */
+		memset(&ssup, 0, sizeof(ssup));
+		ssup.ssup_cxt = CurrentMemoryContext;
+
+		/* We always use the default collation for statistics */
+		ssup.ssup_collation = DEFAULT_COLLATION_OID;
+		ssup.ssup_nulls_first = false;
+
+		PrepareSortSupportFromOrderingOp(type->lt_opr, &ssup);
+
+		/* can't split NULL-only dimension */
+		if (bucket->nullsonly[i])
+			continue;
+
+		/* can't split dimension with a single ndistinct value */
+		if (bucket->ndistincts[i] <= 1)
+			continue;
+
+		/* search for min boundary in the distinct list */
+		a = (Datum *) bsearch_arg(&bucket->min[i],
+								  distvalues[i], ndistvalues[i],
+								  sizeof(Datum), compare_scalars_simple, &ssup);
+
+		b = (Datum *) bsearch_arg(&bucket->max[i],
+								  distvalues[i], ndistvalues[i],
+								  sizeof(Datum), compare_scalars_simple, &ssup);
+
+		/* if this dimension is 'larger' then partition by it */
+		if (((b - a) * 1.0 / ndistvalues[i]) > delta)
+		{
+			delta = ((b - a) * 1.0 / ndistvalues[i]);
+			dimension = i;
+		}
+	}
+
+	/*
+	 * If we haven't found a dimension here, we've done something wrong in
+	 * select_bucket_to_partition.
+	 */
+	Assert(dimension != -1);
+
+	/*
+	 * Walk through the selected dimension, collect and sort the values and
+	 * then choose the value to use as the new boundary.
+	 */
+	type = lookup_type_cache(stats[dimension]->attrtypid, TYPECACHE_LT_OPR);
+
+	/* initialize sort support, etc. */
+	memset(&ssup, 0, sizeof(ssup));
+	ssup.ssup_cxt = CurrentMemoryContext;
+
+	/* We always use the default collation for statistics */
+	ssup.ssup_collation = DEFAULT_COLLATION_OID;
+	ssup.ssup_nulls_first = false;
+
+	PrepareSortSupportFromOrderingOp(type->lt_opr, &ssup);
+
+	attnums = build_attnums(attrs);
+
+	for (i = 0; i < bucket->numrows; i++)
+	{
+		/*
+		 * remember the index of the sample row, to make the partitioning
+		 * simpler
+		 */
+		values[nvalues].value = heap_getattr(bucket->rows[i], attnums[dimension],
+											 stats[dimension]->tupDesc, &isNull);
+		values[nvalues].tupno = i;
+
+		/* no NULL values allowed here (we never split null-only dimension) */
+		Assert(!isNull);
+
+		nvalues++;
+	}
+
+	/* sort the array of values */
+	qsort_arg((void *) values, nvalues, sizeof(ScalarItem),
+			  compare_scalars_partition, (void *) &ssup);
+
+	/*
+	 * We know there are bucket->ndistincts[dimension] distinct values in this
+	 * dimension, and we want to split this into half, so walk through the
+	 * array and stop once we see (ndistinct/2) values.
+	 *
+	 * We always choose the "next" value, i.e. (n/2+1)-th distinct value, and
+	 * use it as an exclusive upper boundary (and inclusive lower boundary).
+	 *
+	 * TODO Maybe we should use "average" of the two middle distinct values
+	 * (at least for even distinct counts), but that would require being able
+	 * to do an average (which does not work for non-numeric types).
+	 *
+	 * TODO Another option is to look for a split that'd give about 50% tuples
+	 * (not distinct values) in each partition. That might work better when
+	 * there are a few very frequent values, and many rare ones.
+	 */
+	delta = bucket->numrows;
+	split_value = values[0].value;
+
+	for (i = 1; i < bucket->numrows; i++)
+	{
+		if (values[i].value != values[i - 1].value)
+		{
+			/* are we closer to splitting the bucket in half? */
+			if (fabs(i - bucket->numrows / 2.0) < delta)
+			{
+				/* let's assume we'll use this value for the split */
+				split_value = values[i].value;
+				delta = fabs(i - bucket->numrows / 2.0);
+				nrows = i;
+			}
+		}
+	}
+
+	Assert(nrows > 0);
+	Assert(nrows < bucket->numrows);
+
+	/*
+	 * create the new bucket as a (incomplete) copy of the one being
+	 * partitioned.
+	 */
+	new_bucket = copy_ext_bucket(bucket, numattrs);
+
+	/*
+	 * Do the actual split of the chosen dimension, using the split value as
+	 * the upper bound for the existing bucket, and lower bound for the new
+	 * one.
+	 */
+	bucket->max[dimension] = split_value;
+	new_bucket->min[dimension] = split_value;
+
+	/*
+	 * We also treat only one side of the new boundary as inclusive, in the
+	 * bucket where it happens to be the upper boundary. We never set the
+	 * min_inclusive[] to false anywhere, but we set it to true anyway.
+	 */
+	bucket->max_inclusive[dimension] = false;
+	new_bucket->min_inclusive[dimension] = true;
+
+	/*
+	 * Redistribute the sample tuples using the 'ScalarItem->tupno' index. We
+	 * know 'nrows' rows should remain in the original bucket and the rest
+	 * goes to the new one.
+	 */
+	bucket->numrows = nrows;
+	new_bucket->numrows = (oldnrows - nrows);
+
+	bucket->rows = (HeapTuple *) palloc0(bucket->numrows * sizeof(HeapTuple));
+	new_bucket->rows = (HeapTuple *) palloc0(new_bucket->numrows * sizeof(HeapTuple));
+
+	/*
+	 * The first nrows should go to the first bucket, the rest should go to
+	 * the new one. Use the tupno field to get the actual HeapTuple row from
+	 * the original array of sample rows.
+	 */
+	for (i = 0; i < nrows; i++)
+		memcpy(&bucket->rows[i], &oldrows[values[i].tupno], sizeof(HeapTuple));
+
+	for (i = nrows; i < oldnrows; i++)
+		memcpy(&new_bucket->rows[i - nrows], &oldrows[values[i].tupno], sizeof(HeapTuple));
+
+	/* update ndistinct values for the buckets (total and per dimension) */
+	update_bucket_ndistinct(bucket, attrs, stats);
+	update_bucket_ndistinct(new_bucket, attrs, stats);
+
+	/*
+	 * TODO We don't need to do this for the dimension we used for split,
+	 * because we know how many distinct values went to each partition.
+	 */
+	for (i = 0; i < numattrs; i++)
+	{
+		update_dimension_ndistinct(bucket, i, attrs, stats, false);
+		update_dimension_ndistinct(new_bucket, i, attrs, stats, false);
+	}
+
+	pfree(oldrows);
+	pfree(values);
+
+	return new_bucket;
+}
+
+/*
+ * Copy a histogram bucket. The copy does not include the build-time data, i.e.
+ * sampled rows etc.
+ */
+static MVBucketBuild *
+copy_ext_bucket(MVBucketBuild * bucket, uint32 ndimensions)
+{
+	/* TODO allocate as a single piece (including all the fields) */
+	MVBucketBuild *new_bucket = (MVBucketBuild *) palloc0(sizeof(MVBucketBuild));
+
+	/*
+	 * Copy only the attributes that will stay the same after the split, and
+	 * we'll recompute the rest after the split.
+	 */
+
+	/* allocate the per-dimension arrays */
+	new_bucket->nullsonly = (bool *) palloc0(ndimensions * sizeof(bool));
+
+	/* inclusiveness boundaries - lower/upper bounds */
+	new_bucket->min_inclusive = (bool *) palloc0(ndimensions * sizeof(bool));
+	new_bucket->max_inclusive = (bool *) palloc0(ndimensions * sizeof(bool));
+
+	/* lower/upper boundaries */
+	new_bucket->min = (Datum *) palloc0(ndimensions * sizeof(Datum));
+	new_bucket->max = (Datum *) palloc0(ndimensions * sizeof(Datum));
+
+	/* copy data */
+	memcpy(new_bucket->nullsonly, bucket->nullsonly, ndimensions * sizeof(bool));
+
+	memcpy(new_bucket->min_inclusive, bucket->min_inclusive, ndimensions * sizeof(bool));
+	memcpy(new_bucket->min, bucket->min, ndimensions * sizeof(Datum));
+
+	memcpy(new_bucket->max_inclusive, bucket->max_inclusive, ndimensions * sizeof(bool));
+	memcpy(new_bucket->max, bucket->max, ndimensions * sizeof(Datum));
+
+	/* allocate and copy the interesting part of the build data */
+	new_bucket->ndistincts = (uint32 *) palloc0(ndimensions * sizeof(uint32));
+
+	return new_bucket;
+}
+
+/*
+ * Counts the number of distinct values in the bucket. This just copies the
+ * Datum values into a simple array, and sorts them using memcmp-based
+ * comparator. That means it only works for pass-by-value data types (assuming
+ * they don't use collations etc.)
+ */
+static void
+update_bucket_ndistinct(MVBucketBuild * bucket, Bitmapset *attrs, VacAttrStats **stats)
+{
+	int			i;
+	int			numattrs = bms_num_members(attrs);
+	int			numrows = bucket->numrows;
+
+	MultiSortSupport mss = multi_sort_init(numattrs);
+	int		   *attnums;
+	SortItem   *items;
+
+	attnums = build_attnums(attrs);
+
+	/* prepare the sort function for the first dimension */
+	for (i = 0; i < numattrs; i++)
+	{
+		VacAttrStats *colstat = stats[i];
+		TypeCacheEntry *type;
+
+		type = lookup_type_cache(colstat->attrtypid, TYPECACHE_LT_OPR);
+		if (type->lt_opr == InvalidOid) /* shouldn't happen */
+			elog(ERROR, "cache lookup failed for ordering operator for type %u",
+				 colstat->attrtypid);
+
+		multi_sort_add_dimension(mss, i, type->lt_opr);
+	}
+
+	/*
+	 * build an array of SortItem(s) sorted using the multi-sort support
+	 *
+	 * XXX This relies on all stats entries pointing to the same tuple
+	 * descriptor. Not sure if that might not be the case.
+	 */
+	items = build_sorted_items(numrows, bucket->rows, stats[0]->tupDesc,
+							   mss, numattrs, attnums);
+
+	bucket->ndistinct = 1;
+
+	for (i = 1; i < numrows; i++)
+		if (multi_sort_compare(&items[i], &items[i - 1], mss) != 0)
+			bucket->ndistinct += 1;
+
+	pfree(items);
+}
+
+/*
+ * Count distinct values per bucket dimension.
+ */
+static void
+update_dimension_ndistinct(MVBucketBuild * bucket, int dimension, Bitmapset *attrs,
+						   VacAttrStats **stats, bool update_boundaries)
+{
+	int			j;
+	int			nvalues = 0;
+	bool		isNull;
+	Datum	   *values;
+	SortSupportData ssup;
+	TypeCacheEntry *type;
+	int		   *attnums;
+
+	values = (Datum *) palloc0(bucket->numrows * sizeof(Datum));
+	type = lookup_type_cache(stats[dimension]->attrtypid, TYPECACHE_LT_OPR);
+
+	/* we may already know this is a NULL-only dimension */
+	if (bucket->nullsonly[dimension])
+		bucket->ndistincts[dimension] = 1;
+
+	memset(&ssup, 0, sizeof(ssup));
+	ssup.ssup_cxt = CurrentMemoryContext;
+
+	/* We always use the default collation for statistics */
+	ssup.ssup_collation = DEFAULT_COLLATION_OID;
+	ssup.ssup_nulls_first = false;
+
+	PrepareSortSupportFromOrderingOp(type->lt_opr, &ssup);
+
+	attnums = build_attnums(attrs);
+
+	for (j = 0; j < bucket->numrows; j++)
+	{
+		values[nvalues] = heap_getattr(bucket->rows[j], attnums[dimension],
+									   stats[dimension]->tupDesc, &isNull);
+
+		/* ignore NULL values */
+		if (!isNull)
+			nvalues++;
+	}
+
+	/* there's always at least 1 distinct value (may be NULL) */
+	bucket->ndistincts[dimension] = 1;
+
+	/*
+	 * if there are only NULL values in the column, mark it so and continue
+	 * with the next one
+	 */
+	if (nvalues == 0)
+	{
+		pfree(values);
+		bucket->nullsonly[dimension] = true;
+		return;
+	}
+
+	/* sort the array (pass-by-value datum */
+	qsort_arg((void *) values, nvalues, sizeof(Datum),
+			  compare_scalars_simple, (void *) &ssup);
+
+	/*
+	 * Update min/max boundaries to the smallest bounding box. Generally, this
+	 * needs to be done only when constructing the initial bucket.
+	 */
+	if (update_boundaries)
+	{
+		/* store the min/max values */
+		bucket->min[dimension] = values[0];
+		bucket->min_inclusive[dimension] = true;
+
+		bucket->max[dimension] = values[nvalues - 1];
+		bucket->max_inclusive[dimension] = true;
+	}
+
+	/*
+	 * Walk through the array and count distinct values by comparing
+	 * succeeding values.
+	 */
+	for (j = 1; j < nvalues; j++)
+	{
+		if (compare_datums_simple(values[j - 1], values[j], &ssup) != 0)
+			bucket->ndistincts[dimension] += 1;
+	}
+
+	pfree(values);
+}
+
+/*
+ * A properly built histogram must not contain buckets mixing NULL and non-NULL
+ * values in a single dimension. Each dimension may either be marked as 'nulls
+ * only', and thus containing only NULL values, or it must not contain any NULL
+ * values.
+ *
+ * Therefore, if the sample contains NULL values in any of the columns, it's
+ * necessary to build those NULL-buckets. This is done in an iterative way
+ * using this algorithm, operating on a single bucket:
+ *
+ *	   (1) Check that all dimensions are well-formed (not mixing NULL and
+ *		   non-NULL values).
+ *
+ *	   (2) If all dimensions are well-formed, terminate.
+ *
+ *	   (3) If the dimension contains only NULL values, but is not marked as
+ *		   NULL-only, mark it as NULL-only and run the algorithm again (on
+ *		   this bucket).
+ *
+ *	   (4) If the dimension mixes NULL and non-NULL values, split the bucket
+ *		   into two parts - one with NULL values, one with non-NULL values
+ *		   (replacing the current one). Then run the algorithm on both buckets.
+ *
+ * This is executed in a recursive manner, but the number of executions should
+ * be quite low - limited by the number of NULL-buckets. Also, in each branch
+ * the number of nested calls is limited by the number of dimensions
+ * (attributes) of the histogram.
+ *
+ * At the end, there should be buckets with no mixed dimensions. The number of
+ * buckets produced by this algorithm is rather limited - with N dimensions,
+ * there may be only 2^N such buckets (each dimension may be either NULL or
+ * non-NULL). So with 8 dimensions (current value of STATS_MAX_DIMENSIONS)
+ * there may be only 256 such buckets.
+ *
+ * After this, a 'regular' bucket-split algorithm shall run, further optimizing
+ * the histogram.
+ */
+static void
+create_null_buckets(MVHistogramBuild * histogram, int bucket_idx,
+					Bitmapset *attrs, VacAttrStats **stats)
+{
+	int			i,
+				j;
+	int			null_dim = -1;
+	int			null_count = 0;
+	bool		null_found = false;
+	MVBucketBuild *bucket,
+			   *null_bucket;
+	int			null_idx,
+				curr_idx;
+	int		   *attnums;
+
+	/* remember original values from the bucket */
+	int			numrows;
+	HeapTuple  *oldrows = NULL;
+
+	Assert(bucket_idx < histogram->nbuckets);
+	Assert(histogram->ndimensions == bms_num_members(attrs));
+
+	bucket = histogram->buckets[bucket_idx];
+
+	numrows = bucket->numrows;
+	oldrows = bucket->rows;
+
+	attnums = build_attnums(attrs);
+
+	/*
+	 * Walk through all rows / dimensions, and stop once we find NULL in a
+	 * dimension not yet marked as NULL-only.
+	 */
+	for (i = 0; i < bucket->numrows; i++)
+	{
+		for (j = 0; j < histogram->ndimensions; j++)
+		{
+			/* Is this a NULL-only dimension? If yes, skip. */
+			if (bucket->nullsonly[j])
+				continue;
+
+			/* found a NULL in that dimension? */
+			if (heap_attisnull(bucket->rows[i], attnums[j],
+							   stats[j]->tupDesc))
+			{
+				null_found = true;
+				null_dim = j;
+				break;
+			}
+		}
+
+		/* terminate if we found attribute with NULL values */
+		if (null_found)
+			break;
+	}
+
+	/* no regular dimension contains NULL values => we're done */
+	if (!null_found)
+		return;
+
+	/* walk through the rows again, count NULL values in 'null_dim' */
+	for (i = 0; i < bucket->numrows; i++)
+	{
+		if (heap_attisnull(bucket->rows[i], attnums[null_dim],
+						   stats[null_dim]->tupDesc))
+			null_count += 1;
+	}
+
+	Assert(null_count <= bucket->numrows);
+
+	/*
+	 * If (null_count == numrows) the dimension already is NULL-only, but is
+	 * not yet marked like that. It's enough to mark it and repeat the process
+	 * recursively (until we run out of dimensions).
+	 */
+	if (null_count == bucket->numrows)
+	{
+		bucket->nullsonly[null_dim] = true;
+		create_null_buckets(histogram, bucket_idx, attrs, stats);
+		return;
+	}
+
+	/*
+	 * We have to split the bucket into two - one with NULL values in the
+	 * dimension, one with non-NULL values. We don't need to sort the data or
+	 * anything, but otherwise it's similar to what partition_bucket() does.
+	 */
+
+	/* create bucket with NULL-only dimension 'dim' */
+	null_bucket = copy_ext_bucket(bucket, histogram->ndimensions);
+
+	/* remember the current array info */
+	oldrows = bucket->rows;
+	numrows = bucket->numrows;
+
+	/* we'll keep non-NULL values in the current bucket */
+	bucket->numrows = (numrows - null_count);
+	bucket->rows
+		= (HeapTuple *) palloc0(bucket->numrows * sizeof(HeapTuple));
+
+	/* and the NULL values will go to the new one */
+	null_bucket->numrows = null_count;
+	null_bucket->rows
+		= (HeapTuple *) palloc0(null_bucket->numrows * sizeof(HeapTuple));
+
+	/* mark the dimension as NULL-only (in the new bucket) */
+	null_bucket->nullsonly[null_dim] = true;
+
+	/* walk through the sample rows and distribute them accordingly */
+	null_idx = 0;
+	curr_idx = 0;
+	for (i = 0; i < numrows; i++)
+	{
+		if (heap_attisnull(oldrows[i], attnums[null_dim],
+						   stats[null_dim]->tupDesc))
+			/* NULL => copy to the new bucket */
+			memcpy(&null_bucket->rows[null_idx++], &oldrows[i],
+				   sizeof(HeapTuple));
+		else
+			memcpy(&bucket->rows[curr_idx++], &oldrows[i],
+				   sizeof(HeapTuple));
+	}
+
+	/* update ndistinct values for the buckets (total and per dimension) */
+	update_bucket_ndistinct(bucket, attrs, stats);
+	update_bucket_ndistinct(null_bucket, attrs, stats);
+
+	/*
+	 * TODO We don't need to do this for the dimension we used for split,
+	 * because we know how many distinct values went to each bucket (NULL is
+	 * not a value, so NULL buckets get 0, and the other bucket got all the
+	 * distinct values).
+	 */
+	for (i = 0; i < histogram->ndimensions; i++)
+	{
+		update_dimension_ndistinct(bucket, i, attrs, stats, false);
+		update_dimension_ndistinct(null_bucket, i, attrs, stats, false);
+	}
+
+	pfree(oldrows);
+
+	/* add the NULL bucket to the histogram */
+	histogram->buckets[histogram->nbuckets++] = null_bucket;
+
+	/*
+	 * And now run the function recursively on both buckets (the new one
+	 * first, because the call may change number of buckets, and it's used as
+	 * an index).
+	 */
+	create_null_buckets(histogram, (histogram->nbuckets - 1), attrs, stats);
+	create_null_buckets(histogram, bucket_idx, attrs, stats);
+}
+
+/*
+ * SRF with details about buckets of a histogram:
+ *
+ * - bucket ID (0...nbuckets)
+ * - min values (string array)
+ * - max values (string array)
+ * - nulls only (boolean array)
+ * - min inclusive flags (boolean array)
+ * - max inclusive flags (boolean array)
+ * - frequency (double precision)
+ *
+ * The input is the OID of the statistics, and there are no rows returned if the
+ * statistics contains no histogram (or if there's no statistics for the OID).
+ *
+ * The second parameter (type) determines what values will be returned
+ * in the (minvals,maxvals). There are three possible values:
+ *
+ * 0 (actual values)
+ * -----------------
+ *	  - prints actual values
+ *	  - using the output function of the data type (as string)
+ *	  - handy for investigating the histogram
+ *
+ * 1 (distinct index)
+ * ------------------
+ *	  - prints index of the distinct value (into the serialized array)
+ *	  - makes it easier to spot neighbor buckets, etc.
+ *	  - handy for plotting the histogram
+ *
+ * 2 (normalized distinct index)
+ * -----------------------------
+ *	  - prints index of the distinct value, but normalized into [0,1]
+ *	  - similar to 1, but shows how 'long' the bucket range is
+ *	  - handy for plotting the histogram
+ *
+ * When plotting the histogram, be careful as the (1) and (2) options skew the
+ * lengths by distributing the distinct values uniformly. For data types
+ * without a clear meaning of 'distance' (e.g. strings) that is not a big deal,
+ * but for numbers it may be confusing.
+ */
+PG_FUNCTION_INFO_V1(pg_histogram_buckets);
+
+#define OUTPUT_FORMAT_RAW		0
+#define OUTPUT_FORMAT_INDEXES	1
+#define OUTPUT_FORMAT_DISTINCT	2
+
+Datum
+pg_histogram_buckets(PG_FUNCTION_ARGS)
+{
+	FuncCallContext *funcctx;
+	int			call_cntr;
+	int			max_calls;
+	TupleDesc	tupdesc;
+	AttInMetadata *attinmeta;
+
+	int			otype = PG_GETARG_INT32(1);
+
+	if ((otype < 0) || (otype > 2))
+		elog(ERROR, "invalid output type specified");
+
+	/* stuff done only on the first call of the function */
+	if (SRF_IS_FIRSTCALL())
+	{
+		MemoryContext oldcontext;
+		MVHistogram *histogram;
+
+		/* create a function context for cross-call persistence */
+		funcctx = SRF_FIRSTCALL_INIT();
+
+		/* switch to memory context appropriate for multiple function calls */
+		oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
+
+		histogram = statext_histogram_deserialize(PG_GETARG_BYTEA_P(0));
+
+		funcctx->user_fctx = histogram;
+
+		/* total number of tuples to be returned */
+		funcctx->max_calls = 0;
+		if (funcctx->user_fctx != NULL)
+			funcctx->max_calls = histogram->nbuckets;
+
+		/* Build a tuple descriptor for our result type */
+		if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
+			ereport(ERROR,
+					(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
+					 errmsg("function returning record called in context "
+							"that cannot accept type record")));
+
+		/*
+		 * generate attribute metadata needed later to produce tuples from raw
+		 * C strings
+		 */
+		attinmeta = TupleDescGetAttInMetadata(tupdesc);
+		funcctx->attinmeta = attinmeta;
+
+		MemoryContextSwitchTo(oldcontext);
+	}
+
+	/* stuff done on every call of the function */
+	funcctx = SRF_PERCALL_SETUP();
+
+	call_cntr = funcctx->call_cntr;
+	max_calls = funcctx->max_calls;
+	attinmeta = funcctx->attinmeta;
+
+	if (call_cntr < max_calls)	/* do when there is more left to send */
+	{
+		char	  **values;
+		HeapTuple	tuple;
+		Datum		result;
+		double		bucket_volume = 1.0;
+		StringInfo	bufs;
+
+		char	   *format;
+		int			i;
+
+		Oid		   *outfuncs;
+		FmgrInfo   *fmgrinfo;
+
+		MVHistogram *histogram;
+		MVBucket   *bucket;
+
+		histogram = (MVHistogram *) funcctx->user_fctx;
+
+		Assert(call_cntr < histogram->nbuckets);
+
+		bucket = histogram->buckets[call_cntr];
+
+		/*
+		 * The scalar values will be formatted directly, using snprintf.
+		 *
+		 * The 'array' values will be formatted through StringInfo.
+		 */
+		values = (char **) palloc0(9 * sizeof(char *));
+		bufs = (StringInfo) palloc0(9 * sizeof(StringInfoData));
+
+		values[0] = (char *) palloc(64 * sizeof(char));
+
+		initStringInfo(&bufs[1]);	/* lower boundaries */
+		initStringInfo(&bufs[2]);	/* upper boundaries */
+		initStringInfo(&bufs[3]);	/* nulls-only */
+		initStringInfo(&bufs[4]);	/* lower inclusive */
+		initStringInfo(&bufs[5]);	/* upper inclusive */
+
+		values[6] = (char *) palloc(64 * sizeof(char));
+		values[7] = (char *) palloc(64 * sizeof(char));
+		values[8] = (char *) palloc(64 * sizeof(char));
+
+		/* we need to do this only when printing the actual values */
+		outfuncs = (Oid *) palloc0(sizeof(Oid) * histogram->ndimensions);
+		fmgrinfo = (FmgrInfo *) palloc0(sizeof(FmgrInfo) * histogram->ndimensions);
+
+		/*
+		 * lookup output functions for all histogram dimensions
+		 *
+		 * XXX This might be one in the first call and stored in user_fctx.
+		 */
+		for (i = 0; i < histogram->ndimensions; i++)
+		{
+			bool		isvarlena;
+
+			getTypeOutputInfo(histogram->types[i], &outfuncs[i], &isvarlena);
+
+			fmgr_info(outfuncs[i], &fmgrinfo[i]);
+		}
+
+		snprintf(values[0], 64, "%d", call_cntr);	/* bucket ID */
+
+		/*
+		 * for the arrays of lower/upper boundaries, formated according to
+		 * otype
+		 */
+		for (i = 0; i < histogram->ndimensions; i++)
+		{
+			Datum	   *vals = histogram->values[i];
+
+			uint16		minidx = bucket->min[i];
+			uint16		maxidx = bucket->max[i];
+
+			int			d = 1;
+
+			/*
+			 * compute bucket volume, using distinct values as a measure
+			 *
+			 * XXX Not really sure what to do for NULL dimensions or
+			 * dimensions with just a single value here, so let's simply count
+			 * them as 1. They will not affect the volume anyway.
+			 */
+			if (histogram->nvalues[i] > 1)
+				d = (histogram->nvalues[i] - 1);
+
+			bucket_volume *= (double) (maxidx - minidx + 1) / d;
+
+			if (i == 0)
+				format = "{%s"; /* fist dimension */
+			else if (i < (histogram->ndimensions - 1))
+				format = ", %s";	/* medium dimensions */
+			else
+				format = ", %s}";	/* last dimension */
+
+			appendStringInfo(&bufs[3], format, bucket->nullsonly[i] ? "t" : "f");
+			appendStringInfo(&bufs[4], format, bucket->min_inclusive[i] ? "t" : "f");
+			appendStringInfo(&bufs[5], format, bucket->max_inclusive[i] ? "t" : "f");
+
+			/*
+			 * for NULL-only  dimension, simply put there the NULL and
+			 * continue
+			 */
+			if (bucket->nullsonly[i])
+			{
+				if (i == 0)
+					format = "{%s";
+				else if (i < (histogram->ndimensions - 1))
+					format = ", %s";
+				else
+					format = ", %s}";
+
+				appendStringInfo(&bufs[1], format, "NULL");
+				appendStringInfo(&bufs[2], format, "NULL");
+
+				continue;
+			}
+
+			/* otherwise we really need to format the value */
+			switch (otype)
+			{
+				case OUTPUT_FORMAT_RAW: /* actual boundary values */
+
+					if (i == 0)
+						format = "{%s";
+					else if (i < (histogram->ndimensions - 1))
+						format = ", %s";
+					else
+						format = ", %s}";
+
+					appendStringInfo(&bufs[1], format,
+									 FunctionCall1(&fmgrinfo[i], vals[minidx]));
+
+					appendStringInfo(&bufs[2], format,
+									 FunctionCall1(&fmgrinfo[i], vals[maxidx]));
+
+					break;
+
+				case OUTPUT_FORMAT_INDEXES: /* indexes into deduplicated
+											 * arrays */
+
+					if (i == 0)
+						format = "{%d";
+					else if (i < (histogram->ndimensions - 1))
+						format = ", %d";
+					else
+						format = ", %d}";
+
+					appendStringInfo(&bufs[1], format, minidx);
+					appendStringInfo(&bufs[2], format, maxidx);
+
+					break;
+
+				case OUTPUT_FORMAT_DISTINCT:	/* distinct arrays as measure */
+
+					if (i == 0)
+						format = "{%f";
+					else if (i < (histogram->ndimensions - 1))
+						format = ", %f";
+					else
+						format = ", %f}";
+
+					appendStringInfo(&bufs[1], format, (minidx * 1.0 / d));
+					appendStringInfo(&bufs[2], format, (maxidx * 1.0 / d));
+
+					break;
+
+				default:
+					elog(ERROR, "unknown output type: %d", otype);
+			}
+		}
+
+		values[1] = bufs[1].data;
+		values[2] = bufs[2].data;
+		values[3] = bufs[3].data;
+		values[4] = bufs[4].data;
+		values[5] = bufs[5].data;
+
+		snprintf(values[6], 64, "%f", bucket->frequency);	/* frequency */
+		snprintf(values[7], 64, "%f", bucket->frequency / bucket_volume);	/* density */
+		snprintf(values[8], 64, "%f", bucket_volume);	/* volume (as a
+														 * fraction) */
+
+		/* build a tuple */
+		tuple = BuildTupleFromCStrings(attinmeta, values);
+
+		/* make the tuple into a datum */
+		result = HeapTupleGetDatum(tuple);
+
+		/* clean up (this is not really necessary) */
+		pfree(values[0]);
+		pfree(values[6]);
+		pfree(values[7]);
+		pfree(values[8]);
+
+		resetStringInfo(&bufs[1]);
+		resetStringInfo(&bufs[2]);
+		resetStringInfo(&bufs[3]);
+		resetStringInfo(&bufs[4]);
+		resetStringInfo(&bufs[5]);
+
+		pfree(bufs);
+		pfree(values);
+
+		SRF_RETURN_NEXT(funcctx, result);
+	}
+	else						/* do when there is no more left */
+	{
+		SRF_RETURN_DONE(funcctx);
+	}
+}
+
+/*
+ * pg_histogram_in		- input routine for type pg_histogram.
+ *
+ * pg_histogram is real enough to be a table column, but it has no operations
+ * of its own, and disallows input too
+ */
+Datum
+pg_histogram_in(PG_FUNCTION_ARGS)
+{
+	/*
+	 * pg_histogram stores the data in binary form and parsing text input is
+	 * not needed, so disallow this.
+	 */
+	ereport(ERROR,
+			(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
+			 errmsg("cannot accept a value of type %s", "pg_histogram")));
+
+	PG_RETURN_VOID();			/* keep compiler quiet */
+}
+
+/*
+ * pg_histogram_out		- output routine for type pg_histogram.
+ *
+ * histograms are serialized into a bytea value, so we simply call byteaout()
+ * to serialize the value into text. But it'd be nice to serialize that into
+ * a meaningful representation (e.g. for inspection by people).
+ *
+ * XXX This should probably return something meaningful, similar to what
+ * pg_dependencies_out does. Not sure how to deal with the deduplicated
+ * values, though - do we want to expand that or not?
+ */
+Datum
+pg_histogram_out(PG_FUNCTION_ARGS)
+{
+	return byteaout(fcinfo);
+}
+
+/*
+ * pg_histogram_recv		- binary input routine for type pg_histogram.
+ */
+Datum
+pg_histogram_recv(PG_FUNCTION_ARGS)
+{
+	ereport(ERROR,
+			(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
+			 errmsg("cannot accept a value of type %s", "pg_histogram")));
+
+	PG_RETURN_VOID();			/* keep compiler quiet */
+}
+
+/*
+ * pg_histogram_send		- binary output routine for type pg_histogram.
+ *
+ * Histograms are serialized in a bytea value (although the type is named
+ * differently), so let's just send that.
+ */
+Datum
+pg_histogram_send(PG_FUNCTION_ARGS)
+{
+	return byteasend(fcinfo);
+}
+
+/*
+ * selectivity estimation
+ */
+
+/*
+ * When evaluating conditions on the histogram, we can leverage the fact that
+ * each bucket boundary value is used by many buckets (each bucket split
+ * introduces a single new value, duplicating all the other values). That
+ * allows us to significantly reduce the number of function calls by caching
+ * the results.
+ *
+ * This is one of the reasons why we keep the histogram in partially serialized
+ * form, with deduplicated values. This allows us to maintain a simple array
+ * of results indexed by uint16 values.
+ *
+ * We only need 2 bits per value, but we allocate a full char as it's more
+ * convenient and there's not much to gain. 0 means 'unknown' as the function
+ * was not executed for this value yet.
+ */
+
+#define HIST_CACHE_FALSE			0x01
+#define HIST_CACHE_TRUE				0x03
+#define HIST_CACHE_MASK				0x02
+
+/*
+ * bucket_contains_value
+ *		Decide if the bucket (a range of values in a particular dimension) may
+ *		contain the supplied value.
+ *
+ * The function does not simply return true/false, but a "match level" (none,
+ * partial, full), just like other similar functions. In fact, thise function
+ * only returns "partial" or "none" levels, as a range can never match exactly
+ * a value (we never generate histograms with "collapsed" dimensions).
+ *
+ * FIXME Should use a better estimate than DEFAULT_EQ_SEL, e.g. derived
+ * from ndistinct for the variable. But for histograms we shouldn't really
+ * get here, because equalities are handled as conditions (i.e. we'll get
+ * here when deciding which buckets match the conditions, but the fraction
+ * value does not really matter, we only care about the match flag).
+ */
+static bool
+bucket_contains_value(FmgrInfo ltproc, Datum constvalue,
+					  Datum min_value, Datum max_value,
+					  int min_index, int max_index,
+					  bool min_include, bool max_include,
+					  char *callcache, double *fraction)
+{
+	bool		a,
+				b;
+
+	char		min_cached = callcache[min_index];
+	char		max_cached = callcache[max_index];
+
+	/*
+	 * First some quick checks on equality - if any of the boundaries equals,
+	 * we have a partial match (so no need to call the comparator).
+	 */
+	if (((min_value == constvalue) && (min_include)) ||
+		((max_value == constvalue) && (max_include)))
+	{
+		*fraction = DEFAULT_EQ_SEL;
+		return true;
+	}
+
+	/* Keep the values 0/1 because of the XOR at the end. */
+	a = ((min_cached & HIST_CACHE_MASK) >> 1);
+	b = ((max_cached & HIST_CACHE_MASK) >> 1);
+
+	/*
+	 * If result for the bucket lower bound not in cache, evaluate the
+	 * function and store the result in the cache.
+	 */
+	if (!min_cached)
+	{
+		a = DatumGetBool(FunctionCall2Coll(&ltproc,
+										   DEFAULT_COLLATION_OID,
+										   constvalue, min_value));
+		/* remember the result */
+		callcache[min_index] = (a) ? HIST_CACHE_TRUE : HIST_CACHE_FALSE;
+	}
+
+	/* And do the same for the upper bound. */
+	if (!max_cached)
+	{
+		b = DatumGetBool(FunctionCall2Coll(&ltproc,
+										   DEFAULT_COLLATION_OID,
+										   constvalue, max_value));
+		/* remember the result */
+		callcache[max_index] = (b) ? HIST_CACHE_TRUE : HIST_CACHE_FALSE;
+	}
+
+	*fraction = (a ^ b) ? DEFAULT_EQ_SEL : 0.0;
+
+	return (a ^ b) ? true : false;
+}
+
+/*
+ * bucket_is_smaller_than_value
+ *		Decide if the bucket (a range of values in a particular dimension) is
+ *		smaller than the supplied value.
+ *
+ * The function does not simply return true/false, but a "match level" (none,
+ * partial, full), just like other similar functions.
+ *
+ * Unlike bucket_contains_value this may return all three match levels, i.e.
+ * "full" (e.g. [10,20] < 30), "partial" (e.g. [10,20] < 15) and "none"
+ * (e.g. [10,20] < 5).
+ *
+ * FIXME Use a better estimate, instead of DEFAULT_INEQ_SEL, i.e. something
+ * derived in a way similar to convert_to_scalar.
+ */
+static bool
+bucket_is_smaller_than_value(FmgrInfo opproc, Oid typeoid, Datum constvalue,
+							 Datum min_value, Datum max_value,
+							 int min_index, int max_index,
+							 bool min_include, bool max_include,
+							 char *callcache, bool isgt,
+							 double *fraction)
+{
+	char		min_cached = callcache[min_index];
+	char		max_cached = callcache[max_index];
+
+	/* Keep the values 0/1 because of the XOR at the end. */
+	bool		a = ((min_cached & HIST_CACHE_MASK) >> 1);
+	bool		b = ((max_cached & HIST_CACHE_MASK) >> 1);
+
+	if (!min_cached)
+	{
+		a = DatumGetBool(FunctionCall2Coll(&opproc,
+										   DEFAULT_COLLATION_OID,
+										   min_value,
+										   constvalue));
+		/* remember the result */
+		callcache[min_index] = (a) ? HIST_CACHE_TRUE : HIST_CACHE_FALSE;
+	}
+
+	if (!max_cached)
+	{
+		b = DatumGetBool(FunctionCall2Coll(&opproc,
+										   DEFAULT_COLLATION_OID,
+										   max_value,
+										   constvalue));
+		/* remember the result */
+		callcache[max_index] = (b) ? HIST_CACHE_TRUE : HIST_CACHE_FALSE;
+	}
+
+	/*
+	 * Now, we need to combine both results into the final answer, and we need
+	 * to be careful about the 'isgt' variable which kinda inverts the
+	 * meaning.
+	 *
+	 * First, we handle the case when each boundary returns different results.
+	 * In that case the outcome can only be 'partial' match, and the fraction
+	 * is computed using convert_to_scalar, just like for 1D histograms.
+	 */
+	if (a != b)
+	{
+		double	val, high, low, binfrac;
+
+		if (convert_to_scalar(constvalue, typeoid, &val, min_value, max_value,
+							  typeoid, &low, &high))
+		{
+
+			/* shamelessly copied from ineq_histogram_selectivity */
+			if (high <= low)
+			{
+				/* cope if bin boundaries appear identical */
+				binfrac = 0.5;
+			}
+			else if (val <= low)
+				binfrac = 0.0;
+			else if (val >= high)
+				binfrac = 1.0;
+			else
+			{
+				binfrac = (val - low) / (high - low);
+
+				/*
+				 * Watch out for the possibility that we got a NaN or
+				 * Infinity from the division.  This can happen
+				 * despite the previous checks, if for example "low"
+				 * is -Infinity.
+				 */
+				if (isnan(binfrac) ||
+					binfrac < 0.0 || binfrac > 1.0)
+					binfrac = 0.5;
+			}
+		}
+		else
+			binfrac = 0.5;
+
+		*fraction = (isgt) ? binfrac : (1-binfrac);
+		return true;
+	}
+
+	/*
+	 * When the results are the same, then it depends on the 'isgt' value.
+	 * There are four options:
+	 *
+	 * isgt=false a=b=true	=> full match isgt=false a=b=false => empty
+	 * isgt=true  a=b=true	=> empty isgt=true	a=b=false => full match
+	 *
+	 * We'll cheat a bit, because we know that (a=b) so we'll use just one of
+	 * them.
+	 */
+	if (isgt)
+	{
+		*fraction = (!a) ? 1.0 : 0.0;
+		return (!a);
+	}
+	else
+	{
+		*fraction = (a) ? 1.0 : 0.0;
+		return a;
+	}
+}
+
+/*
+ * Evaluate clauses using the histogram, and update the match bitmap.
+ *
+ * The bitmap may be already partially set, so this is really a way to
+ * combine results of several clause lists - either when computing
+ * conditional probability P(A|B) or a combination of AND/OR clauses.
+ *
+ * Note: This is not a simple bitmap in the sense that there are three
+ * possible values for each item - no match, partial match and full match.
+ * So we need at least 2 bits per item.
+ *
+ * TODO: This works with 'bitmap' where each item is represented as a
+ * char, which is slightly wasteful. Instead, we could use a bitmap
+ * with 2 bits per item, reducing the size to ~1/4. By using values
+ * 0, 1 and 3 (instead of 0, 1 and 2), the operations (merging etc.)
+ * might be performed just like for simple bitmap by using & and |,
+ * which might be faster than min/max.
+ */
+static void
+histogram_update_match_bitmap(PlannerInfo *root, List *clauses,
+							  Bitmapset *stakeys,
+							  MVHistogram * histogram,
+							  bucket_match *matches, bool is_or)
+{
+	int			i;
+	ListCell   *l;
+
+	/*
+	 * Used for caching function calls, only once per deduplicated value.
+	 *
+	 * We know may have up to (2 * nbuckets) values per dimension. It's
+	 * probably overkill, but let's allocate that once for all clauses, to
+	 * minimize overhead.
+	 *
+	 * Also, we only need two bits per value, but this allocates byte per
+	 * value. Might be worth optimizing.
+	 *
+	 * 0x00 - not yet called 0x01 - called, result is 'false' 0x03 - called,
+	 * result is 'true'
+	 */
+	char	   *callcache = palloc(histogram->nbuckets);
+
+	Assert(histogram != NULL);
+	Assert(histogram->nbuckets > 0);
+
+	Assert(clauses != NIL);
+	Assert(list_length(clauses) >= 1);
+
+	/* loop through the clauses and do the estimation */
+	foreach(l, clauses)
+	{
+		Node	   *clause = (Node *) lfirst(l);
+
+		/* if it's a RestrictInfo, then extract the clause */
+		if (IsA(clause, RestrictInfo))
+			clause = (Node *) ((RestrictInfo *) clause)->clause;
+
+		/* it's either OpClause, or NullTest */
+		if (is_opclause(clause))
+		{
+			OpExpr	   *expr = (OpExpr *) clause;
+			bool		varonleft = true;
+			bool		ok;
+
+			FmgrInfo	opproc; /* operator */
+
+			fmgr_info(get_opcode(expr->opno), &opproc);
+
+			/* reset the cache (per clause) */
+			memset(callcache, 0, histogram->nbuckets);
+
+			ok = (NumRelids(clause) == 1) &&
+				(is_pseudo_constant_clause(lsecond(expr->args)) ||
+				 (varonleft = false,
+				  is_pseudo_constant_clause(linitial(expr->args))));
+
+			if (ok)
+			{
+				FmgrInfo	ltproc;
+				RegProcedure oprrest = get_oprrest(expr->opno);
+
+				Var		   *var = (varonleft) ? linitial(expr->args) : lsecond(expr->args);
+				Const	   *cst = (varonleft) ? lsecond(expr->args) : linitial(expr->args);
+				bool		isgt = (!varonleft);
+
+				TypeCacheEntry *typecache
+				= lookup_type_cache(var->vartype, TYPECACHE_LT_OPR);
+
+				/* lookup dimension for the attribute */
+				int			idx = bms_member_index(stakeys, var->varattno);
+
+				fmgr_info(get_opcode(typecache->lt_opr), &ltproc);
+
+				/*
+				 * Check this for all buckets that still have "true" in the
+				 * bitmap
+				 *
+				 * We already know the clauses use suitable operators (because
+				 * that's how we filtered them).
+				 */
+				for (i = 0; i < histogram->nbuckets; i++)
+				{
+					bool		res;
+					double		fraction;
+
+					MVBucket   *bucket = histogram->buckets[i];
+
+					/* histogram boundaries */
+					Datum		minval,
+								maxval;
+					bool		mininclude,
+								maxinclude;
+					int			minidx,
+								maxidx;
+
+					/*
+					 * For AND-lists, we can also mark NULL buckets as 'no
+					 * match' (and then skip them). For OR-lists this is not
+					 * possible.
+					 */
+					if ((!is_or) && bucket->nullsonly[idx])
+						matches[i].match = false;
+
+					/*
+					 * XXX There used to be logic to skip buckets that can't
+					 * possibly match, depending on the is_or flag (either
+					 * fully matching or elimated). Once we abandoned the
+					 * concept of NONE/PARTIAL/FULL matches and switched to
+					 * a bool flag + fraction that does not seem possible.
+					 * But maybe we can make it work somehow?
+					 */
+
+					/* lookup the values and cache of function calls */
+					minidx = bucket->min[idx];
+					maxidx = bucket->max[idx];
+
+					minval = histogram->values[idx][bucket->min[idx]];
+					maxval = histogram->values[idx][bucket->max[idx]];
+
+					mininclude = bucket->min_inclusive[idx];
+					maxinclude = bucket->max_inclusive[idx];
+
+					/*
+					 * If it's not a "<" or ">" or "=" operator, just ignore
+					 * the clause. Otherwise note the relid and attnum for the
+					 * variable.
+					 *
+					 * TODO I'm really unsure the handling of 'isgt' flag
+					 * (that is, clauses with reverse order of
+					 * variable/constant) is correct. I wouldn't be surprised
+					 * if there was some mixup. Using the lt/gt operators
+					 * instead of messing with the opproc could make it
+					 * simpler. It would however be using a different operator
+					 * than the query, although it's not any shadier than
+					 * using the selectivity function as is done currently.
+					 */
+					switch (oprrest)
+					{
+						case F_SCALARLTSEL: /* Var < Const */
+						case F_SCALARLESEL: /* Var <= Const */
+						case F_SCALARGTSEL: /* Var > Const */
+						case F_SCALARGESEL: /* Var >= Const */
+
+							res = bucket_is_smaller_than_value(opproc, var->vartype,
+															   cst->constvalue,
+															   minval, maxval,
+															   minidx, maxidx,
+															   mininclude, maxinclude,
+															   callcache, isgt, &fraction);
+
+							break;
+
+						case F_EQSEL:
+						case F_NEQSEL:
+
+							/*
+							 * We only check whether the value is within the
+							 * bucket, using the lt operator, and we also
+							 * check for equality with the boundaries.
+							 */
+
+							res = bucket_contains_value(ltproc, cst->constvalue,
+														minval, maxval,
+														minidx, maxidx,
+														mininclude, maxinclude,
+														callcache, &fraction);
+
+							break;
+					}
+
+					/*
+					 * Merge the result into the bitmap, depending on type
+					 * of the current clause (AND or OR).
+					 */
+					if (is_or)
+					{
+						Selectivity s1, s2;
+
+						/* OR follows the Max() semantics */
+						matches[i].match |= res;
+
+						/*
+						 * Selectivities for an OR clause are combined as s1+s2 - s1*s2
+						 * to account for the probable overlap of selected tuple sets.
+						 * This is the same formula as in clause_selectivity, because
+						 * the fraction is computed assuming independence (but then we
+						 * also apply geometric mean).
+						 */
+						s1 = matches[i].fraction;
+						s2 = fraction;
+
+						matches[i].fraction = s1 + s2 - s1 * s2;
+
+						CLAMP_PROBABILITY(matches[i].fraction);
+					}
+					else
+					{
+						/* AND follows Min() semantics */
+						matches[i].match &= res;
+						matches[i].fraction *= fraction;
+					}
+				}
+			}
+		}
+		else if (IsA(clause, NullTest))
+		{
+			NullTest   *expr = (NullTest *) clause;
+			Var		   *var = (Var *) (expr->arg);
+
+			/* lookup index of attribute in the statistics */
+			int			idx = bms_member_index(stakeys, var->varattno);
+
+			/*
+			 * Walk through the buckets and evaluate the current clause. We
+			 * can skip items that were already ruled out, and terminate if
+			 * there are no remaining buckets that might possibly match.
+			 */
+			for (i = 0; i < histogram->nbuckets; i++)
+			{
+				char		match = false;
+				MVBucket   *bucket = histogram->buckets[i];
+
+				/*
+				 * Skip buckets that were already eliminated - this is
+				 * impotant considering how we update the info (we only lower
+				 * the match)
+				 */
+				if ((!is_or) && (!matches[i].match))
+					continue;
+				else if (is_or && (matches[i].match))
+					continue;
+
+				switch (expr->nulltesttype)
+				{
+					case IS_NULL:
+						match = (bucket->nullsonly[idx]) ? true : match;
+						break;
+
+					case IS_NOT_NULL:
+						match = (!bucket->nullsonly[idx]) ? true : match;
+						break;
+				}
+
+				/* now, update the match bitmap, depending on OR/AND type */
+				if (is_or)
+				{
+					matches[i].match |= match;
+					matches[i].fraction = (match) ? 1.0 : matches[i].fraction;
+				}
+				else
+				{
+					matches[i].match &= match;
+					matches[i].fraction = (match) ? matches[i].fraction : 0.0;
+				}
+			}
+		}
+		else if (or_clause(clause) || and_clause(clause))
+		{
+			/*
+			 * AND/OR clause, with all sub-clauses compatible with the stats
+			 */
+
+			int			i;
+			BoolExpr   *bool_clause = ((BoolExpr *) clause);
+			List	   *bool_clauses = bool_clause->args;
+
+			/* match/mismatch bitmap for each bucket */
+			bucket_match   *bool_matches = NULL;
+
+			Assert(bool_clauses != NIL);
+			Assert(list_length(bool_clauses) >= 2);
+
+			/* by default none of the buckets matches the clauses */
+			bool_matches = palloc0(sizeof(bucket_match) * histogram->nbuckets);
+
+			if (or_clause(clause))
+			{
+				/* OR clauses assume nothing matches, initially */
+				for (i = 0; i < histogram->nbuckets; i++)
+				{
+					bool_matches[i].match = false;
+					bool_matches[i].fraction = 0.0;
+				}
+			}
+			else
+			{
+				/* AND clauses assume nothing matches, initially */
+				for (i = 0; i < histogram->nbuckets; i++)
+				{
+					bool_matches[i].match = true;
+					bool_matches[i].fraction = 1.0;
+				}
+			}
+
+			/* build the match bitmap for the OR-clauses */
+			histogram_update_match_bitmap(root, bool_clauses,
+										  stakeys, histogram,
+										  bool_matches, or_clause(clause));
+
+			/*
+			 * Merge the bitmap produced by histogram_update_match_bitmap into
+			 * the current one. We need to consider if we're evaluating AND or
+			 * OR condition when merging the results.
+			 */
+			for (i = 0; i < histogram->nbuckets; i++)
+			{
+				/* Is this OR or AND clause? */
+				if (is_or)
+				{
+					Selectivity	s1, s2;
+
+					matches[i].match |= bool_matches[i].match;
+
+					/*
+					 * Selectivities for an OR clause are combined as s1+s2 - s1*s2
+					 * to account for the probable overlap of selected tuple sets.
+					 * This is the same formula as in clause_selectivity, because
+					 * the fraction is computed assuming independence (but then we
+					 * also apply geometric mean).
+					 */
+					s1 = matches[i].fraction;
+					s2 = bool_matches[i].fraction;
+
+					matches[i].fraction = s1 + s2 - s1 * s2;
+
+					CLAMP_PROBABILITY(matches[i].fraction);
+				}
+				else
+				{
+					matches[i].match &= bool_matches[i].match;
+					matches[i].fraction *= bool_matches[i].fraction;
+				}
+			}
+
+			pfree(bool_matches);
+
+		}
+		else if (not_clause(clause))
+		{
+			/* NOT clause, with all subclauses compatible */
+
+			int			i;
+			BoolExpr   *not_clause = ((BoolExpr *) clause);
+			List	   *not_args = not_clause->args;
+
+			/* match/mismatch bitmap for each MCV item */
+			bucket_match   *not_matches = NULL;
+
+			Assert(not_args != NIL);
+			Assert(list_length(not_args) == 1);
+
+			/* by default none of the MCV items matches the clauses */
+			not_matches = palloc0(sizeof(bucket_match) * histogram->nbuckets);
+
+			/* NOT clauses assume nothing matches, initially
+			 *
+			 * FIXME The comment seems to disagree with the code - not sure
+			 * if nothing should match (code is wrong) or everything should
+			 * match (comment is wrong) by default.
+			 */
+			for (i = 0; i < histogram->nbuckets; i++)
+			{
+				not_matches[i].match = true;
+				not_matches[i].fraction = 1.0;
+			}
+
+			/* build the match bitmap for the OR-clauses */
+			histogram_update_match_bitmap(root, not_args,
+										  stakeys, histogram,
+										  not_matches, false);
+
+			/*
+			 * Merge the bitmap produced by histogram_update_match_bitmap into
+			 * the current one.
+			 *
+			 * This is similar to what mcv_update_match_bitmap does, but we
+			 * need to be a tad more careful here, as histograms also track
+			 * what fraction of a bucket matches.
+			 */
+			for (i = 0; i < histogram->nbuckets; i++)
+			{
+				/*
+				 * When handling a NOT clause, invert the result before
+				 * merging it into the global result. We don't care about
+				 * partial matches here (those invert to partial).
+				 */
+				not_matches[i].match = (!not_matches[i].match);
+
+				/* Is this OR or AND clause? */
+				if (is_or)
+				{
+					Selectivity s1, s2;
+
+					matches[i].match |= not_matches[i].match;
+
+					/*
+					 * Selectivities for an OR clause are combined as s1+s2 - s1*s2
+					 * to account for the probable overlap of selected tuple sets.
+					 * This is the same formula as in clause_selectivity, because
+					 * the fraction is computed assuming independence (but then we
+					 * also apply geometric mean).
+					 */
+					s1 = matches[i].fraction;
+					s2 = not_matches[i].fraction;
+
+					matches[i].fraction = s1 + s2 - s1 * s2;
+
+					CLAMP_PROBABILITY(matches[i].fraction);
+				}
+				else
+				{
+					matches[i].match &= not_matches[i].match;
+					matches[i].fraction *= not_matches[i].fraction;
+				}
+			}
+
+			pfree(not_matches);
+		}
+		else if (IsA(clause, Var))
+		{
+			/* Var (has to be a boolean Var, possibly from below NOT) */
+
+			Var		   *var = (Var *) (clause);
+
+			/* match the attribute to a dimension of the statistic */
+			int			idx = bms_member_index(stakeys, var->varattno);
+
+			Assert(var->vartype == BOOLOID);
+
+			/*
+			 * Walk through the buckets and evaluate the current clause.
+			 */
+			for (i = 0; i < histogram->nbuckets; i++)
+			{
+				MVBucket   *bucket = histogram->buckets[i];
+				bool	match = false;
+				double	fraction = 0.0;
+
+				/*
+				 * If the bucket is NULL, it's a mismatch. Otherwise check
+				 * if lower/upper boundaries match and choose partial/full
+				 * match accordingly.
+				 */
+				if (!bucket->nullsonly[idx])
+				{
+					int		minidx = bucket->min[idx];
+					int		maxidx = bucket->max[idx];
+
+					bool 	a = DatumGetBool(histogram->values[idx][minidx]);
+					bool	b = DatumGetBool(histogram->values[idx][maxidx]);
+
+					/* How many boundary values match? */
+					if (a && b)
+					{
+						/* both values match - the whole bucket matches */
+						match = true;
+						fraction = 1.0;
+					}
+					else if (a || b)
+					{
+						/* one value matches - assume half the bucket matches */
+						match = true;
+						fraction = 0.5;
+					}
+				}
+
+				/* now, update the match bitmap, depending on OR/AND type */
+				if (is_or)
+				{
+					Selectivity	s1, s2;
+
+					matches[i].match |= match;
+
+					/*
+					 * Selectivities for an OR clause are combined as s1+s2 - s1*s2
+					 * to account for the probable overlap of selected tuple sets.
+					 * This is the same formula as in clause_selectivity, because
+					 * the fraction is computed assuming independence (but then we
+					 * also apply geometric mean).
+					 */
+					s1 = matches[i].fraction;
+					s2 = fraction;
+
+					matches[i].fraction = s1 + s2 - s1 * s2;
+
+					CLAMP_PROBABILITY(matches[i].fraction);
+				}
+				else
+				{
+					matches[i].match &= match;
+					matches[i].fraction *= fraction;
+				}
+			}
+		}
+		else
+			elog(ERROR, "unknown clause type: %d", clause->type);
+	}
+
+	/* free the call cache */
+	pfree(callcache);
+}
+
+/*
+ * Estimate selectivity of clauses using a histogram.
+ *
+ * If there's no histogram for the stats, the function returns 0.0.
+ *
+ * The general idea of this method is similar to how MCV lists are
+ * processed, except that this introduces the concept of a partial
+ * match (MCV only works with full match / mismatch).
+ *
+ * The algorithm works like this:
+ *
+ *	 1) mark all buckets as 'full match'
+ *	 2) walk through all the clauses
+ *	 3) for a particular clause, walk through all the buckets
+ *	 4) skip buckets that are already 'no match'
+ *	 5) check clause for buckets that still match (at least partially)
+ *	 6) sum frequencies for buckets to get selectivity
+ *
+ * Unlike MCV lists, histograms have a concept of a partial match. In
+ * that case we use 1/2 the bucket, to minimize the average error. The
+ * MV histograms are usually less detailed than the per-column ones,
+ * meaning the sum is often quite high (thanks to combining a lot of
+ * "partially hit" buckets).
+ *
+ * Maybe we could use per-bucket information with number of distinct
+ * values it contains (for each dimension), and then use that to correct
+ * the estimate (so with 10 distinct values, we'd use 1/10 of the bucket
+ * frequency). We might also scale the value depending on the actual
+ * ndistinct estimate (not just the values observed in the sample).
+ *
+ * Another option would be to multiply the selectivities, i.e. if we get
+ * 'partial match' for a bucket for multiple conditions, we might use
+ * 0.5^k (where k is the number of conditions), instead of 0.5. This
+ * probably does not minimize the average error, though.
+ *
+ * TODO: This might use a similar shortcut to MCV lists - count buckets
+ * marked as partial/full match, and terminate once this drop to 0.
+ * Not sure if it's really worth it - for MCV lists a situation like
+ * this is not uncommon, but for histograms it's not that clear.
+ */
+Selectivity
+histogram_clauselist_selectivity(PlannerInfo *root, StatisticExtInfo *stat,
+								 List *clauses, List *conditions,
+								 int varRelid, JoinType jointype,
+								 SpecialJoinInfo *sjinfo, RelOptInfo *rel)
+{
+	int			i;
+	MVHistogram *histogram;
+	Selectivity	s = 0.0;
+	Selectivity	total_sel = 0.0;
+	Size		len;
+	int			nclauses;
+
+	/* match/mismatch bitmap for each MCV item */
+	bucket_match   *matches = NULL;
+	bucket_match   *condition_matches = NULL;
+
+	nclauses = list_length(clauses);
+
+	/* load the histogram stored in the statistics object */
+	histogram = statext_histogram_load(stat->statOid);
+
+	/* size of the match "bitmap" */
+	len = sizeof(bucket_match) * histogram->nbuckets;
+
+	/* by default all the histogram buckets match the clauses fully */
+	matches = palloc0(len);
+
+	/* by default all buckets match fully */
+	for (i = 0; i < histogram->nbuckets; i++)
+	{
+		matches[i].match = true;
+		matches[i].fraction = 1.0;
+	}
+
+	histogram_update_match_bitmap(root, clauses, stat->keys,
+								  histogram, matches, false);
+
+	/* if there are condition clauses, build a match bitmap for them */
+	if (conditions)
+	{
+		/* match bitmap for conditions, by default all buckets match */
+		condition_matches = palloc0(len);
+
+		/* by default all buckets match fully */
+		for (i = 0; i < histogram->nbuckets; i++)
+		{
+			condition_matches[i].match = true;
+			condition_matches[i].fraction = 1.0;
+		}
+
+		histogram_update_match_bitmap(root, conditions, stat->keys,
+									  histogram, condition_matches, false);
+	}
+
+	/* now, walk through the buckets and sum the selectivities */
+	for (i = 0; i < histogram->nbuckets; i++)
+	{
+		double fraction;
+
+		/* skip buckets that don't satisfy the conditions */
+		if (conditions && (!condition_matches[i].match))
+			continue;
+
+		/* compute selectivity for buckets matching conditions */
+		total_sel += histogram->buckets[i]->frequency;
+
+		/* geometric mean of the bucket fraction */
+		fraction = pow(matches[i].fraction, 1.0 / nclauses);
+
+		if (matches[i].match)
+			s += histogram->buckets[i]->frequency * fraction;
+	}
+
+	/* conditional selectivity P(clauses|conditions) */
+	if (total_sel > 0.0)
+		return (s / total_sel);
+
+	return 0.0;
+}
diff --git a/src/backend/statistics/mcv.c b/src/backend/statistics/mcv.c
index 533fbdc037..79e6f24deb 100644
--- a/src/backend/statistics/mcv.c
+++ b/src/backend/statistics/mcv.c
@@ -85,7 +85,8 @@ static int count_distinct_groups(int numrows, SortItem *items,
  */
 MCVList *
 statext_mcv_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
-				  VacAttrStats **stats, double totalrows)
+				  VacAttrStats **stats, HeapTuple **rows_filtered,
+				  int *numrows_filtered, double totalrows)
 {
 	int			i,
 				j,
@@ -96,6 +97,7 @@ statext_mcv_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
 	double		stadistinct;
 	int		   *mcv_counts;
 	int			f1;
+	int			numrows_mcv;
 
 	int		   *attnums = build_attnums(attrs);
 
@@ -111,6 +113,9 @@ statext_mcv_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
 	/* transform the sorted rows into groups (sorted by frequency) */
 	SortItem   *groups = build_distinct_groups(numrows, items, mss, &ngroups);
 
+	/* Either we have both pointers or none of them. */
+	Assert((rows_filtered && numrows_filtered) || (!rows_filtered && !numrows_filtered));
+
 	/*
 	 * Maximum number of MCV items to store, based on the attribute with the
 	 * largest stats target (and the number of groups we have available).
@@ -167,6 +172,9 @@ statext_mcv_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
 								  numrows, totalrows);
 	}
 
+	/* number of rows represented by MCV items */
+	numrows_mcv = 0;
+
 	/*
 	 * At this point we know the number of items for the MCV list. There might
 	 * be none (for uniform distribution with many groups), and in that case
@@ -243,9 +251,93 @@ statext_mcv_build(int numrows, HeapTuple *rows, Bitmapset *attrs,
 
 				item->base_frequency *= (double) count / numrows;
 			}
+
+			/* update the number of sampled rows represented by the MCV list */
+			numrows_mcv += groups[i].count;
 		}
 	}
 
+	/* Assume we're not returning any filtered rows by default. */
+	if (numrows_filtered)
+		*numrows_filtered = 0;
+
+	if (rows_filtered)
+		*rows_filtered = NULL;
+
+	/*
+	 * Produce an array with only tuples not covered by the MCV list. This is
+	 * needed when building MCV+histogram pair, where MCV covers the most
+	 * common combinations and histogram covers the remaining part.
+	 *
+	 * We will first sort the groups by the keys (not by count) and then use
+	 * binary search in the group array to check which rows are covered by the
+	 * MCV items.
+	 *
+	 * Do not modify the array in place, as there may be additional stats on
+	 * the table and we need to keep the original array for them.
+	 *
+	 * We only do this when requested by passing non-NULL rows_filtered, and
+	 * when there are rows not covered by the MCV list (that is, when
+	 * numrows_mcv < numrows), or also (nitems < ngroups).
+	 */
+	if (rows_filtered && numrows_filtered && (nitems < ngroups))
+	{
+		int			i,
+					j;
+
+		/* used to build the filtered array of tuples */
+		HeapTuple  *filtered;
+		int			nfiltered;
+
+		/* used for the searches */
+		SortItem	key;
+
+		/* We do know how many rows we expect (total - MCV rows). */
+		nfiltered = (numrows - numrows_mcv);
+		filtered = (HeapTuple *) palloc(nfiltered * sizeof(HeapTuple));
+
+		/* wfill this with data from the rows */
+		key.values = (Datum *) palloc0(numattrs * sizeof(Datum));
+		key.isnull = (bool *) palloc0(numattrs * sizeof(bool));
+
+		/*
+		 * Sort the groups for bsearch_r (but only the items that actually
+		 * made it to the MCV list).
+		 */
+		qsort_arg((void *) groups, nitems, sizeof(SortItem),
+				  multi_sort_compare, mss);
+
+		/* walk through the tuples, compare the values to MCV items */
+		nfiltered = 0;
+		for (i = 0; i < numrows; i++)
+		{
+			/* collect the key values from the row */
+			for (j = 0; j < numattrs; j++)
+				key.values[j]
+					= heap_getattr(rows[i], attnums[j],
+								   stats[j]->tupDesc, &key.isnull[j]);
+
+			/* if not included in the MCV list, keep it in the array */
+			if (bsearch_arg(&key, groups, nitems, sizeof(SortItem),
+							multi_sort_compare, mss) == NULL)
+				filtered[nfiltered++] = rows[i];
+
+			/* do not overflow the array */
+			Assert(nfiltered <= (numrows - numrows_mcv));
+		}
+
+		/* expect to get the right number of remaining rows exactly */
+		Assert(nfiltered + numrows_mcv == numrows);
+
+		/* pass the filtered tuples up */
+		*numrows_filtered = nfiltered;
+		*rows_filtered = filtered;
+
+		/* free all the data used here */
+		pfree(key.values);
+		pfree(key.isnull);
+	}
+
 	pfree(items);
 	pfree(groups);
 	pfree(mcv_counts);
diff --git a/src/backend/utils/adt/ruleutils.c b/src/backend/utils/adt/ruleutils.c
index c941c3310b..837660950e 100644
--- a/src/backend/utils/adt/ruleutils.c
+++ b/src/backend/utils/adt/ruleutils.c
@@ -1505,6 +1505,7 @@ pg_get_statisticsobj_worker(Oid statextid, bool missing_ok)
 	bool		ndistinct_enabled;
 	bool		dependencies_enabled;
 	bool		mcv_enabled;
+	bool		histogram_enabled;
 	int			i;
 
 	statexttup = SearchSysCache1(STATEXTOID, ObjectIdGetDatum(statextid));
@@ -1541,6 +1542,7 @@ pg_get_statisticsobj_worker(Oid statextid, bool missing_ok)
 	ndistinct_enabled = false;
 	dependencies_enabled = false;
 	mcv_enabled = false;
+	histogram_enabled = false;
 
 	for (i = 0; i < ARR_DIMS(arr)[0]; i++)
 	{
@@ -1550,6 +1552,8 @@ pg_get_statisticsobj_worker(Oid statextid, bool missing_ok)
 			dependencies_enabled = true;
 		if (enabled[i] == STATS_EXT_MCV)
 			mcv_enabled = true;
+		if (enabled[i] == STATS_EXT_HISTOGRAM)
+			histogram_enabled = true;
 	}
 
 	/*
@@ -1578,7 +1582,13 @@ pg_get_statisticsobj_worker(Oid statextid, bool missing_ok)
 		}
 
 		if (mcv_enabled)
+		{
 			appendStringInfo(&buf, "%smcv", gotone ? ", " : "");
+			gotone = true;
+		}
+
+		if (histogram_enabled)
+			appendStringInfo(&buf, "%shistogram", gotone ? ", " : "");
 
 		appendStringInfoChar(&buf, ')');
 	}
diff --git a/src/backend/utils/adt/selfuncs.c b/src/backend/utils/adt/selfuncs.c
index fdfc0d6a1b..bd1cbe8e50 100644
--- a/src/backend/utils/adt/selfuncs.c
+++ b/src/backend/utils/adt/selfuncs.c
@@ -172,9 +172,6 @@ static double eqjoinsel_semi(Oid operator,
 			   RelOptInfo *inner_rel);
 static bool estimate_multivariate_ndistinct(PlannerInfo *root,
 								RelOptInfo *rel, List **varinfos, double *ndistinct);
-static bool convert_to_scalar(Datum value, Oid valuetypid, double *scaledvalue,
-				  Datum lobound, Datum hibound, Oid boundstypid,
-				  double *scaledlobound, double *scaledhibound);
 static double convert_numeric_to_scalar(Datum value, Oid typid, bool *failure);
 static void convert_string_to_scalar(char *value,
 						 double *scaledvalue,
@@ -4104,7 +4101,7 @@ estimate_multivariate_ndistinct(PlannerInfo *root, RelOptInfo *rel,
 		int			nshared;
 
 		/* skip statistics of other kinds */
-		if (info->kind != STATS_EXT_NDISTINCT)
+		if ((info->kinds & STATS_EXT_INFO_NDISTINCT) == 0)
 			continue;
 
 		/* compute attnums shared by the vars and the statistics object */
@@ -4213,7 +4210,7 @@ estimate_multivariate_ndistinct(PlannerInfo *root, RelOptInfo *rel,
  * The several datatypes representing relative times (intervals) are all
  * converted to measurements expressed in seconds.
  */
-static bool
+bool
 convert_to_scalar(Datum value, Oid valuetypid, double *scaledvalue,
 				  Datum lobound, Datum hibound, Oid boundstypid,
 				  double *scaledlobound, double *scaledhibound)
diff --git a/src/bin/psql/describe.c b/src/bin/psql/describe.c
index 394301caa8..929aacf906 100644
--- a/src/bin/psql/describe.c
+++ b/src/bin/psql/describe.c
@@ -2519,7 +2519,8 @@ describeOneTableDetails(const char *schemaname,
 							  "        a.attnum = s.attnum AND NOT attisdropped)) AS columns,\n"
 							  "  'd' = any(stxkind) AS ndist_enabled,\n"
 							  "  'f' = any(stxkind) AS deps_enabled,\n"
-							  "  'm' = any(stxkind) AS mcv_enabled\n"
+							  "  'm' = any(stxkind) AS mcv_enabled,\n"
+							  "  'h' = any(stxkind) AS histogram_enabled\n"
 							  "FROM pg_catalog.pg_statistic_ext stat "
 							  "WHERE stxrelid = '%s'\n"
 							  "ORDER BY 1;",
@@ -2562,6 +2563,12 @@ describeOneTableDetails(const char *schemaname,
 					if (strcmp(PQgetvalue(result, i, 7), "t") == 0)
 					{
 						appendPQExpBuffer(&buf, "%smcv", gotone ? ", " : "");
+						gotone = true;
+					}
+
+					if (strcmp(PQgetvalue(result, i, 8), "t") == 0)
+					{
+						appendPQExpBuffer(&buf, "%shistogram", gotone ? ", " : "");
 					}
 
 					appendPQExpBuffer(&buf, ") ON %s FROM %s",
diff --git a/src/include/catalog/pg_cast.dat b/src/include/catalog/pg_cast.dat
index dff3a9a08a..99d250bb14 100644
--- a/src/include/catalog/pg_cast.dat
+++ b/src/include/catalog/pg_cast.dat
@@ -330,6 +330,10 @@
 { castsource => 'pg_mcv_list', casttarget => 'text', castfunc => '0',
   castcontext => 'i', castmethod => 'i' },
 
+# pg_histogram can be coerced to, but not from, bytea
+{ castsource => 'pg_histogram', casttarget => 'bytea', castfunc => '0',
+  castcontext => 'i', castmethod => 'b' },
+
 # Datetime category
 { castsource => 'abstime', casttarget => 'date', castfunc => 'date(abstime)',
   castcontext => 'a', castmethod => 'f' },
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
index 3cfcafcb11..cbe2b36923 100644
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -5097,6 +5097,30 @@
   proargnames => '{mcv_list,index,values,nulls,frequency,base_frequency}',
   prosrc => 'pg_stats_ext_mcvlist_items' },
 
+{ oid => '3426', descr => 'I/O',
+  proname => 'pg_histogram_in', prorettype => 'pg_histogram',
+  proargtypes => 'cstring', prosrc => 'pg_histogram_in' },
+{ oid => '3427', descr => 'I/O',
+  proname => 'pg_histogram_out', prorettype => 'cstring',
+  proargtypes => 'pg_histogram', prosrc => 'pg_histogram_out' },
+{ oid => '3428', descr => 'I/O',
+  proname => 'pg_histogram_recv', provolatile => 's',
+  prorettype => 'pg_histogram', proargtypes => 'internal',
+  prosrc => 'pg_histogram_recv' },
+{ oid => '3429', descr => 'I/O',
+  proname => 'pg_histogram_send', provolatile => 's', prorettype => 'bytea',
+  proargtypes => 'pg_histogram', prosrc => 'pg_histogram_send' },
+
+{ oid => '3430',
+  descr => 'details about histogram buckets',
+  proname => 'pg_histogram_buckets', prorows => '1000', proisstrict => 'f',
+  proretset => 't', provolatile => 's', prorettype => 'record',
+  proargtypes => 'pg_histogram int4',
+  proallargtypes => '{pg_histogram,int4,int4,_text,_text,_bool,_bool,_bool,float8,float8,float8}',
+  proargmodes => '{i,i,o,o,o,o,o,o,o,o,o}',
+  proargnames => '{histogram,otype,index,minvals,maxvals,nullsonly,mininclusive,maxinclusive,frequency,density,bucket_volume}',
+  prosrc => 'pg_histogram_buckets' },
+
 { oid => '1928', descr => 'statistics: number of scans done for table/index',
   proname => 'pg_stat_get_numscans', provolatile => 's', proparallel => 'r',
   prorettype => 'int8', proargtypes => 'oid',
diff --git a/src/include/catalog/pg_statistic_ext.h b/src/include/catalog/pg_statistic_ext.h
index 7ddbee63c9..514f4230c9 100644
--- a/src/include/catalog/pg_statistic_ext.h
+++ b/src/include/catalog/pg_statistic_ext.h
@@ -48,6 +48,7 @@ CATALOG(pg_statistic_ext,3381,StatisticExtRelationId)
 	pg_ndistinct stxndistinct;	/* ndistinct coefficients (serialized) */
 	pg_dependencies stxdependencies;	/* dependencies (serialized) */
 	pg_mcv_list stxmcv;			/* MCV (serialized) */
+	pg_histogram stxhistogram;	/* MV histogram (serialized) */
 #endif
 
 } FormData_pg_statistic_ext;
@@ -64,6 +65,7 @@ typedef FormData_pg_statistic_ext *Form_pg_statistic_ext;
 #define STATS_EXT_NDISTINCT			'd'
 #define STATS_EXT_DEPENDENCIES		'f'
 #define STATS_EXT_MCV				'm'
+#define STATS_EXT_HISTOGRAM			'h'
 
 #endif							/* EXPOSE_TO_CLIENT_CODE */
 
diff --git a/src/include/catalog/pg_type.dat b/src/include/catalog/pg_type.dat
index 0ff25e87a7..15e3ab9f93 100644
--- a/src/include/catalog/pg_type.dat
+++ b/src/include/catalog/pg_type.dat
@@ -172,6 +172,13 @@
   typoutput => 'pg_mcv_list_out', typreceive => 'pg_mcv_list_recv',
   typsend => 'pg_mcv_list_send', typalign => 'i', typstorage => 'x',
   typcollation => '100' },
+{ oid => '3425', oid_symbol => 'PGHISTOGRAMOID',
+  descr => 'multivariate histogram',
+  typname => 'pg_histogram', typlen => '-1', typbyval => 'f',
+  typcategory => 'S', typinput => 'pg_histogram_in',
+  typoutput => 'pg_histogram_out', typreceive => 'pg_histogram_recv',
+  typsend => 'pg_histogram_send', typalign => 'i', typstorage => 'x',
+  typcollation => '100' },
 { oid => '32', oid_symbol => 'PGDDLCOMMANDOID',
   descr => 'internal type for passing CollectedCommand',
   typname => 'pg_ddl_command', typlen => 'SIZEOF_POINTER', typbyval => 't',
diff --git a/src/include/nodes/relation.h b/src/include/nodes/relation.h
index adb4265047..5849a1d0d4 100644
--- a/src/include/nodes/relation.h
+++ b/src/include/nodes/relation.h
@@ -859,10 +859,15 @@ typedef struct StatisticExtInfo
 
 	Oid			statOid;		/* OID of the statistics row */
 	RelOptInfo *rel;			/* back-link to statistic's table */
-	char		kind;			/* statistic kind of this entry */
+	int			kinds;			/* statistic kinds of this entry */
 	Bitmapset  *keys;			/* attnums of the columns covered */
 } StatisticExtInfo;
 
+#define STATS_EXT_INFO_NDISTINCT			1
+#define STATS_EXT_INFO_DEPENDENCIES			2
+#define STATS_EXT_INFO_MCV					4
+#define STATS_EXT_INFO_HISTOGRAM			8
+
 /*
  * EquivalenceClasses
  *
diff --git a/src/include/statistics/extended_stats_internal.h b/src/include/statistics/extended_stats_internal.h
index 11159b58ee..ddb0066378 100644
--- a/src/include/statistics/extended_stats_internal.h
+++ b/src/include/statistics/extended_stats_internal.h
@@ -69,10 +69,16 @@ extern MVDependencies *statext_dependencies_deserialize(bytea *data);
 
 extern MCVList * statext_mcv_build(int numrows, HeapTuple *rows,
 								   Bitmapset *attrs, VacAttrStats **stats,
+								   HeapTuple **rows_filtered, int *numrows_filtered,
 								   double totalrows);
 extern bytea *statext_mcv_serialize(MCVList * mcv, VacAttrStats **stats);
 extern MCVList * statext_mcv_deserialize(bytea *data);
 
+extern MVHistogram * statext_histogram_build(int numrows, HeapTuple *rows,
+											 Bitmapset *attrs, VacAttrStats **stats,
+											 int numrows_total);
+extern MVHistogram * statext_histogram_deserialize(bytea *data);
+
 extern MultiSortSupport multi_sort_init(int ndims);
 extern void multi_sort_add_dimension(MultiSortSupport mss, int sortdim,
 						 Oid oper);
@@ -83,6 +89,7 @@ extern int multi_sort_compare_dims(int start, int end, const SortItem *a,
 						const SortItem *b, MultiSortSupport mss);
 extern int	compare_scalars_simple(const void *a, const void *b, void *arg);
 extern int	compare_datums_simple(Datum a, Datum b, SortSupport ssup);
+extern int	compare_scalars_partition(const void *a, const void *b, void *arg);
 
 extern void *bsearch_arg(const void *key, const void *base,
 			size_t nmemb, size_t size,
@@ -109,4 +116,12 @@ extern Selectivity mcv_clauselist_selectivity(PlannerInfo *root,
 						   Selectivity *basesel,
 						   Selectivity *totalsel);
 
+extern Selectivity histogram_clauselist_selectivity(PlannerInfo *root,
+								 StatisticExtInfo *stat,
+								 List *clauses, List *conditions,
+								 int varRelid,
+								 JoinType jointype,
+								 SpecialJoinInfo *sjinfo,
+								 RelOptInfo *rel);
+
 #endif							/* EXTENDED_STATS_INTERNAL_H */
diff --git a/src/include/statistics/statistics.h b/src/include/statistics/statistics.h
index e69d6a0232..1d276f0b6d 100644
--- a/src/include/statistics/statistics.h
+++ b/src/include/statistics/statistics.h
@@ -119,9 +119,68 @@ typedef struct MCVList
 	MCVItem   **items;			/* array of MCV items */
 }			MCVList;
 
+
+/* used to flag stats serialized to bytea */
+#define STATS_HIST_MAGIC       0x7F8C5670	/* marks serialized bytea */
+#define STATS_HIST_TYPE_BASIC  1	/* basic histogram type */
+
+/* max buckets in a histogram (mostly arbitrary number) */
+#define STATS_HIST_MAX_BUCKETS 16384
+
+/*
+ * Histogram in a partially serialized form, with deduplicated boundary
+ * values etc.
+ */
+typedef struct MVBucket
+{
+	/* Frequencies of this bucket. */
+	float		frequency;
+
+	/*
+	 * Information about dimensions being NULL-only. Not yet used.
+	 */
+	bool	   *nullsonly;
+
+	/* lower boundaries - values and information about the inequalities */
+	uint16	   *min;
+	bool	   *min_inclusive;
+
+	/*
+	 * indexes of upper boundaries - values and information about the
+	 * inequalities (exclusive vs. inclusive)
+	 */
+	uint16	   *max;
+	bool	   *max_inclusive;
+}			MVBucket;
+
+typedef struct MVHistogram
+{
+	/* varlena header (do not touch directly!) */
+	int32		vl_len_;
+	uint32		magic;			/* magic constant marker */
+	uint32		type;			/* type of histogram (BASIC) */
+	uint32		nbuckets;		/* number of buckets (buckets array) */
+	uint32		ndimensions;	/* number of dimensions */
+	Oid			types[STATS_MAX_DIMENSIONS];	/* OIDs of data types */
+
+	/*
+	 * keep this the same with MVHistogram, because of deserialization (same
+	 * offset)
+	 */
+	MVBucket  **buckets;		/* array of buckets */
+
+	/*
+	 * serialized boundary values, one array per dimension, deduplicated (the
+	 * min/max indexes point into these arrays)
+	 */
+	int		   *nvalues;
+	Datum	  **values;
+}			MVHistogram;
+
 extern MVNDistinct *statext_ndistinct_load(Oid mvoid);
 extern MVDependencies *statext_dependencies_load(Oid mvoid);
 extern MCVList * statext_mcv_load(Oid mvoid);
+extern MVHistogram * statext_histogram_load(Oid mvoid);
 
 extern void BuildRelationExtStatistics(Relation onerel, double totalrows,
 						   int numrows, HeapTuple *rows,
@@ -141,8 +200,8 @@ extern Selectivity statext_clauselist_selectivity(PlannerInfo *root,
 							   SpecialJoinInfo *sjinfo,
 							   RelOptInfo *rel,
 							   Bitmapset **estimatedclauses);
-extern bool has_stats_of_kind(List *stats, char requiredkind);
+extern bool has_stats_of_kind(List *stats, int requiredkinds);
 extern StatisticExtInfo *choose_best_statistics(List *stats,
-					   Bitmapset *attnums, char requiredkind);
+					   Bitmapset *attnums, int requiredkinds);
 
 #endif							/* STATISTICS_H */
diff --git a/src/include/utils/selfuncs.h b/src/include/utils/selfuncs.h
index 4e9aaca6b5..6c9dd1d85a 100644
--- a/src/include/utils/selfuncs.h
+++ b/src/include/utils/selfuncs.h
@@ -222,6 +222,10 @@ extern void genericcostestimate(PlannerInfo *root, IndexPath *path,
 					List *qinfos,
 					GenericCosts *costs);
 
+extern bool convert_to_scalar(Datum value, Oid valuetypid, double *scaledvalue,
+				  Datum lobound, Datum hibound, Oid boundstypid,
+				  double *scaledlobound, double *scaledhibound);
+
 /* Functions in array_selfuncs.c */
 
 extern Selectivity scalararraysel_containment(PlannerInfo *root,
diff --git a/src/test/regress/expected/create_table_like.out b/src/test/regress/expected/create_table_like.out
index 0f97355165..fc672f6cd6 100644
--- a/src/test/regress/expected/create_table_like.out
+++ b/src/test/regress/expected/create_table_like.out
@@ -243,7 +243,7 @@ Indexes:
 Check constraints:
     "ctlt1_a_check" CHECK (length(a) > 2)
 Statistics objects:
-    "public"."ctlt_all_a_b_stat" (ndistinct, dependencies, mcv) ON a, b FROM ctlt_all
+    "public"."ctlt_all_a_b_stat" (ndistinct, dependencies, mcv, histogram) ON a, b FROM ctlt_all
 
 SELECT c.relname, objsubid, description FROM pg_description, pg_index i, pg_class c WHERE classoid = 'pg_class'::regclass AND objoid = i.indexrelid AND c.oid = i.indexrelid AND i.indrelid = 'ctlt_all'::regclass ORDER BY c.relname, objsubid;
     relname     | objsubid | description 
diff --git a/src/test/regress/expected/opr_sanity.out b/src/test/regress/expected/opr_sanity.out
index d13f0928d0..174efe3ce0 100644
--- a/src/test/regress/expected/opr_sanity.out
+++ b/src/test/regress/expected/opr_sanity.out
@@ -903,11 +903,12 @@ WHERE c.castmethod = 'b' AND
  pg_ndistinct      | bytea             |        0 | i
  pg_dependencies   | bytea             |        0 | i
  pg_mcv_list       | bytea             |        0 | i
+ pg_histogram      | bytea             |        0 | i
  cidr              | inet              |        0 | i
  xml               | text              |        0 | a
  xml               | character varying |        0 | a
  xml               | character         |        0 | a
-(10 rows)
+(11 rows)
 
 -- **************** pg_conversion ****************
 -- Look for illegal values in pg_conversion fields.
diff --git a/src/test/regress/expected/stats_ext.out b/src/test/regress/expected/stats_ext.out
index 5d05962c04..67975c91d3 100644
--- a/src/test/regress/expected/stats_ext.out
+++ b/src/test/regress/expected/stats_ext.out
@@ -58,7 +58,7 @@ ALTER TABLE ab1 DROP COLUMN a;
  b      | integer |           |          | 
  c      | integer |           |          | 
 Statistics objects:
-    "public"."ab1_b_c_stats" (ndistinct, dependencies, mcv) ON b, c FROM ab1
+    "public"."ab1_b_c_stats" (ndistinct, dependencies, mcv, histogram) ON b, c FROM ab1
 
 -- Ensure statistics are dropped when table is
 SELECT stxname FROM pg_statistic_ext WHERE stxname LIKE 'ab1%';
@@ -204,9 +204,9 @@ CREATE STATISTICS s10 ON a, b, c FROM ndistinct;
 ANALYZE ndistinct;
 SELECT stxkind, stxndistinct
   FROM pg_statistic_ext WHERE stxrelid = 'ndistinct'::regclass;
- stxkind |                      stxndistinct                       
----------+---------------------------------------------------------
- {d,f,m} | {"3, 4": 301, "3, 6": 301, "4, 6": 301, "3, 4, 6": 301}
+  stxkind  |                      stxndistinct                       
+-----------+---------------------------------------------------------
+ {d,f,m,h} | {"3, 4": 301, "3, 6": 301, "4, 6": 301, "3, 4, 6": 301}
 (1 row)
 
 -- Hash Aggregate, thanks to estimates improved by the statistic
@@ -270,9 +270,9 @@ INSERT INTO ndistinct (a, b, c, filler1)
 ANALYZE ndistinct;
 SELECT stxkind, stxndistinct
   FROM pg_statistic_ext WHERE stxrelid = 'ndistinct'::regclass;
- stxkind |                        stxndistinct                         
----------+-------------------------------------------------------------
- {d,f,m} | {"3, 4": 2550, "3, 6": 800, "4, 6": 1632, "3, 4, 6": 10000}
+  stxkind  |                        stxndistinct                         
+-----------+-------------------------------------------------------------
+ {d,f,m,h} | {"3, 4": 2550, "3, 6": 800, "4, 6": 1632, "3, 4, 6": 10000}
 (1 row)
 
 -- plans using Group Aggregate, thanks to using correct esimates
@@ -758,7 +758,6 @@ EXPLAIN (COSTS OFF)
          Index Cond: ((a IS NULL) AND (b IS NULL))
 (5 rows)
 
-RESET random_page_cost;
 -- mcv with arrays
 CREATE TABLE mcv_lists_arrays (
     a TEXT[],
@@ -822,3 +821,197 @@ EXPLAIN (COSTS OFF) SELECT * FROM mcv_lists_bool WHERE NOT a AND b AND NOT c;
    Filter: ((NOT a) AND b AND (NOT c))
 (3 rows)
 
+RESET random_page_cost;
+-- histograms
+CREATE TABLE histograms (
+    filler1 TEXT,
+    filler2 NUMERIC,
+    a INT,
+    b TEXT,
+    filler3 DATE,
+    c INT,
+    d TEXT
+);
+SET random_page_cost = 1.2;
+CREATE INDEX histograms_ab_idx ON mcv_lists (a, b);
+CREATE INDEX histograms_abc_idx ON histograms (a, b, c);
+-- random data (we still get histogram, but as the columns are not
+-- correlated, the estimates remain about the same)
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT mod(i,37), mod(i,41), mod(i,43), mod(i,47) FROM generate_series(1,5000) s(i);
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5';
+                    QUERY PLAN                     
+---------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 5) AND (b < '5'::text))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 5) AND (b < '5'::text))
+(4 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5' AND c < 5;
+                          QUERY PLAN                           
+---------------------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 5) AND (b < '5'::text) AND (c < 5))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 5) AND (b < '5'::text) AND (c < 5))
+(4 rows)
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5';
+                    QUERY PLAN                     
+---------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 5) AND (b < '5'::text))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 5) AND (b < '5'::text))
+(4 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5' AND c < 5;
+                          QUERY PLAN                           
+---------------------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 5) AND (b < '5'::text) AND (c < 5))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 5) AND (b < '5'::text) AND (c < 5))
+(4 rows)
+
+-- values correlated along the diagonal
+TRUNCATE histograms;
+DROP STATISTICS histograms_stats;
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT mod(i,100), mod(i,100) + mod(i,7), mod(i,100) + mod(i,11), i FROM generate_series(1,5000) s(i);
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND c < 3;
+                    QUERY PLAN                     
+---------------------------------------------------
+ Index Scan using histograms_abc_idx on histograms
+   Index Cond: ((a < 3) AND (c < 3))
+(2 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND b > '2' AND c < 3;
+                       QUERY PLAN                        
+---------------------------------------------------------
+ Index Scan using histograms_abc_idx on histograms
+   Index Cond: ((a < 3) AND (b > '2'::text) AND (c < 3))
+(2 rows)
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND c < 3;
+                  QUERY PLAN                   
+-----------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 3) AND (c < 3))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 3) AND (c < 3))
+(4 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND b > '2' AND c < 3;
+                          QUERY PLAN                           
+---------------------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a < 3) AND (b > '2'::text) AND (c < 3))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a < 3) AND (b > '2'::text) AND (c < 3))
+(4 rows)
+
+-- almost 5000 unique combinations with NULL values
+TRUNCATE histograms;
+DROP STATISTICS histograms_stats;
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT
+         (CASE WHEN mod(i,100) =  0 THEN NULL ELSE mod(i,100) END),
+         (CASE WHEN mod(i,100) <= 1 THEN NULL ELSE mod(i,100) + mod(i,7)  END),
+         (CASE WHEN mod(i,100) <= 2 THEN NULL ELSE mod(i,100) + mod(i,11) END),
+         i
+     FROM generate_series(1,5000) s(i);
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+                    QUERY PLAN                     
+---------------------------------------------------
+ Index Scan using histograms_abc_idx on histograms
+   Index Cond: ((a IS NULL) AND (b IS NULL))
+(2 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL AND c IS NULL;
+                         QUERY PLAN                          
+-------------------------------------------------------------
+ Index Scan using histograms_abc_idx on histograms
+   Index Cond: ((a IS NULL) AND (b IS NULL) AND (c IS NULL))
+(2 rows)
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+                    QUERY PLAN                     
+---------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a IS NULL) AND (b IS NULL))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a IS NULL) AND (b IS NULL))
+(4 rows)
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL AND c IS NULL;
+                            QUERY PLAN                             
+-------------------------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a IS NULL) AND (b IS NULL) AND (c IS NULL))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a IS NULL) AND (b IS NULL) AND (c IS NULL))
+(4 rows)
+
+-- check change of column type resets the histogram statistics
+ALTER TABLE histograms ALTER COLUMN c TYPE numeric;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+                    QUERY PLAN                     
+---------------------------------------------------
+ Index Scan using histograms_abc_idx on histograms
+   Index Cond: ((a IS NULL) AND (b IS NULL))
+(2 rows)
+
+ANALYZE histograms;
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+                    QUERY PLAN                     
+---------------------------------------------------
+ Bitmap Heap Scan on histograms
+   Recheck Cond: ((a IS NULL) AND (b IS NULL))
+   ->  Bitmap Index Scan on histograms_abc_idx
+         Index Cond: ((a IS NULL) AND (b IS NULL))
+(4 rows)
+
+-- histograms with arrays
+CREATE TABLE histograms_arrays (
+    a TEXT[],
+    b NUMERIC[],
+    c INT[]
+);
+INSERT INTO histograms_arrays (a, b, c)
+     SELECT
+         ARRAY[md5(i::text), md5((i-1)::text), md5((i+1)::text)],
+         ARRAY[(i-1)::numeric/1000, i::numeric/1000, (i+1)::numeric/1000],
+         ARRAY[(i-1), i, (i+1)]
+     FROM generate_series(1,5000) s(i);
+CREATE STATISTICS histogram_array_stats (histogram) ON a, b, c
+  FROM histograms_arrays;
+ANALYZE histograms_arrays;
+RESET random_page_cost;
diff --git a/src/test/regress/expected/type_sanity.out b/src/test/regress/expected/type_sanity.out
index a56d6c5231..97c292f6f9 100644
--- a/src/test/regress/expected/type_sanity.out
+++ b/src/test/regress/expected/type_sanity.out
@@ -73,8 +73,9 @@ WHERE p1.typtype not in ('c','d','p') AND p1.typname NOT LIKE E'\\_%'
  3361 | pg_ndistinct
  3402 | pg_dependencies
  4001 | pg_mcv_list
+ 3425 | pg_histogram
   210 | smgr
-(5 rows)
+(6 rows)
 
 -- Make sure typarray points to a varlena array type of our own base
 SELECT p1.oid, p1.typname as basetype, p2.typname as arraytype,
diff --git a/src/test/regress/sql/stats_ext.sql b/src/test/regress/sql/stats_ext.sql
index ad1f103217..a949c7e6d1 100644
--- a/src/test/regress/sql/stats_ext.sql
+++ b/src/test/regress/sql/stats_ext.sql
@@ -414,8 +414,6 @@ EXPLAIN (COSTS OFF)
 EXPLAIN (COSTS OFF)
  SELECT * FROM mcv_lists WHERE a IS NULL AND b IS NULL AND c IS NULL;
 
-RESET random_page_cost;
-
 -- mcv with arrays
 CREATE TABLE mcv_lists_arrays (
     a TEXT[],
@@ -463,3 +461,134 @@ EXPLAIN (COSTS OFF) SELECT * FROM mcv_lists_bool WHERE NOT a AND b AND c;
 EXPLAIN (COSTS OFF) SELECT * FROM mcv_lists_bool WHERE NOT a AND NOT b AND c;
 
 EXPLAIN (COSTS OFF) SELECT * FROM mcv_lists_bool WHERE NOT a AND b AND NOT c;
+
+RESET random_page_cost;
+
+-- histograms
+CREATE TABLE histograms (
+    filler1 TEXT,
+    filler2 NUMERIC,
+    a INT,
+    b TEXT,
+    filler3 DATE,
+    c INT,
+    d TEXT
+);
+
+SET random_page_cost = 1.2;
+
+CREATE INDEX histograms_ab_idx ON mcv_lists (a, b);
+CREATE INDEX histograms_abc_idx ON histograms (a, b, c);
+
+-- random data (we still get histogram, but as the columns are not
+-- correlated, the estimates remain about the same)
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT mod(i,37), mod(i,41), mod(i,43), mod(i,47) FROM generate_series(1,5000) s(i);
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5';
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5' AND c < 5;
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5';
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 5 AND b < '5' AND c < 5;
+
+-- values correlated along the diagonal
+TRUNCATE histograms;
+DROP STATISTICS histograms_stats;
+
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT mod(i,100), mod(i,100) + mod(i,7), mod(i,100) + mod(i,11), i FROM generate_series(1,5000) s(i);
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND c < 3;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND b > '2' AND c < 3;
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND c < 3;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a < 3 AND b > '2' AND c < 3;
+
+-- almost 5000 unique combinations with NULL values
+TRUNCATE histograms;
+DROP STATISTICS histograms_stats;
+
+INSERT INTO histograms (a, b, c, filler1)
+     SELECT
+         (CASE WHEN mod(i,100) =  0 THEN NULL ELSE mod(i,100) END),
+         (CASE WHEN mod(i,100) <= 1 THEN NULL ELSE mod(i,100) + mod(i,7)  END),
+         (CASE WHEN mod(i,100) <= 2 THEN NULL ELSE mod(i,100) + mod(i,11) END),
+         i
+     FROM generate_series(1,5000) s(i);
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL AND c IS NULL;
+
+-- create statistics
+CREATE STATISTICS histograms_stats (histogram) ON a, b, c FROM histograms;
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL AND c IS NULL;
+
+-- check change of column type resets the histogram statistics
+ALTER TABLE histograms ALTER COLUMN c TYPE numeric;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+
+ANALYZE histograms;
+
+EXPLAIN (COSTS OFF)
+ SELECT * FROM histograms WHERE a IS NULL AND b IS NULL;
+
+-- histograms with arrays
+CREATE TABLE histograms_arrays (
+    a TEXT[],
+    b NUMERIC[],
+    c INT[]
+);
+
+INSERT INTO histograms_arrays (a, b, c)
+     SELECT
+         ARRAY[md5(i::text), md5((i-1)::text), md5((i+1)::text)],
+         ARRAY[(i-1)::numeric/1000, i::numeric/1000, (i+1)::numeric/1000],
+         ARRAY[(i-1), i, (i+1)]
+     FROM generate_series(1,5000) s(i);
+
+CREATE STATISTICS histogram_array_stats (histogram) ON a, b, c
+  FROM histograms_arrays;
+
+ANALYZE histograms_arrays;
+
+RESET random_page_cost;
-- 
2.13.6

