Implemented the zeta and gamma functions

Author: Sekqies

src/shaders/embedded_shaders.h

@@ -53,6 +53,7 @@ const float TWO_OVER_PI = 2.0f / PI;
 
 const vec2 CPI = vec2(PI,0.0f);
 const vec2 ONE = vec2(1.0f,0.0f);
+const vec2 MINUS_ONE = vec2(-1.0f,0.0f);
 const vec2 I = vec2(0.0f,1.0f);
 
 #define END_CONSTANT_DEFINITIONS HERE
@@ -260,6 +261,82 @@ vec2 im(vec2 z){
     return vec2(z.y,0.0f);
 }
 
+// Non-elementary functions
+
+const int GAMMA_PRECISION = 7;
+
+const float GAMMA_COEFFICIENTS[7] = float[](
+    1.000000000190015, 76.18009172947146, -86.50532032941677,
+    24.01409824083091, -1.231739572450155, 1.208650973866179e-3, -5.395239384953e-6
+);
+
+const vec2 SQUARE_ROOT_TWO_PI = vec2(sqrt(PI * 2.0f),0.0f);
+
+vec2 cgamma(vec2 z){
+    vec2 sum = vec2(GAMMA_COEFFICIENTS[0],0.0);
+    for(int i = 1; i < 7; ++i) {
+        sum = cadd(sum, cdiv(vec2(GAMMA_COEFFICIENTS[i], 0.0), cadd(z, vec2(float(i), 0.0))));
+    }
+    vec2 t = cadd(z, vec2(5.5, 0.0));
+    vec2 exp_term = cmult(cpow(t, cadd(z, vec2(0.5, 0.0))), cexp(-t));
+    vec2 result_pos = cdiv(cmult(cmult(exp_term, SQUARE_ROOT_TWO_PI), sum),z);
+
+    vec2 z_neg = cneg(z);
+    vec2 sum_neg = vec2(GAMMA_COEFFICIENTS[0],0.0);
+    for(int i = 1; i < 7; ++i) {
+        sum_neg = cadd(sum_neg, cdiv(vec2(GAMMA_COEFFICIENTS[i], 0.0), cadd(z_neg, vec2(float(i), 0.0))));
+    }
+    vec2 t_neg = cadd(z_neg,vec2(5.5,0.0));
+    vec2 exp_term_neg = cmult(cpow(t_neg, cadd(z_neg, vec2(0.5, 0.0))), cexp(-t_neg));
+    vec2 gamma_1_minus_z = cmult(cmult(exp_term_neg, SQUARE_ROOT_TWO_PI), sum_neg);
+    vec2 pi_z = cmult(z,vec2(PI,0.0));
+    vec2 result_neg = cdiv(vec2(PI,0.0),cmult(csin(pi_z),gamma_1_minus_z));
+
+    float is_positive = step(0.5,z.x);
+    return mix(result_neg,result_pos,is_positive);
+}
+
+const int ZETA_PRECISION = 10;
+const float ZETA_COEFFICIENTS [10] = float[](1.0, 513.0, 23041.0, 263169.0, 1153025.0, 2306049.0, 2191361.0, 950273.0, 171537.0, 9217.0);
+
+vec2 czeta_main_branch(vec2 s){
+    vec2 sum = vec2(0.0, 0.0);
+    float sign = 1.0;
+    for(int n = 0; n < 10; ++n){
+        vec2 n_to_neg_s = cpow(vec2(float(n+1),0.0), cmult(s, MINUS_ONE));
+
+        n_to_neg_s = cmult(n_to_neg_s, vec2(sign * ZETA_COEFFICIENTS[ ZETA_PRECISION - 1 - n], 0.0));
+        sum = cadd(sum,n_to_neg_s);
+        sign = -sign;
+    }
+    vec2 eta = cmult(sum,vec2(1.0/ZETA_COEFFICIENTS[ZETA_PRECISION - 1], 0.0));
+    vec2 one_minus_s = csub(ONE,s);
+    vec2 two_pow = cpow(vec2(2.0,0.0),one_minus_s);
+    return cdiv(eta,csub(ONE,two_pow));
+}
+
+vec2 czeta_negative_branch(vec2 s){
+    vec2 one_minus_s = csub(ONE, s);
+    vec2 term1 = cpow(vec2(2.0, 0.0), s);
+    vec2 term2 = cpow(vec2(PI, 0.0), csub(s, ONE));
+    vec2 term3 = csin(cmult(s, vec2(PI / 2.0, 0.0)));
+    vec2 term4 = cgamma(one_minus_s);
+    vec2 term5 = czeta_main_branch(one_minus_s);
+    return cmult(term1, cmult(term2, cmult(term3, cmult(term4, term5))));
+
+}
+
+vec2 czeta(vec2 s){
+    if (s.x >= 0.0) {
+        return czeta_main_branch(s);
+    } 
+    else {
+        return czeta_negative_branch(s);
+    }
+}
+
+
+
 #define END_FUNCTION_DEFINITIONS HERE
 
 

