diff --git a/CHANGELOG.md b/CHANGELOG.md
index 96eeb249..54440da9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ## [Unreleased]
 
+## [3.1.6] - 2024-02-28
+### Fixed
+- Fix invalid calculation of `getPeriodicBezierInterpolation`
+
 ## [3.1.5] - 2024-02-27
 ### Fixed
 - Fix invalid calculation of `getBezierInterpolation`
diff --git a/package.json b/package.json
index 8eb902e7..5cba5dfb 100644
--- a/package.json
+++ b/package.json
@@ -1,6 +1,6 @@
 {
   "name": "okageo",
-  "version": "3.1.5",
+  "version": "3.1.6",
   "description": "parse SVG to polygons",
   "main": "./dist/okageo.js",
   "module": "./dist/okageo.mjs",
diff --git a/src/geo.ts b/src/geo.ts
index c9af6f45..7af5285a 100644
--- a/src/geo.ts
+++ b/src/geo.ts
@@ -1702,7 +1702,7 @@ export function getPeriodicBezierInterpolation(
   for (let i = 0; i < points.length - 2; i++) {
     B[i] = sub(multi(points[i + 1], 2), A[i + 1])
   }
-  B[points.length - 2] = sub(multi(points[0], 2), A[0])
+  B[points.length - 2] = sub(multi(points[points.length - 1], 2), A[0])
 
   return A.map((a, i) => [a, B[i]])
 }
@@ -1713,23 +1713,28 @@ export function getPeriodicBezierInterpolationA(points: IVec2[]): IVec2[] {
 
   const values: IVec2[] = []
   for (let i = 0; i < points.length - 1; i++) {
-    values.push(multi(add(multi(points[i], 2), points[i + 1]), 2))
+    values.push({
+      x: 4 * points[i].x + 2 * points[i + 1].x,
+      y: 4 * points[i].y + 2 * points[i + 1].y,
+    })
   }
   const y = solvePeriodicBezierInterpolationEquations(values, gamma)
 
-  const u = points.map(() => ({ x: 0, y: 0 }))
+  const u: IVec2[] = []
   u[0] = { x: gamma, y: gamma }
-  u[paramSize - 1] = { x: 1, y: 1 }
+  for (let i = 1; i < points.length - 2; i++) {
+    u.push({ x: 0, y: 0 })
+  }
+  u.push({ x: 1, y: 1 })
   const q = solvePeriodicBezierInterpolationEquations(u, gamma)
 
-  const v = u
   const vy = {
-    x: v[0].x * y[0].x + v[paramSize - 1].x * y[paramSize - 1].x,
-    y: v[0].y * y[0].y + v[paramSize - 1].y * y[paramSize - 1].y,
+    x: y[0].x + (1 / gamma) * y[paramSize - 1].x,
+    y: y[0].y + (1 / gamma) * y[paramSize - 1].y,
   }
   const vq = {
-    x: v[0].x * q[0].x + v[paramSize - 1].x * q[paramSize - 1].x,
-    y: v[0].y * q[0].y + v[paramSize - 1].y * q[paramSize - 1].y,
+    x: q[0].x + (1 / gamma) * q[paramSize - 1].x,
+    y: q[0].y + (1 / gamma) * q[paramSize - 1].y,
   }
 
   const A: IVec2[] = []
diff --git a/test/geo.test.ts b/test/geo.test.ts
index c8be1b31..9a5f9a71 100644
--- a/test/geo.test.ts
+++ b/test/geo.test.ts
@@ -2045,15 +2045,36 @@ describe('getPeriodicBezierInterpolation', () => {
     const ret0 = geo.getPeriodicBezierInterpolation(points)
     expect(ret0).toHaveLength(4)
     expect(ret0[0][0].x).toBeCloseTo(2.5)
-    expect(ret0[0][0].y).toBeCloseTo(-2.583)
+    expect(ret0[0][0].y).toBeCloseTo(-2.5)
     expect(ret0[0][1].x).toBeCloseTo(7.5)
-    expect(ret0[0][1].y).toBeCloseTo(-2.541)
+    expect(ret0[0][1].y).toBeCloseTo(-2.5)
     expect(ret0[1][0].x).toBeCloseTo(12.5)
-    expect(ret0[1][0].y).toBeCloseTo(2.542)
+    expect(ret0[1][0].y).toBeCloseTo(2.5)
     expect(ret0[1][1].x).toBeCloseTo(12.5)
-    expect(ret0[1][1].y).toBeCloseTo(7.583)
+    expect(ret0[1][1].y).toBeCloseTo(7.5)
     expect(ret0[3][1].x).toBeCloseTo(-2.5)
-    expect(ret0[3][1].y).toBeCloseTo(2.583)
+    expect(ret0[3][1].y).toBeCloseTo(2.5)
+  })
+
+  it('should return bezier control points: non-zero origin', () => {
+    const points = [
+      { x: 0, y: 0 },
+      { x: 10, y: 0 },
+      { x: 10, y: 10 },
+      { x: 0, y: 10 },
+      { x: 0, y: 0 },
+    ]
+    const ret0 = geo.getPeriodicBezierInterpolation(points)
+    const origin = { x: 100, y: 200 }
+    const shift = geo.getPeriodicBezierInterpolation(
+      points.map((p) => geo.add(p, origin))
+    )
+    ret0.forEach(([c1, c2], i) => {
+      expect(c1.x).toBeCloseTo(shift[i][0].x - origin.x)
+      expect(c1.y).toBeCloseTo(shift[i][0].y - origin.y)
+      expect(c2.x).toBeCloseTo(shift[i][1].x - origin.x)
+      expect(c2.y).toBeCloseTo(shift[i][1].y - origin.y)
+    })
   })
 })