diff --git a/Tests/images/palette_negative.png b/Tests/images/palette_negative.png
index 7fcfd29a0db..2bce8d74952 100644
Binary files a/Tests/images/palette_negative.png and b/Tests/images/palette_negative.png differ
diff --git a/Tests/images/palette_sepia.png b/Tests/images/palette_sepia.png
index 9e7d6b0345b..77cec03de61 100644
Binary files a/Tests/images/palette_sepia.png and b/Tests/images/palette_sepia.png differ
diff --git a/Tests/images/palette_wedge.png b/Tests/images/palette_wedge.png
index 4b3d9ff3a21..613c34daa77 100644
Binary files a/Tests/images/palette_wedge.png and b/Tests/images/palette_wedge.png differ
diff --git a/Tests/test_image_entropy.py b/Tests/test_image_entropy.py
index c1dbb879b0b..a0230dff5c6 100644
--- a/Tests/test_image_entropy.py
+++ b/Tests/test_image_entropy.py
@@ -11,7 +11,7 @@ def entropy(mode: str) -> float:
     assert round(abs(entropy("L") - 7.063008716585465), 7) == 0
     assert round(abs(entropy("I") - 7.063008716585465), 7) == 0
     assert round(abs(entropy("F") - 7.063008716585465), 7) == 0
-    assert round(abs(entropy("P") - 5.082506854662517), 7) == 0
+    assert round(abs(entropy("P") - 5.085630558679487), 7) == 0
     assert round(abs(entropy("RGB") - 8.821286587714319), 7) == 0
     assert round(abs(entropy("RGBA") - 7.42724306524488), 7) == 0
     assert round(abs(entropy("CMYK") - 7.4272430652448795), 7) == 0
diff --git a/Tests/test_image_histogram.py b/Tests/test_image_histogram.py
index 436eb78a26e..864af4d8423 100644
--- a/Tests/test_image_histogram.py
+++ b/Tests/test_image_histogram.py
@@ -14,7 +14,7 @@ def histogram(mode: str) -> tuple[int, int, int]:
     assert histogram("La") == (512, 0, 16384)
     assert histogram("I") == (256, 0, 662)
     assert histogram("F") == (256, 0, 662)
-    assert histogram("P") == (256, 0, 1551)
+    assert histogram("P") == (256, 0, 1557)
     assert histogram("PA") == (512, 0, 16384)
     assert histogram("RGB") == (768, 4, 675)
     assert histogram("RGBA") == (1024, 0, 16384)
diff --git a/src/libImaging/Convert.c b/src/libImaging/Convert.c
index 330e5325c33..0c99ae1ccaf 100644
--- a/src/libImaging/Convert.c
+++ b/src/libImaging/Convert.c
@@ -1246,19 +1246,18 @@ topalette(
             /* Map each pixel to the nearest palette entry */
             ImagingSectionEnter(&cookie);
             for (y = 0; y < imIn->ysize; y++) {
-                int r, r0, r1, r2;
-                int g, g0, g1, g2;
-                int b, b0, b1, b2;
+                int r, r0, r1;
+                int g, g0, g1;
+                int b, b0, b1;
                 UINT8 *in = (UINT8 *)imIn->image[y];
                 UINT8 *out = alpha ? (UINT8 *)imOut->image32[y] : imOut->image8[y];
                 int *e = errors;
 
                 r = r0 = r1 = 0;
                 g = g0 = g1 = 0;
-                b = b0 = b1 = b2 = 0;
+                b = b0 = b1 = 0;
 
                 for (x = 0; x < imIn->xsize; x++, in += 4) {
-                    int d2;
                     INT16 *cache;
 
                     r = CLIP8(in[0] + (r + e[3 + 0]) / 16);
@@ -1281,38 +1280,28 @@ topalette(
                     g -= (int)palette->palette[cache[0] * 4 + 1];
                     b -= (int)palette->palette[cache[0] * 4 + 2];
 
-                    /* propagate errors (don't ask ;-) */
-                    r2 = r;
-                    d2 = r + r;
-                    r += d2;
-                    e[0] = r + r0;
-                    r += d2;
-                    r0 = r + r1;
-                    r1 = r2;
-                    r += d2;
-                    g2 = g;
-                    d2 = g + g;
-                    g += d2;
-                    e[1] = g + g0;
-                    g += d2;
-                    g0 = g + g1;
-                    g1 = g2;
-                    g += d2;
-                    b2 = b;
-                    d2 = b + b;
-                    b += d2;
-                    e[2] = b + b0;
-                    b += d2;
-                    b0 = b + b1;
-                    b1 = b2;
-                    b += d2;
+                    /* propagate errors */
+                    e[0] = 3 * r + r0;
+                    r0 = 5 * r + r1;
+                    r1 = r;
+                    r = 7 * r;
+
+                    e[1] = 3 * g + g0;
+                    g0 = 5 * g + g1;
+                    g1 = g;
+                    g = 7 * g;
+
+                    e[2] = 3 * b + b0;
+                    b0 = 5 * b + b1;
+                    b1 = b;
+                    b = 7 * b;
 
                     e += 3;
                 }
 
-                e[0] = b0;
-                e[1] = b1;
-                e[2] = b2;
+                e[0] = r0;
+                e[1] = g0;
+                e[2] = b0;
             }
             ImagingSectionLeave(&cookie);
             free(errors);
@@ -1385,7 +1374,7 @@ tobilevel(Imaging imOut, Imaging imIn) {
         /* map each pixel to black or white, using error diffusion */
         ImagingSectionEnter(&cookie);
         for (y = 0; y < imIn->ysize; y++) {
-            int l, l0, l1, l2, d2;
+            int l, l0, l1;
             UINT8 *in = (UINT8 *)imIn->image[y];
             UINT8 *out = imOut->image8[y];
 
@@ -1398,14 +1387,10 @@ tobilevel(Imaging imOut, Imaging imIn) {
 
                 /* propagate errors */
                 l -= (int)out[x];
-                l2 = l;
-                d2 = l + l;
-                l += d2;
-                errors[x] = l + l0;
-                l += d2;
-                l0 = l + l1;
-                l1 = l2;
-                l += d2;
+                errors[x] = 3 * l + l0;
+                l0 = 5 * l + l1;
+                l1 = l;
+                l = 7 * l;
             }
 
             errors[x] = l0;
@@ -1416,7 +1401,7 @@ tobilevel(Imaging imOut, Imaging imIn) {
         /* map each pixel to black or white, using error diffusion */
         ImagingSectionEnter(&cookie);
         for (y = 0; y < imIn->ysize; y++) {
-            int l, l0, l1, l2, d2;
+            int l, l0, l1;
             UINT8 *in = (UINT8 *)imIn->image[y];
             UINT8 *out = imOut->image8[y];
 
@@ -1429,14 +1414,10 @@ tobilevel(Imaging imOut, Imaging imIn) {
 
                 /* propagate errors */
                 l -= (int)out[x];
-                l2 = l;
-                d2 = l + l;
-                l += d2;
-                errors[x] = l + l0;
-                l += d2;
-                l0 = l + l1;
-                l1 = l2;
-                l += d2;
+                errors[x] = 3 * l + l0;
+                l0 = 5 * l + l1;
+                l1 = l;
+                l = 7 * l;
             }
 
             errors[x] = l0;