src/shaders/plotter.frag

@@ -15,6 +15,7 @@ const float TWO_OVER_PI = 2.0f / PI;
 
 const vec2 CPI = vec2(PI,0.0f);
 const vec2 ONE = vec2(1.0f,0.0f);
+const vec2 MINUS_ONE = vec2(-1.0f,0.0f);
 const vec2 I = vec2(0.0f,1.0f);
 
 #define END_CONSTANT_DEFINITIONS HERE
@@ -222,6 +223,82 @@ vec2 im(vec2 z){
     return vec2(z.y,0.0f);
 }
 
+// Non-elementary functions
+
+const int GAMMA_PRECISION = 7;
+
+const float GAMMA_COEFFICIENTS[7] = float[](
+    1.000000000190015, 76.18009172947146, -86.50532032941677,
+    24.01409824083091, -1.231739572450155, 1.208650973866179e-3, -5.395239384953e-6
+);
+
+const vec2 SQUARE_ROOT_TWO_PI = vec2(sqrt(PI * 2.0f),0.0f);
+
+vec2 cgamma(vec2 z){
+    vec2 sum = vec2(GAMMA_COEFFICIENTS[0],0.0);
+    for(int i = 1; i < 7; ++i) {
+        sum = cadd(sum, cdiv(vec2(GAMMA_COEFFICIENTS[i], 0.0), cadd(z, vec2(float(i), 0.0))));
+    }
+    vec2 t = cadd(z, vec2(5.5, 0.0));
+    vec2 exp_term = cmult(cpow(t, cadd(z, vec2(0.5, 0.0))), cexp(-t));
+    vec2 result_pos = cdiv(cmult(cmult(exp_term, SQUARE_ROOT_TWO_PI), sum),z);
+
+    vec2 z_neg = cneg(z);
+    vec2 sum_neg = vec2(GAMMA_COEFFICIENTS[0],0.0);
+    for(int i = 1; i < 7; ++i) {
+        sum_neg = cadd(sum_neg, cdiv(vec2(GAMMA_COEFFICIENTS[i], 0.0), cadd(z_neg, vec2(float(i), 0.0))));
+    }
+    vec2 t_neg = cadd(z_neg,vec2(5.5,0.0));
+    vec2 exp_term_neg = cmult(cpow(t_neg, cadd(z_neg, vec2(0.5, 0.0))), cexp(-t_neg));
+    vec2 gamma_1_minus_z = cmult(cmult(exp_term_neg, SQUARE_ROOT_TWO_PI), sum_neg);
+    vec2 pi_z = cmult(z,vec2(PI,0.0));
+    vec2 result_neg = cdiv(vec2(PI,0.0),cmult(csin(pi_z),gamma_1_minus_z));
+
+    float is_positive = step(0.5,z.x);
+    return mix(result_neg,result_pos,is_positive);
+}
+
+const int ZETA_PRECISION = 10;
+const float ZETA_COEFFICIENTS [10] = float[](1.0, 513.0, 23041.0, 263169.0, 1153025.0, 2306049.0, 2191361.0, 950273.0, 171537.0, 9217.0);
+
+vec2 czeta_main_branch(vec2 s){
+    vec2 sum = vec2(0.0, 0.0);
+    float sign = 1.0;
+    for(int n = 0; n < 10; ++n){
+        vec2 n_to_neg_s = cpow(vec2(float(n+1),0.0), cmult(s, MINUS_ONE));
+
+        n_to_neg_s = cmult(n_to_neg_s, vec2(sign * ZETA_COEFFICIENTS[ ZETA_PRECISION - 1 - n], 0.0));
+        sum = cadd(sum,n_to_neg_s);
+        sign = -sign;
+    }
+    vec2 eta = cmult(sum,vec2(1.0/ZETA_COEFFICIENTS[ZETA_PRECISION - 1], 0.0));
+    vec2 one_minus_s = csub(ONE,s);
+    vec2 two_pow = cpow(vec2(2.0,0.0),one_minus_s);
+    return cdiv(eta,csub(ONE,two_pow));
+}
+
+vec2 czeta_negative_branch(vec2 s){
+    vec2 one_minus_s = csub(ONE, s);
+    vec2 term1 = cpow(vec2(2.0, 0.0), s);
+    vec2 term2 = cpow(vec2(PI, 0.0), csub(s, ONE));
+    vec2 term3 = csin(cmult(s, vec2(PI / 2.0, 0.0)));
+    vec2 term4 = cgamma(one_minus_s);
+    vec2 term5 = czeta_main_branch(one_minus_s);
+    return cmult(term1, cmult(term2, cmult(term3, cmult(term4, term5))));
+
+}
+
+vec2 czeta(vec2 s){
+    if (s.x >= 0.0) {
+        return czeta_main_branch(s);
+    } 
+    else {
+        return czeta_negative_branch(s);
+    }
+}
+
+
+
 #define END_FUNCTION_DEFINITIONS HERE
 
 

