summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/convex.c93
-rw-r--r--lib/convex.h21
-rw-r--r--lib/wolff.h1
3 files changed, 115 insertions, 0 deletions
diff --git a/lib/convex.c b/lib/convex.c
new file mode 100644
index 0000000..33af1b4
--- /dev/null
+++ b/lib/convex.c
@@ -0,0 +1,93 @@
+
+#include "convex.h"
+
+double slope(point_t *P, point_t *Q) {
+ return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x));
+}
+
+double *get_convex_minorant(uint64_t n, double *Gammas) {
+ list_t *L = (list_t *)calloc(1, sizeof(list_t));
+ L->p = (point_t *)calloc(1, sizeof(point_t));
+ L->p->x = 0;
+ L->p->y = Gammas[0];
+
+ list_t *pos = L;
+
+ for (uint64_t i = 1; i < n; i++) {
+ pos->next = (list_t *)calloc(1, sizeof(list_t));
+ pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+ pos->next->p->x = i;
+ pos->next->p->y = Gammas[i];
+ pos->next->prev = pos;
+ pos = pos->next;
+ }
+
+ pos->next = (list_t *)calloc(1, sizeof(list_t));
+ pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+ pos->next->p->x = n;
+ pos->next->p->y = 0;
+ pos->next->prev = pos;
+
+ list_t *X = L;
+ list_t *Y = L->next;
+ list_t *Z = Y->next;
+
+ while (true) {
+ if (slope(X->p, Y->p) <= slope(Y->p, Z->p)) {
+ X = Y;
+ Y = Z;
+ if (Z->next == NULL) {
+ break;
+ } else {
+ Z = Z->next;
+ }
+ } else {
+ Y->prev->next = Y->next;
+ Y->next->prev = Y->prev;
+ free(Y->p);
+ free(Y);
+ if (X->prev != NULL) {
+ Y = X;
+ X = X->prev;
+ } else {
+ if (Z->next != NULL) {
+ Y = Z;
+ Z = Z->next;
+ } else {
+ break;
+ }
+ }
+ }
+ }
+
+ pos = L;
+
+ double *g = (double *)calloc(n, sizeof(double));
+ double rho = 0;
+
+ for (uint64_t i = 0; i < n; i++) {
+ if (i > pos->next->p->x) {
+ pos = pos->next;
+ }
+
+ g[i] = pos->p->y + ((double)i - (double)(pos->p->x)) * (pos->next->p->y - pos->p->y) / ((double)(pos->next->p->x) - (double)(pos->p->x));
+
+ if (Gammas[i] - g[i] > rho) {
+ rho = Gammas[i] - g[i];
+ }
+ }
+
+ for (uint64_t i = 0; i < n; i++) {
+ g[i] += rho / 2;
+ }
+
+ while (L != NULL) {
+ free(L->p);
+ list_t *L_save = L;
+ L = L->next;
+ free(L);
+ }
+
+ return g;
+}
+
diff --git a/lib/convex.h b/lib/convex.h
new file mode 100644
index 0000000..29f4ae0
--- /dev/null
+++ b/lib/convex.h
@@ -0,0 +1,21 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef struct {
+ uint64_t x;
+ double y;
+} point_t;
+
+typedef struct list_tag {
+ struct list_tag *prev;
+ struct list_tag *next;
+ point_t *p;
+} list_t;
+
+double *get_convex_minorant(uint64_t n, double *Gammas);
+
diff --git a/lib/wolff.h b/lib/wolff.h
index 66ce493..dfe7b1c 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -17,6 +17,7 @@
#include <jst/rand.h>
#include "queue.h"
+#include "convex.h"
typedef enum {
WOLFF,