From 2194a785da379cd7dfe64eadccdc817eedd368eb Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 2 Apr 2024 17:54:43 +0200
Subject: [PATCH] efficient simplicial numbers in miscfunc

---
 miscfunc.cc | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++
 miscfunc.h  |  1 +
 t.cc        | 10 +++++++++-
 3 files changed, 61 insertions(+), 1 deletion(-)

diff --git a/miscfunc.cc b/miscfunc.cc
index 9a735c8..aa095f2 100644
--- a/miscfunc.cc
+++ b/miscfunc.cc
@@ -758,5 +758,56 @@ return(bern[n]);
 
 
 
+//there is a closed form s(d,n) = binom(n+d-1,d)
+//but that would be inefficient
+//we use recurrence s(d,n) = s(d,n-1) + s(d-1,n)
+//we cannot also use the above efficient binom() since here we need laerge n while d is small
+//
+//enlarge storage size if needed
+//
+#define SIMPLICIAL_MAXD 8
+#define SIMPLICIAL_MAXN 1024
+
+unsigned long long simplicial(int d, int n)
+{
+static unsigned long long stored[SIMPLICIAL_MAXD-2][SIMPLICIAL_MAXN]={{0,1,3}};
+static int maxd=2;
+static int maxn=2;
+
+if(n<=0) return 0;
+if(d==0||n==1) return 1;
+if(d==1) return n;
+
+//resize precomputed part as needed
+if(n>maxn)
+	{
+	if(n>=SIMPLICIAL_MAXN) laerror("storage overflow in simplicial");
+	int newdim=2*n; 
+	if(newdim<2*maxn) newdim=2*maxn;
+	if(newdim>=SIMPLICIAL_MAXN) newdim=SIMPLICIAL_MAXN-1;
+	for(int i=maxn+1; i<=newdim;++i) stored[0][i] = i*(i+1)/2;
+	for(int dd=1; dd<=maxd-2; ++dd)
+		{
+		for(int i=maxn+1; i<=newdim;++i) stored[dd][i] =  stored[dd][i-1] + stored[dd-1][i];
+		}
+	maxn=newdim;
+	}
+
+if(d>maxd)
+	{
+	if(d>=SIMPLICIAL_MAXD) laerror("storage overflow in simplicial");
+	for(int dd=maxd+1; dd<=d; ++dd)
+		{
+		stored[dd-2][0]=0;
+		stored[dd-2][1]=1;
+		for(int i=2; i<=maxn; ++i) stored[dd-2][i] = stored[dd-3][i] + stored[dd-2][i-1];
+		}
+	maxd=d;
+	}
+
+return stored[d-2][n];
+}
+
+
 
 }//namespace
diff --git a/miscfunc.h b/miscfunc.h
index cc1f50d..96009d7 100644
--- a/miscfunc.h
+++ b/miscfunc.h
@@ -35,6 +35,7 @@ extern double zeta(double x);
 extern double gamma(double x);
 extern double gammln(double x);
 extern double bernoulli_number(int n);
+extern unsigned long long simplicial(int d, int n); //simplicial numbers in a precomputed efficient way
 
 
 
diff --git a/t.cc b/t.cc
index 23a2d6b..f8a2a58 100644
--- a/t.cc
+++ b/t.cc
@@ -28,6 +28,7 @@
 #include "simple.h"
 #include "graph.h"
 #include "numbers.h"
+#include "miscfunc.h"
 
 
 using namespace std;
@@ -3161,7 +3162,7 @@ cout <<a.rsum()<<endl;
 cout <<a.csum()<<endl;
 }
 
-if(1)
+if(0)
 {
 int n;
 cin>>n;
@@ -3175,5 +3176,12 @@ for(int l=1; l<k; ++l)
 cout <<count<<" "<<binom(n,4)<<endl;
 }
 
+if(1)
+{
+int d,n;
+cin>>d>>n;
+cout <<simplicial(d,n)<<" "<<binom(n+d-1,d)<<endl;
+}
+
 
 }