src/types/types.h

@@ -65,6 +65,9 @@ enum Operator {
 	ASECH,
 	ACOTH,
 	CONJ,
+	ZETA,
+	GAMMA,
+	FACTORIAL,
 	RE,
 	IM,
 	DERIVATIVE,
@@ -155,6 +158,10 @@ inline std::vector<FullOperator> full_operators = {
 	{ {Arity::UNARY, Associativity::RIGHT, Operator::ASECH, "asech", 5}, 0, "casech", "ASECH"},
 	{ {Arity::UNARY, Associativity::RIGHT, Operator::ACOTH, "acoth", 5}, 0, "cacoth", "ACOTH"},
 
+	// NON ELEMANTARY
+	{ {Arity::UNARY, Associativity::RIGHT, Operator::GAMMA, "gamma", 5}, 0, "cgamma", "GAMMA"},
+	{ {Arity::UNARY, Associativity::RIGHT, Operator::ZETA, "zeta", 5}, 0, "czeta", "ZETA"},
+
 	//CONSTANTS
 	{ {Arity::NULLARY, Associativity::NONE, Operator::CONSTANT, "CONSTANT", 0, glm::vec2(0.0f)}, 0, "", "CONSTANT"},
 	{ {Arity::NULLARY, Associativity::NONE, Operator::CONSTANT, "i", 0, glm::vec2(0,1)}, 0, ""},
@@ -167,13 +174,14 @@ inline std::vector<FullOperator> full_operators = {
 	{ {Arity::NULLARY, Associativity::NONE, Operator::VARIABLEX, "x", 0}, 0, "vec2(z.x,0.0f)", "VARIABLEX"},
 	{ {Arity::NULLARY, Associativity::NONE, Operator::VARIABLEY, "y", 0}, 0, "vec2(z.y,0.0f)", "VARIABLEY"},
 	{ {Arity::NULLARY, Associativity::NONE, Operator::VARIABLET, "t", 0}, 0, "vec2(time,0.0f)", "VARIABLET"},
-	
+
 	//PARSER STUFF
 	{ {Arity::NULLARY, Associativity::NONE, Operator::VARIABLEPLACEHOLDER, "uplaceholder", 0}, 0},
 	{ {Arity::NULLARY, Associativity::NONE, Operator::SECONDVARIABLEPLACEHOLDER, "vplaceholder", 0}, 0},
 	{ {Arity::PAREN, Associativity::NONE, Operator::LPAREN, "(", 0}, 0, "" },
 	{ {Arity::PAREN, Associativity::NONE, Operator::RPAREN, ")", 0}, 0, "" },
 
+
 	//HIGHER ORDER
 	{ {Arity::UNARY, Associativity::RIGHT, Operator::DERIVATIVE, "derivative", 0,glm::vec2(0.0f),true}, 0, ""},