From 5d1d33d2d00b4509d0a9ca1084625e797bf46eae Mon Sep 17 00:00:00 2001
From: Wetitpig <winsto003@hotmail.com>
Date: Sat, 5 Sep 2020 14:06:09 +0800
Subject: [PATCH] Calculate on the go

---
 include/geodesic.h |  1 +
 src/math.c         | 13 +++++++++++++
 src/mpblock.c      | 15 +--------------
 src/vincenty.c     |  8 ++------
 4 files changed, 17 insertions(+), 20 deletions(-)

diff --git a/include/geodesic.h b/include/geodesic.h
index 499eaa1..26f49ba 100644
--- a/include/geodesic.h
+++ b/include/geodesic.h
@@ -28,6 +28,7 @@ struct Vector {
 };
 
 long double sqr(long double operand);
+long double double_fac(int x);
 long double atan2_modified(long double y, long double x);
 long double normalise_a(long double x);
 long double normalise_c(long double x);
diff --git a/src/math.c b/src/math.c
index 1e00cf7..8d5d193 100644
--- a/src/math.c
+++ b/src/math.c
@@ -6,6 +6,19 @@ long double sqr(long double operand)
 	return operand * operand;
 }
 
+long double double_fac(int x)
+{
+	if (x == -3)
+		return -1;
+	else {
+		long double y = 1;
+		while (x > 0) {
+			y *= x;
+			x -= 2;
+		}
+		return y;
+	}
+}
 
 long double atan2_modified(long double y, long double x)
 {
diff --git a/src/mpblock.c b/src/mpblock.c
index 36b8c79..610594c 100644
--- a/src/mpblock.c
+++ b/src/mpblock.c
@@ -1,21 +1,8 @@
 #include <math.h>
+#include "math.h"
 #include "io.h"
 #include "mpblock.h"
 
-long double double_fac(int x)
-{
-	if (x == -3)
-		return -1;
-	else {
-		long double y = 1;
-		while (x > 0) {
-			y *= x;
-			x -= 2;
-		}
-		return y;
-	}
-}
-
 int isblock(struct Coordinates *vertex)
 {
 	int block = 1;
diff --git a/src/vincenty.c b/src/vincenty.c
index b32786e..67c9e85 100644
--- a/src/vincenty.c
+++ b/src/vincenty.c
@@ -2,8 +2,6 @@
 #include "geodesic.h"
 #include "vincenty.h"
 
-const long double Alookup[5] = {25.0l/16384.0l, 1.0l/256.0l, 1.0l/64.0l, 1.0l/4.0l, 1.0l};
-
 const long double Blookup[8][4] = {
 {-1.0l/2.0l,3.0l/16.0l,-1.0l/32.0l,19.0l/2048.0l},
 {-1.0l/16.0l,1.0l/32.0l,-9.0l/2048.0l,7.0l/4096.0l},
@@ -64,10 +62,8 @@ long double helmertA(long double calp)
 	usq = calp * ECC2;
 	k1 = sqrtl(1.0l + usq);
 	k1 = (k1 - 1.0l) / (k1 + 1.0l);
-	for (k = 0; k < 5; k++) {
-		A *= sqr(k1);
-		A += Alookup[k];
-	}
+	for (k = 0; k < 11; k++)
+		A += sqr(double_fac(2 * k - 3) / double_fac(2 * k) * powl(k1, k));
 	A /= (1 - k1);
 	return A;
 }