Add files via upload

main.cpp

@@ -11,29 +11,88 @@
 
 using namespace std; 
 
-const long MAXN = 1e4;
+//const size_t MAXN = (size_t)1e3;                 //size of YT
+const double EPS = 1e-8;
 long long factorial[20];
 //const unsigned int N = 5;
-long double branch_prob[MAXN];
+//long double branch_prob[MAXN];
 
 
 int main() {
-	//freopen("out.txt", "w", stdout);
-	freopen("err.txt", "w", stderr);
+	freopen("out.txt", "w", stdout);
+	//freopen("err.txt", "w", stderr);
+	cout.precision(6);
 
-	long b[10];
-	for (int i = 0; i < 10; ++i)
+	/*
+	double add = 1 / (double)MAXN;
+	double r = 0;
+	for (size_t i = 0; i < MAXN; ++i)
+	{
+		ra[i] = r;
+		r += add;
+	}
+
+	rand_permute_self_inverse(ra, MAXN);
+	*/
+
+	const long MAN = (long)2e2;
+	long *func = new long[MAN];
+	long *height = new long[MAN];
+	long *length = new long[MAN];
+	for (size_t MAXN = 1; MAXN < MAN; ++MAXN)
 	{
-		b[i] = i;
-		cout << i << ' ';
+		cerr << MAXN << "\n";
+		double *ra = 0;
+		ra = new double[MAXN];
+		rand_arr(ra, MAXN);
+		PQ_tableaux<double> PQ;
+		PQ.push(ra, MAXN);
+		height[MAXN] = PQ.height();
+		length[MAXN] = PQ.length();
+
+		double *max_ra = new double[MAXN];
+		for (size_t cnt = 0; cnt * cnt < 16 * MAXN; ++cnt)
+		{
+			for (size_t i = 0; i < MAXN; ++i)
+			{
+				double tmp = 0;
+				for (size_t j = 0; j < i; ++j)
+				{
+					tmp = max(tmp, (ra[j] < ra[i]) ? ra[j] : 0);
+				}
+				max_ra[i] = tmp;
+			}
+			PQ_tableaux<double> PQ_max;
+			PQ_max.push(max_ra, MAXN);
+			if (PQ_max.height() == 1)
+			{
+				func[MAXN] = cnt;
+				break;
+			}
+			for (size_t i = 0; i < MAXN; ++i)
+			{
+				ra[i] = max_ra[i];
+			}
+		}
 	}
-	cout << "\n";
-	rand_permute_self_inverse<long>(b, 10);
-	for (int i = 0; i < 10; ++i)
+
+	for (size_t i = 0; i < MAN; ++i)
 	{
-		cout << b[i] << ' ';
+		cout << i << ' ' << func[i] << ' ' << length[i] << ' ' << height[i] << "\n";
 	}
-	cout << "\n";
+
+	/*
+	PQ.showP();
+	cout << "\n\n";
+	PQ.truncate(max_ra[9], ra[9]);
+	PQ.showP();
+	cout << "\n\n";
+	PQ.truncate(max_ra[89], ra[89]);
+	PQ.showP();
+	cout << "\n\n";
+	*/
+
+
 
 	/*
 	int n = 6000;
@@ -97,7 +156,7 @@ int main() {
 		cout << "id_points statistics for all permutation: " << measure << "\n" << "\n";
 	}
 	*/
-	system("pause");
+	//system("pause");
 	return 0;
 	/*
 	PQ_tableaux PQ;

my_vector.h

@@ -0,0 +1,49 @@
+#pragma once
+#include <iostream>
+#include <vector>
+
+using namespace std;
+
+template<typename T>
+ostream& operator <<(ostream&, const vector<T>&);        //Output operator for vector
+template<typename T>
+vector<double> operator +(
+	const vector<T> &,
+	const vector<T> &
+	);
+template<typename T>
+T sum(const vector<T> &);
+
+
+
+//implementation
+
+template<typename T>
+ostream& operator <<(ostream& os, const vector<T>& out)
+{
+	for (auto i : out)
+	{
+		os << i << ' ';
+	}
+	return os;
+}
+template<typename T>
+vector<double> operator+(const vector<T> & vec1, const vector<T> & vec2)
+{
+	vector<double> vec_res(max(vec1.size(), vec2.size()));
+	for (int i = 0; i < (int)min(vec1.size(), vec2.size()); ++i)
+	{
+		vec_res[i] = vec1[i] + vec2[i];
+	}
+	return vec_res;
+}
+template<typename T>
+T sum(const vector<T>& vec)
+{
+	T sum = 0;
+	for (auto i : vec)
+	{
+		sum += i;
+	}
+	return sum;
+}

permutations.cpp

@@ -0,0 +1,15 @@
+#include <vector>
+#include "permutations.h"
+#include <algorithm>
+
+using namespace std;
+
+void init_involution_branch_prob(long double * b, long n)  //count branch prob. for involutions
+{
+	b[1] = 1;
+	for (int i = 1; i < n; ++i)
+	{
+		b[i + 1] = 1 / (1 + i * b[i]);
+	}
+	return;
+}

permutations.h

@@ -0,0 +1,93 @@
+#pragma once
+
+#include <vector>
+#include "random.h"
+
+using namespace std;
+
+template<typename T>
+void permutate(vector<T> &, const vector<int> &);        //Apply permutation to vector
+template<typename T>
+int count_identity_points(vector<T> &);                  //Count the number of identity points
+template<typename T>
+inline void rand_permute_self_inverse(
+	T*, long, long *,
+	long double *, bool);                                //Return involution of length n
+void init_involution_branch_prob(long double*, long);    //Count branch probabilities for generating random involution 
+
+
+
+//implementation
+
+
+template<typename T>
+void permutate(vector<T>& value, const vector<int>& permutation)
+{
+	auto new_value = value;
+	for (int i = 0; i < (int)value.size(); ++i)
+	{
+		new_value[permutation[i]] = value[i];
+	}
+	for (int i = 0; i < (int)value.size(); ++i)
+	{
+		value[i] = new_value[i];
+	}
+	return;
+}
+template<typename T>
+int count_identity_points(vector<T>& vec)
+{
+	vector<pair<T, int>> new_vec(vec.size());
+	for (int i = 0; i < (int)vec.size(); ++i)
+	{
+		new_vec[i] = { vec[i], i };
+	}
+	sort(new_vec.begin(), new_vec.end());
+	int identity_points = 0;
+	for (int i = 0; i < (int)vec.size(); ++i)
+	{
+		identity_points += new_vec[i].second == i;
+	}
+	return identity_points;
+}
+template<typename T>                                    //not finished yet, destructor needs some treatment
+void rand_permute_self_inverse(T *f, long n,
+	long *tr = 0,
+	long double *tb = 0, bool bi = false)
+	// Permute the elements of f by a random involution.
+	// Set bi:=true to signal that the branch probabilities in tb[]
+	// have been precomputed (via init_involution_branch_ratios()).
+{
+	long *r = tr;
+	if (tr == 0) r = new long[n];
+	for (long k = 0; k < n; ++k) r[k] = k;
+	long nr = n;                                         // number of elements available
+														 // available positions are r[0], ..., r[nr-1]
+
+	long double *b = tb;
+	if (tb == 0) { b = new long double[n]; bi = false; }
+	if (!bi) init_involution_branch_prob(b, n);
+
+	while (nr >= 2)
+	{
+		const long x1 = nr - 1;                           // choose last element
+		const long r1 = r[x1];                            // available position
+														  // remove from set:
+		--nr;                                             // no swap needed if x1==last
+
+		const long double rat = b[nr];                    // probability to choose fixed point
+
+		const long double t = rand01();                   // 0 <= t < 1
+		if (t > rat)                                      // 2-cycle
+		{
+			const long x2 = randN(nr);
+			const long r2 = r[x2];                        // random available position != r1
+			--nr; swap(r[x2], r[nr]);                     // remove from set
+			swap(f[r1], f[r2]);                           // create a 2-cycle
+		}
+		// else fixed point, nothing to do
+	}
+
+	//if (tr == 0) delete[] r;
+	//if (tb == 0) delete[] b;
+}

random.cpp

@@ -0,0 +1,30 @@
+
+#include <random>
+#include <ctime>
+#include "random.h"
+
+using namespace std;
+
+long randN(long N)
+{
+	long tmp = rand() ^ (rand() << 16);
+	return tmp % N;
+}
+double rand01(double acc)
+{
+	static bool f = false;
+	if (!f)
+	{
+		double tme = (double)time(NULL);
+		srand(tme);
+		f = true;
+	}
+	double rand_num = ((rand() ^ (rand() << 16)) % (int)acc) / acc;
+	return rand_num;
+}
+void rand_arr(double *ra, size_t n)
+{
+	for (size_t i = 0; i < n; i++)
+		ra[i] = rand01();
+	return;
+}

random.h

@@ -0,0 +1,7 @@
+#pragma once
+
+
+
+double rand01(double acc = 1e6);                             //Return real random number from [0, 1]
+long randN(long);                            //Return integer random number from [0, N-1], where N is an argument     
+void rand_arr(double *, size_t);             //Return an array of real random numbers form [0, 1]