From a5b99dd302cd8c84a28f2b45ba02002540129fd7 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Sun, 1 May 2016 17:31:54 +0300
Subject: [PATCH 01/21] move checks before mallocs to prevent memory leaks

---
 libImaging/Resample.c | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 64e90eed082..144ed6cba74 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -139,11 +139,11 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
     kmax = (int) ceil(support) * 2 + 1;
 
     // check for overflow
-    if (kmax > 0 && xsize > SIZE_MAX / kmax)
+    if (xsize > SIZE_MAX / (kmax * sizeof(float)))
         return (Imaging) ImagingError_MemoryError();
 
-    // sizeof(float) should be greater than 0
-    if (xsize * kmax > SIZE_MAX / sizeof(float))
+    // sizeof(int) should be greater than 0 as well
+    if (xsize > SIZE_MAX / (2 * sizeof(int)))
         return (Imaging) ImagingError_MemoryError();
 
     /* coefficient buffer */
@@ -151,10 +151,6 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
     if ( ! kk)
         return (Imaging) ImagingError_MemoryError();
 
-    // sizeof(int) should be greater than 0 as well
-    if (xsize > SIZE_MAX / (2 * sizeof(int)))
-        return (Imaging) ImagingError_MemoryError();
-
     xbounds = malloc(xsize * 2 * sizeof(int));
     if ( ! xbounds) {
         free(kk);

From 0e2a8e43232ad93eef0b26971fbf8687472e8d2d Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Sun, 1 May 2016 17:48:13 +0300
Subject: [PATCH 02/21] intermediate not normalized buffer for coefficients

---
 libImaging/Resample.c | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 144ed6cba74..26740281f7f 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -104,7 +104,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
     float center, ww, ss, ss0, ss1, ss2, ss3;
     int xx, yy, x, kmax, xmin, xmax;
     int *xbounds;
-    float *k, *kk;
+    float *k, *kk, *kw;
 
     /* check filter */
     switch (filter) {
@@ -151,14 +151,21 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
     if ( ! kk)
         return (Imaging) ImagingError_MemoryError();
 
+    /* intermediate not normalized buffer for coefficients */ 
+    kw = malloc(kmax * sizeof(float));
+    if ( ! kw) {
+        free(kk);
+        return (Imaging) ImagingError_MemoryError();
+    }
+
     xbounds = malloc(xsize * 2 * sizeof(int));
     if ( ! xbounds) {
         free(kk);
+        free(kw);
         return (Imaging) ImagingError_MemoryError();
     }
 
     for (xx = 0; xx < xsize; xx++) {
-        k = &kk[xx * kmax];
         center = (xx + 0.5) * scale;
         ww = 0.0;
         ss = 1.0 / filterscale;
@@ -169,18 +176,21 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
         if (xmax > imIn->xsize)
             xmax = imIn->xsize;
         for (x = xmin; x < xmax; x++) {
-            float w = filterp->filter((x - center + 0.5) * ss) * ss;
-            k[x - xmin] = w;
+            float w = filterp->filter((x - center + 0.5) * ss);
+            kw[x - xmin] = w;
             ww += w;
         }
+        k = &kk[xx * kmax];
         for (x = 0; x < xmax - xmin; x++) {
             if (ww != 0.0)
-                k[x] /= ww;
+                k[x] = kw[x] / ww;
         }
         xbounds[xx * 2 + 0] = xmin;
         xbounds[xx * 2 + 1] = xmax;
     }
 
+    free(kw);
+
     imOut = ImagingNew(imIn->mode, xsize, imIn->ysize);
     if ( ! imOut) {
         free(kk);

From 3ee407e8590de584e1607c8bbd2ffd207aa42bdc Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Sun, 1 May 2016 17:54:26 +0300
Subject: [PATCH 03/21] pass struct filter to ImagingResampleHorizontal

---
 libImaging/Resample.c | 42 +++++++++++++++++++++---------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 26740281f7f..edb5c5ba525 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -95,34 +95,16 @@ static float inline i2f(int v) { return (float) v; }
 
 
 Imaging
-ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
+ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
 {
     ImagingSectionCookie cookie;
     Imaging imOut;
-    struct filter *filterp;
     float support, scale, filterscale;
     float center, ww, ss, ss0, ss1, ss2, ss3;
     int xx, yy, x, kmax, xmin, xmax;
     int *xbounds;
     float *k, *kk, *kw;
 
-    /* check filter */
-    switch (filter) {
-    case IMAGING_TRANSFORM_LANCZOS:
-        filterp = &LANCZOS;
-        break;
-    case IMAGING_TRANSFORM_BILINEAR:
-        filterp = &BILINEAR;
-        break;
-    case IMAGING_TRANSFORM_BICUBIC:
-        filterp = &BICUBIC;
-        break;
-    default:
-        return (Imaging) ImagingError_ValueError(
-            "unsupported resampling filter"
-            );
-    }
-
     /* prepare for horizontal stretch */
     filterscale = scale = (float) imIn->xsize / xsize;
 
@@ -302,6 +284,7 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
 {
     Imaging imTemp1, imTemp2, imTemp3;
     Imaging imOut;
+    struct filter *filterp;
 
     if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
         return (Imaging) ImagingError_ModeError();
@@ -309,8 +292,25 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
     if (imIn->type == IMAGING_TYPE_SPECIAL)
         return (Imaging) ImagingError_ModeError();
 
+    /* check filter */
+    switch (filter) {
+    case IMAGING_TRANSFORM_LANCZOS:
+        filterp = &LANCZOS;
+        break;
+    case IMAGING_TRANSFORM_BILINEAR:
+        filterp = &BILINEAR;
+        break;
+    case IMAGING_TRANSFORM_BICUBIC:
+        filterp = &BICUBIC;
+        break;
+    default:
+        return (Imaging) ImagingError_ValueError(
+            "unsupported resampling filter"
+            );
+    }
+
     /* two-pass resize, first pass */
-    imTemp1 = ImagingResampleHorizontal(imIn, xsize, filter);
+    imTemp1 = ImagingResampleHorizontal(imIn, xsize, filterp);
     if ( ! imTemp1)
         return NULL;
 
@@ -321,7 +321,7 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
         return NULL;
 
     /* second pass */
-    imTemp3 = ImagingResampleHorizontal(imTemp2, ysize, filter);
+    imTemp3 = ImagingResampleHorizontal(imTemp2, ysize, filterp);
     ImagingDelete(imTemp2);
     if ( ! imTemp3)
         return NULL;

From 342e6f213b9ed6ff2b37dc5b435d17af9c1fc1cb Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Tue, 3 May 2016 11:31:13 +0200
Subject: [PATCH 04/21] extra line

---
 libImaging/Resample.c | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index edb5c5ba525..8e50c4b1b56 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -107,15 +107,12 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
 
     /* prepare for horizontal stretch */
     filterscale = scale = (float) imIn->xsize / xsize;
-
-    /* determine support size (length of resampling filter) */
-    support = filterp->support;
-
     if (filterscale < 1.0) {
         filterscale = 1.0;
     }
 
-    support = support * filterscale;
+    /* determine support size (length of resampling filter) */
+    support = filterp->support * filterscale;
 
     /* maximum number of coofs */
     kmax = (int) ceil(support) * 2 + 1;

From e13297b396adc24a383862bef14e6fe0023251d8 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Tue, 3 May 2016 12:23:16 +0200
Subject: [PATCH 05/21] make coefficients ints

---
 libImaging/Resample.c | 82 +++++++++++++++++++------------------------
 1 file changed, 37 insertions(+), 45 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 8e50c4b1b56..0098ab4927d 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -38,8 +38,6 @@ static inline float lanczos_filter(float x)
     return 0.0;
 }
 
-static struct filter LANCZOS = { lanczos_filter, 3.0 };
-
 static inline float bilinear_filter(float x)
 {
     if (x < 0.0)
@@ -49,8 +47,6 @@ static inline float bilinear_filter(float x)
     return 0.0;
 }
 
-static struct filter BILINEAR = { bilinear_filter, 1.0 };
-
 static inline float bicubic_filter(float x)
 {
     /* https://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm */
@@ -65,33 +61,27 @@ static inline float bicubic_filter(float x)
 #undef a
 }
 
+static struct filter LANCZOS = { lanczos_filter, 3.0 };
+static struct filter BILINEAR = { bilinear_filter, 1.0 };
 static struct filter BICUBIC = { bicubic_filter, 2.0 };
 
 
-static inline UINT8 clip8(float in)
+
+/* 8 bits for result. Filter can have negative areas.
+   In one cases the sum of the coefficients will be negative,
+   in the other it will be more than 1.0. That is why we need
+   two extra bits for overflow and int type. */
+#define PRECISION_BITS (32 - 8 - 2)
+
+
+static inline UINT8 clip8(int in)
 {
-    int out = (int) in;
-    if (out >= 255)
+    if (in >= (1 << PRECISION_BITS << 8))
        return 255;
-    if (out <= 0)
+    if (in <= 0)
         return 0;
-    return (UINT8) out;
-}
-
-
-/* This is work around bug in GCC prior 4.9 in 64-bit mode.
-   GCC generates code with partial dependency which 3 times slower.
-   See: http://stackoverflow.com/a/26588074/253146 */
-#if defined(__x86_64__) && defined(__SSE__) &&  ! defined(__NO_INLINE__) && \
-    ! defined(__clang__) && defined(GCC_VERSION) && (GCC_VERSION < 40900)
-static float __attribute__((always_inline)) i2f(int v) {
-    float x;
-    __asm__("xorps %0, %0; cvtsi2ss %1, %0" : "=X"(x) : "r"(v) );
-    return x;
+    return (UINT8) (in >> PRECISION_BITS);
 }
-#else
-static float inline i2f(int v) { return (float) v; }
-#endif
 
 
 Imaging
@@ -100,10 +90,12 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
     ImagingSectionCookie cookie;
     Imaging imOut;
     float support, scale, filterscale;
-    float center, ww, ss, ss0, ss1, ss2, ss3;
+    float center, ww, ss;
+    int ss0, ss1, ss2, ss3;
     int xx, yy, x, kmax, xmin, xmax;
     int *xbounds;
-    float *k, *kk, *kw;
+    int *k, *kk;
+    float *kw;
 
     /* prepare for horizontal stretch */
     filterscale = scale = (float) imIn->xsize / xsize;
@@ -118,7 +110,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
     kmax = (int) ceil(support) * 2 + 1;
 
     // check for overflow
-    if (xsize > SIZE_MAX / (kmax * sizeof(float)))
+    if (xsize > SIZE_MAX / (kmax * sizeof(int)))
         return (Imaging) ImagingError_MemoryError();
 
     // sizeof(int) should be greater than 0 as well
@@ -126,7 +118,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
         return (Imaging) ImagingError_MemoryError();
 
     /* coefficient buffer */
-    kk = malloc(xsize * kmax * sizeof(float));
+    kk = malloc(xsize * kmax * sizeof(int));
     if ( ! kk)
         return (Imaging) ImagingError_MemoryError();
 
@@ -162,7 +154,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
         k = &kk[xx * kmax];
         for (x = 0; x < xmax - xmin; x++) {
             if (ww != 0.0)
-                k[x] = kw[x] / ww;
+                k[x] = (int) floor(0.5 + kw[x] / ww * (1 << PRECISION_BITS));
         }
         xbounds[xx * 2 + 0] = xmin;
         xbounds[xx * 2 + 1] = xmax;
@@ -186,10 +178,10 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                 xmin = xbounds[xx * 2 + 0];
                 xmax = xbounds[xx * 2 + 1];
                 k = &kk[xx * kmax];
-                ss = 0.5;
+                ss0 = 0;
                 for (x = xmin; x < xmax; x++)
-                    ss += i2f(imIn->image8[yy][x]) * k[x - xmin];
-                imOut->image8[yy][xx] = clip8(ss);
+                    ss0 += ((UINT8) imIn->image8[yy][x]) * k[x - xmin];
+                imOut->image8[yy][xx] = clip8(ss0);
             }
         } else {
             switch(imIn->type) {
@@ -200,10 +192,10 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = 0.5;
+                        ss0 = ss1 = 0;
                         for (x = xmin; x < xmax; x++) {
-                            ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += i2f((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
+                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                            ss1 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
                         }
                         imOut->image[yy][xx*4 + 0] = clip8(ss0);
                         imOut->image[yy][xx*4 + 3] = clip8(ss1);
@@ -213,11 +205,11 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = 0.5;
+                        ss0 = ss1 = ss2 = 0;
                         for (x = xmin; x < xmax; x++) {
-                            ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += i2f((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                            ss2 += i2f((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
+                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                            ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
+                            ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
                         }
                         imOut->image[yy][xx*4 + 0] = clip8(ss0);
                         imOut->image[yy][xx*4 + 1] = clip8(ss1);
@@ -228,12 +220,12 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = ss3 = 0.5;
+                        ss0 = ss1 = ss2 = ss3 = 0;
                         for (x = xmin; x < xmax; x++) {
-                            ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += i2f((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                            ss2 += i2f((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
-                            ss3 += i2f((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
+                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                            ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
+                            ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
+                            ss3 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
                         }
                         imOut->image[yy][xx*4 + 0] = clip8(ss0);
                         imOut->image[yy][xx*4 + 1] = clip8(ss1);
@@ -250,7 +242,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                     k = &kk[xx * kmax];
                     ss = 0.0;
                     for (x = xmin; x < xmax; x++)
-                        ss += i2f(IMAGING_PIXEL_I(imIn, x, yy)) * k[x - xmin];
+                        ss += (IMAGING_PIXEL_I(imIn, x, yy)) * k[x - xmin];
                     IMAGING_PIXEL_I(imOut, xx, yy) = (int) ss;
                 }
                 break;

From 9618ec98b86929bb521c250e02dd4b08fece0f1f Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Tue, 3 May 2016 13:22:51 +0200
Subject: [PATCH 06/21] test for consistency

---
 Tests/test_image_resample.py | 17 +++++++++++++++++
 libImaging/Resample.c        |  8 ++++----
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 79fc3ffaf74..62e38e002eb 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -150,5 +150,22 @@ def test_enlarge_lanczos(self):
             self.make_sample(data, (12, 12)))
 
 
+class CoreResampleConsistencyTest(PillowTestCase):
+
+    def test_8u(self):
+        im = Image.new('RGB', (512, 9), (0, 64, 255))
+        im = im.resize((9, 512), Image.LANCZOS)
+        r, g, b = im.split()
+
+        for channel, color in [(r, 0), (g, 64), (b, 255)]:
+            px = channel.load()
+            for x in range(channel.size[0]):
+                for y in range(channel.size[1]):
+                    if px[x, y] != color:
+                        message = "{} != {} for pixel {}".format(
+                            px[x, y], color, (x, y))
+                        self.assertEqual(px[x, y], color, message)
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 0098ab4927d..5aea0bb6628 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -178,7 +178,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                 xmin = xbounds[xx * 2 + 0];
                 xmax = xbounds[xx * 2 + 1];
                 k = &kk[xx * kmax];
-                ss0 = 0;
+                ss0 = 1 << (PRECISION_BITS -1);
                 for (x = xmin; x < xmax; x++)
                     ss0 += ((UINT8) imIn->image8[yy][x]) * k[x - xmin];
                 imOut->image8[yy][xx] = clip8(ss0);
@@ -192,7 +192,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = 0;
+                        ss0 = ss1 = 1 << (PRECISION_BITS -1);
                         for (x = xmin; x < xmax; x++) {
                             ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
                             ss1 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
@@ -205,7 +205,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = 0;
+                        ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
                         for (x = xmin; x < xmax; x++) {
                             ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
                             ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
@@ -220,7 +220,7 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                         xmin = xbounds[xx * 2 + 0];
                         xmax = xbounds[xx * 2 + 1];
                         k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = ss3 = 0;
+                        ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
                         for (x = xmin; x < xmax; x++) {
                             ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
                             ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];

From e47002dec3c3f926cf7288f2ea6e7890ed6a56c4 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 11:53:43 +0300
Subject: [PATCH 07/21] tests for 32bit modes

---
 Tests/test_image_resample.py | 30 +++++++++++++++++++++++++-----
 1 file changed, 25 insertions(+), 5 deletions(-)

diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 62e38e002eb..7d08f7ec4a6 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -152,12 +152,12 @@ def test_enlarge_lanczos(self):
 
 class CoreResampleConsistencyTest(PillowTestCase):
 
-    def test_8u(self):
-        im = Image.new('RGB', (512, 9), (0, 64, 255))
-        im = im.resize((9, 512), Image.LANCZOS)
-        r, g, b = im.split()
+    def make_case(self, mode, fill):
+        im = Image.new(mode, (512, 9), fill)
+        return (im.resize((9, 512), Image.LANCZOS), im.load()[0, 0])
 
-        for channel, color in [(r, 0), (g, 64), (b, 255)]:
+    def run_cases(self, cases):
+        for channel, color in cases:
             px = channel.load()
             for x in range(channel.size[0]):
                 for y in range(channel.size[1]):
@@ -166,6 +166,26 @@ def test_8u(self):
                             px[x, y], color, (x, y))
                         self.assertEqual(px[x, y], color, message)
 
+    def test_8u(self):
+        im, color = self.make_case('RGB', (0, 64, 255))
+        self.run_cases(zip(im.split(), color))
+
+    def test_32i(self):
+        self.run_cases([
+            self.make_case('I', 12),
+            self.make_case('I', 0x7fffffff),
+            self.make_case('I', -12),
+            self.make_case('I', -1 << 31),
+        ])
+
+    def test_32f(self):
+        self.run_cases([
+            self.make_case('F', 1),
+            self.make_case('F', 3.40282306074e+38),
+            self.make_case('F', 1.175494e-38),
+            self.make_case('F', 1.192093e-07),
+        ])
+
 
 if __name__ == '__main__':
     unittest.main()

From 7693da9287f324d5f3a9d072aafef0a64746717f Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 12:53:22 +0300
Subject: [PATCH 08/21] separate 8 and 32 bpc implementations

---
 libImaging/Resample.c | 206 +++++++++++++++++++++++++++++++-----------
 1 file changed, 151 insertions(+), 55 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 5aea0bb6628..7bef3ad8a8e 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -85,7 +85,7 @@ static inline UINT8 clip8(int in)
 
 
 Imaging
-ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
+ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
 {
     ImagingSectionCookie cookie;
     Imaging imOut;
@@ -183,57 +183,140 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                     ss0 += ((UINT8) imIn->image8[yy][x]) * k[x - xmin];
                 imOut->image8[yy][xx] = clip8(ss0);
             }
-        } else {
-            switch(imIn->type) {
-            case IMAGING_TYPE_UINT8:
-                /* n-bit grayscale */
-                if (imIn->bands == 2) {
-                    for (xx = 0; xx < xsize; xx++) {
-                        xmin = xbounds[xx * 2 + 0];
-                        xmax = xbounds[xx * 2 + 1];
-                        k = &kk[xx * kmax];
-                        ss0 = ss1 = 1 << (PRECISION_BITS -1);
-                        for (x = xmin; x < xmax; x++) {
-                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
-                        }
-                        imOut->image[yy][xx*4 + 0] = clip8(ss0);
-                        imOut->image[yy][xx*4 + 3] = clip8(ss1);
+        } else if (imIn->type == IMAGING_TYPE_UINT8) {
+            /* n-bit grayscale */
+            if (imIn->bands == 2) {
+                for (xx = 0; xx < xsize; xx++) {
+                    xmin = xbounds[xx * 2 + 0];
+                    xmax = xbounds[xx * 2 + 1];
+                    k = &kk[xx * kmax];
+                    ss0 = ss1 = 1 << (PRECISION_BITS -1);
+                    for (x = xmin; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                        ss1 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
                     }
-                } else if (imIn->bands == 3) {
-                    for (xx = 0; xx < xsize; xx++) {
-                        xmin = xbounds[xx * 2 + 0];
-                        xmax = xbounds[xx * 2 + 1];
-                        k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
-                        for (x = xmin; x < xmax; x++) {
-                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                            ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
-                        }
-                        imOut->image[yy][xx*4 + 0] = clip8(ss0);
-                        imOut->image[yy][xx*4 + 1] = clip8(ss1);
-                        imOut->image[yy][xx*4 + 2] = clip8(ss2);
+                    imOut->image[yy][xx*4 + 0] = clip8(ss0);
+                    imOut->image[yy][xx*4 + 3] = clip8(ss1);
+                }
+            } else if (imIn->bands == 3) {
+                for (xx = 0; xx < xsize; xx++) {
+                    xmin = xbounds[xx * 2 + 0];
+                    xmax = xbounds[xx * 2 + 1];
+                    k = &kk[xx * kmax];
+                    ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
+                    for (x = xmin; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                        ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
+                        ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
                     }
-                } else {
-                    for (xx = 0; xx < xsize; xx++) {
-                        xmin = xbounds[xx * 2 + 0];
-                        xmax = xbounds[xx * 2 + 1];
-                        k = &kk[xx * kmax];
-                        ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
-                        for (x = xmin; x < xmax; x++) {
-                            ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                            ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                            ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
-                            ss3 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
-                        }
-                        imOut->image[yy][xx*4 + 0] = clip8(ss0);
-                        imOut->image[yy][xx*4 + 1] = clip8(ss1);
-                        imOut->image[yy][xx*4 + 2] = clip8(ss2);
-                        imOut->image[yy][xx*4 + 3] = clip8(ss3);
+                    imOut->image[yy][xx*4 + 0] = clip8(ss0);
+                    imOut->image[yy][xx*4 + 1] = clip8(ss1);
+                    imOut->image[yy][xx*4 + 2] = clip8(ss2);
+                }
+            } else {
+                for (xx = 0; xx < xsize; xx++) {
+                    xmin = xbounds[xx * 2 + 0];
+                    xmax = xbounds[xx * 2 + 1];
+                    k = &kk[xx * kmax];
+                    ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
+                    for (x = xmin; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
+                        ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
+                        ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
+                        ss3 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
                     }
+                    imOut->image[yy][xx*4 + 0] = clip8(ss0);
+                    imOut->image[yy][xx*4 + 1] = clip8(ss1);
+                    imOut->image[yy][xx*4 + 2] = clip8(ss2);
+                    imOut->image[yy][xx*4 + 3] = clip8(ss3);
                 }
-                break;
+            }
+        }
+    }
+    ImagingSectionLeave(&cookie);
+    free(kk);
+    free(xbounds);
+    return imOut;
+}
+
+
+Imaging
+ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
+{
+    ImagingSectionCookie cookie;
+    Imaging imOut;
+    double support, scale, filterscale;
+    double center, ww, ss;
+    int xx, yy, x, kmax, xmin, xmax;
+    int *xbounds;
+    double *k, *kk;
+
+    /* prepare for horizontal stretch */
+    filterscale = scale = (float) imIn->xsize / xsize;
+    if (filterscale < 1.0) {
+        filterscale = 1.0;
+    }
+
+    /* determine support size (length of resampling filter) */
+    support = filterp->support * filterscale;
+
+    /* maximum number of coofs */
+    kmax = (int) ceil(support) * 2 + 1;
+
+    // check for overflow
+    if (xsize > SIZE_MAX / (kmax * sizeof(double)))
+        return (Imaging) ImagingError_MemoryError();
+
+    // sizeof(double) should be greater than 0 as well
+    if (xsize > SIZE_MAX / (2 * sizeof(double)))
+        return (Imaging) ImagingError_MemoryError();
+
+    /* coefficient buffer */
+    kk = malloc(xsize * kmax * sizeof(double));
+    if ( ! kk)
+        return (Imaging) ImagingError_MemoryError();
+
+    xbounds = malloc(xsize * 2 * sizeof(int));
+    if ( ! xbounds) {
+        free(kk);
+        return (Imaging) ImagingError_MemoryError();
+    }
+
+    for (xx = 0; xx < xsize; xx++) {
+        k = &kk[xx * kmax];
+        center = (xx + 0.5) * scale;
+        ww = 0.0;
+        ss = 1.0 / filterscale;
+        xmin = (int) floor(center - support);
+        if (xmin < 0)
+            xmin = 0;
+        xmax = (int) ceil(center + support);
+        if (xmax > imIn->xsize)
+            xmax = imIn->xsize;
+        for (x = xmin; x < xmax; x++) {
+            double w = filterp->filter((x - center + 0.5) * ss);
+            k[x - xmin] = w;
+            ww += w;
+        }
+        for (x = 0; x < xmax - xmin; x++) {
+            if (ww != 0.0)
+                k[x] /= ww;
+        }
+        xbounds[xx * 2 + 0] = xmin;
+        xbounds[xx * 2 + 1] = xmax;
+    }
+
+    imOut = ImagingNew(imIn->mode, xsize, imIn->ysize);
+    if ( ! imOut) {
+        free(kk);
+        free(xbounds);
+        return NULL;
+    }
+
+    ImagingSectionEnter(&cookie);
+    /* horizontal stretch */
+    for (yy = 0; yy < imOut->ysize; yy++) {
+        switch(imIn->type) {
             case IMAGING_TYPE_INT32:
                 /* 32-bit integer */
                 for (xx = 0; xx < xsize; xx++) {
@@ -242,8 +325,8 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                     k = &kk[xx * kmax];
                     ss = 0.0;
                     for (x = xmin; x < xmax; x++)
-                        ss += (IMAGING_PIXEL_I(imIn, x, yy)) * k[x - xmin];
-                    IMAGING_PIXEL_I(imOut, xx, yy) = (int) ss;
+                        ss += IMAGING_PIXEL_I(imIn, x, yy) * k[x - xmin];
+                    IMAGING_PIXEL_I(imOut, xx, yy) = lround(ss);
                 }
                 break;
             case IMAGING_TYPE_FLOAT32:
@@ -258,7 +341,6 @@ ImagingResampleHorizontal(Imaging imIn, int xsize, struct filter *filterp)
                     IMAGING_PIXEL_F(imOut, xx, yy) = ss;
                 }
                 break;
-            }
         }
     }
     ImagingSectionLeave(&cookie);
@@ -274,12 +356,26 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
     Imaging imTemp1, imTemp2, imTemp3;
     Imaging imOut;
     struct filter *filterp;
+    Imaging (*ResampleHorizontal)(Imaging imIn, int xsize, struct filter *filterp);
 
     if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
         return (Imaging) ImagingError_ModeError();
 
-    if (imIn->type == IMAGING_TYPE_SPECIAL)
-        return (Imaging) ImagingError_ModeError();
+    if (imIn->image8) {
+        ResampleHorizontal = ImagingResampleHorizontal_8bpc;
+    } else {
+        switch(imIn->type) {
+            case IMAGING_TYPE_UINT8:
+                ResampleHorizontal = ImagingResampleHorizontal_8bpc;
+                break;
+            case IMAGING_TYPE_INT32:
+            case IMAGING_TYPE_FLOAT32:
+                ResampleHorizontal = ImagingResampleHorizontal_32bpc;
+                break;
+            default:
+                return (Imaging) ImagingError_ModeError();
+        }
+    }
 
     /* check filter */
     switch (filter) {
@@ -299,7 +395,7 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
     }
 
     /* two-pass resize, first pass */
-    imTemp1 = ImagingResampleHorizontal(imIn, xsize, filterp);
+    imTemp1 = ResampleHorizontal(imIn, xsize, filterp);
     if ( ! imTemp1)
         return NULL;
 
@@ -310,7 +406,7 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
         return NULL;
 
     /* second pass */
-    imTemp3 = ImagingResampleHorizontal(imTemp2, ysize, filterp);
+    imTemp3 = ResampleHorizontal(imTemp2, ysize, filterp);
     ImagingDelete(imTemp2);
     if ( ! imTemp3)
         return NULL;

From 023f7adcf6377a3434f88f3d2cbda167a30d04be Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 13:15:54 +0300
Subject: [PATCH 09/21] precompute function

---
 libImaging/Resample.c | 69 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 69 insertions(+)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 7bef3ad8a8e..5d94c29e82c 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -84,6 +84,75 @@ static inline UINT8 clip8(int in)
 }
 
 
+int
+ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
+                  int **xboundsp, double **kkp) {
+    double support, scale, filterscale;
+    double center, ww, ss;
+    int xx, x, kmax, xmin, xmax;
+    int *xbounds;
+    double *kk, *k;
+
+    /* prepare for horizontal stretch */
+    filterscale = scale = (float) inSize / outSize;
+    if (filterscale < 1.0) {
+        filterscale = 1.0;
+    }
+
+    /* determine support size (length of resampling filter) */
+    support = filterp->support * filterscale;
+
+    /* maximum number of coofs */
+    kmax = (int) ceil(support) * 2 + 1;
+
+    // check for overflow
+    if (outSize > SIZE_MAX / (kmax * sizeof(double)))
+        return (int) (size_t) ImagingError_MemoryError();
+
+    // sizeof(double) should be greater than 0 as well
+    if (outSize > SIZE_MAX / (2 * sizeof(double)))
+        return (int) (size_t) ImagingError_MemoryError();
+
+    /* coefficient buffer */
+    kk = malloc(outSize * kmax * sizeof(double));
+    if ( ! kk)
+        return (int) (size_t) ImagingError_MemoryError();
+
+    xbounds = malloc(outSize * 2 * sizeof(int));
+    if ( ! xbounds) {
+        free(kk);
+        return (int) (size_t) ImagingError_MemoryError();
+    }
+
+    for (xx = 0; xx < outSize; xx++) {
+        center = (xx + 0.5) * scale;
+        ww = 0.0;
+        ss = 1.0 / filterscale;
+        xmin = (int) floor(center - support);
+        if (xmin < 0)
+            xmin = 0;
+        xmax = (int) ceil(center + support);
+        if (xmax > inSize)
+            xmax = inSize;
+        k = &kk[xx * kmax];
+        for (x = xmin; x < xmax; x++) {
+            double w = filterp->filter((x - center + 0.5) * ss);
+            k[x - xmin] = w;
+            ww += w;
+        }
+        for (x = 0; x < xmax - xmin; x++) {
+            if (ww != 0.0)
+                k[x] /= ww;
+        }
+        xbounds[xx * 2 + 0] = xmin;
+        xbounds[xx * 2 + 1] = xmax;
+    }
+    *xboundsp = xbounds;
+    *kkp = kk;
+    return kmax;
+}
+
+
 Imaging
 ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
 {

From 04552b015df7b158dbf754eb9474c42b1016ee2b Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 18:03:20 +0300
Subject: [PATCH 10/21] use ImagingPrecompute in
 ImagingResampleHorizontal_32bpc

---
 libImaging/Resample.c | 80 ++++---------------------------------------
 1 file changed, 7 insertions(+), 73 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 5d94c29e82c..def135524e4 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -1,18 +1,3 @@
-/*
- * The Python Imaging Library
- * $Id$
- *
- * Pillow image resampling support
- *
- * history:
- * 2002-03-09 fl  Created (for PIL 1.1.3)
- * 2002-03-10 fl  Added support for mode "F"
- *
- * Copyright (c) 1997-2002 by Secret Labs AB
- *
- * See the README file for information on usage and redistribution.
- */
-
 #include "Imaging.h"
 
 #include <math.h>
@@ -107,21 +92,21 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
 
     // check for overflow
     if (outSize > SIZE_MAX / (kmax * sizeof(double)))
-        return (int) (size_t) ImagingError_MemoryError();
+        return 0;
 
     // sizeof(double) should be greater than 0 as well
     if (outSize > SIZE_MAX / (2 * sizeof(double)))
-        return (int) (size_t) ImagingError_MemoryError();
+        return 0;
 
     /* coefficient buffer */
     kk = malloc(outSize * kmax * sizeof(double));
     if ( ! kk)
-        return (int) (size_t) ImagingError_MemoryError();
+        return 0;
 
     xbounds = malloc(outSize * 2 * sizeof(int));
     if ( ! xbounds) {
         free(kk);
-        return (int) (size_t) ImagingError_MemoryError();
+        return 0;
     }
 
     for (xx = 0; xx < outSize; xx++) {
@@ -314,67 +299,16 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
 {
     ImagingSectionCookie cookie;
     Imaging imOut;
-    double support, scale, filterscale;
-    double center, ww, ss;
+    double ss;
     int xx, yy, x, kmax, xmin, xmax;
     int *xbounds;
     double *k, *kk;
 
-    /* prepare for horizontal stretch */
-    filterscale = scale = (float) imIn->xsize / xsize;
-    if (filterscale < 1.0) {
-        filterscale = 1.0;
-    }
-
-    /* determine support size (length of resampling filter) */
-    support = filterp->support * filterscale;
-
-    /* maximum number of coofs */
-    kmax = (int) ceil(support) * 2 + 1;
-
-    // check for overflow
-    if (xsize > SIZE_MAX / (kmax * sizeof(double)))
-        return (Imaging) ImagingError_MemoryError();
-
-    // sizeof(double) should be greater than 0 as well
-    if (xsize > SIZE_MAX / (2 * sizeof(double)))
-        return (Imaging) ImagingError_MemoryError();
-
-    /* coefficient buffer */
-    kk = malloc(xsize * kmax * sizeof(double));
-    if ( ! kk)
-        return (Imaging) ImagingError_MemoryError();
-
-    xbounds = malloc(xsize * 2 * sizeof(int));
-    if ( ! xbounds) {
-        free(kk);
+    kmax = ImagingPrecompute(imIn->xsize, xsize, filterp, &xbounds, &kk);
+    if ( ! kmax) {
         return (Imaging) ImagingError_MemoryError();
     }
 
-    for (xx = 0; xx < xsize; xx++) {
-        k = &kk[xx * kmax];
-        center = (xx + 0.5) * scale;
-        ww = 0.0;
-        ss = 1.0 / filterscale;
-        xmin = (int) floor(center - support);
-        if (xmin < 0)
-            xmin = 0;
-        xmax = (int) ceil(center + support);
-        if (xmax > imIn->xsize)
-            xmax = imIn->xsize;
-        for (x = xmin; x < xmax; x++) {
-            double w = filterp->filter((x - center + 0.5) * ss);
-            k[x - xmin] = w;
-            ww += w;
-        }
-        for (x = 0; x < xmax - xmin; x++) {
-            if (ww != 0.0)
-                k[x] /= ww;
-        }
-        xbounds[xx * 2 + 0] = xmin;
-        xbounds[xx * 2 + 1] = xmax;
-    }
-
     imOut = ImagingNew(imIn->mode, xsize, imIn->ysize);
     if ( ! imOut) {
         free(kk);

From f3f03bcd675123eb9b93943320ba891a79f64be6 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 18:11:20 +0300
Subject: [PATCH 11/21] unroll 32bpc version

---
 libImaging/Resample.c | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index def135524e4..8513f4d467e 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -318,10 +318,10 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
 
     ImagingSectionEnter(&cookie);
     /* horizontal stretch */
-    for (yy = 0; yy < imOut->ysize; yy++) {
-        switch(imIn->type) {
-            case IMAGING_TYPE_INT32:
-                /* 32-bit integer */
+    
+    switch(imIn->type) {
+        case IMAGING_TYPE_INT32:
+            for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
                     xmax = xbounds[xx * 2 + 1];
@@ -331,9 +331,11 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                         ss += IMAGING_PIXEL_I(imIn, x, yy) * k[x - xmin];
                     IMAGING_PIXEL_I(imOut, xx, yy) = lround(ss);
                 }
-                break;
-            case IMAGING_TYPE_FLOAT32:
-                /* 32-bit float */
+            }
+            break;
+            
+        case IMAGING_TYPE_FLOAT32:
+            for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
                     xmax = xbounds[xx * 2 + 1];
@@ -343,8 +345,8 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                         ss += IMAGING_PIXEL_F(imIn, x, yy) * k[x - xmin];
                     IMAGING_PIXEL_F(imOut, xx, yy) = ss;
                 }
-                break;
-        }
+            }
+            break;
     }
     ImagingSectionLeave(&cookie);
     free(kk);

From 126c3151b5fb6534a04792762366d92b7a1efa29 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 18:25:42 +0300
Subject: [PATCH 12/21] use ImagingPrecompute in ImagingResampleHorizontal_8bpc

---
 libImaging/Resample.c | 78 +++++++++----------------------------------
 1 file changed, 16 insertions(+), 62 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 8513f4d467e..3baa9eedf93 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -143,78 +143,35 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
 {
     ImagingSectionCookie cookie;
     Imaging imOut;
-    float support, scale, filterscale;
-    float center, ww, ss;
     int ss0, ss1, ss2, ss3;
     int xx, yy, x, kmax, xmin, xmax;
     int *xbounds;
     int *k, *kk;
-    float *kw;
+    double *prekk;
 
-    /* prepare for horizontal stretch */
-    filterscale = scale = (float) imIn->xsize / xsize;
-    if (filterscale < 1.0) {
-        filterscale = 1.0;
-    }
-
-    /* determine support size (length of resampling filter) */
-    support = filterp->support * filterscale;
-
-    /* maximum number of coofs */
-    kmax = (int) ceil(support) * 2 + 1;
 
-    // check for overflow
-    if (xsize > SIZE_MAX / (kmax * sizeof(int)))
-        return (Imaging) ImagingError_MemoryError();
-
-    // sizeof(int) should be greater than 0 as well
-    if (xsize > SIZE_MAX / (2 * sizeof(int)))
-        return (Imaging) ImagingError_MemoryError();
-
-    /* coefficient buffer */
-    kk = malloc(xsize * kmax * sizeof(int));
-    if ( ! kk)
-        return (Imaging) ImagingError_MemoryError();
-
-    /* intermediate not normalized buffer for coefficients */ 
-    kw = malloc(kmax * sizeof(float));
-    if ( ! kw) {
-        free(kk);
+    kmax = ImagingPrecompute(imIn->xsize, xsize, filterp, &xbounds, &prekk);
+    if ( ! kmax) {
         return (Imaging) ImagingError_MemoryError();
     }
-
-    xbounds = malloc(xsize * 2 * sizeof(int));
-    if ( ! xbounds) {
-        free(kk);
-        free(kw);
+    
+    kk = malloc(xsize * kmax * sizeof(int));
+    if ( ! kk) {
+        free(xbounds);
+        free(prekk);
         return (Imaging) ImagingError_MemoryError();
     }
 
     for (xx = 0; xx < xsize; xx++) {
-        center = (xx + 0.5) * scale;
-        ww = 0.0;
-        ss = 1.0 / filterscale;
-        xmin = (int) floor(center - support);
-        if (xmin < 0)
-            xmin = 0;
-        xmax = (int) ceil(center + support);
-        if (xmax > imIn->xsize)
-            xmax = imIn->xsize;
-        for (x = xmin; x < xmax; x++) {
-            float w = filterp->filter((x - center + 0.5) * ss);
-            kw[x - xmin] = w;
-            ww += w;
-        }
+        xmin = xbounds[xx * 2 + 0];
+        xmax = xbounds[xx * 2 + 1];
         k = &kk[xx * kmax];
         for (x = 0; x < xmax - xmin; x++) {
-            if (ww != 0.0)
-                k[x] = (int) floor(0.5 + kw[x] / ww * (1 << PRECISION_BITS));
+            k[x] = (int) floor(0.5 + prekk[xx * kmax + x] * (1 << PRECISION_BITS));
         }
-        xbounds[xx * 2 + 0] = xmin;
-        xbounds[xx * 2 + 1] = xmax;
     }
 
-    free(kw);
+    free(prekk);
 
     imOut = ImagingNew(imIn->mode, xsize, imIn->ysize);
     if ( ! imOut) {
@@ -224,10 +181,8 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
     }
 
     ImagingSectionEnter(&cookie);
-    /* horizontal stretch */
     for (yy = 0; yy < imOut->ysize; yy++) {
         if (imIn->image8) {
-            /* 8-bit grayscale */
             for (xx = 0; xx < xsize; xx++) {
                 xmin = xbounds[xx * 2 + 0];
                 xmax = xbounds[xx * 2 + 1];
@@ -238,7 +193,6 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                 imOut->image8[yy][xx] = clip8(ss0);
             }
         } else if (imIn->type == IMAGING_TYPE_UINT8) {
-            /* n-bit grayscale */
             if (imIn->bands == 2) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
@@ -287,6 +241,7 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
             }
         }
     }
+
     ImagingSectionLeave(&cookie);
     free(kk);
     free(xbounds);
@@ -316,9 +271,7 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
         return NULL;
     }
 
-    ImagingSectionEnter(&cookie);
-    /* horizontal stretch */
-    
+    ImagingSectionEnter(&cookie);    
     switch(imIn->type) {
         case IMAGING_TYPE_INT32:
             for (yy = 0; yy < imOut->ysize; yy++) {
@@ -333,7 +286,7 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                 }
             }
             break;
-            
+
         case IMAGING_TYPE_FLOAT32:
             for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
@@ -348,6 +301,7 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
             }
             break;
     }
+
     ImagingSectionLeave(&cookie);
     free(kk);
     free(xbounds);

From d48324bd213685869bfbd52fe6bc95a7b1c7bc3a Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 5 May 2016 18:27:33 +0300
Subject: [PATCH 13/21] unroll 8bpc version

---
 libImaging/Resample.c | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 3baa9eedf93..30cbea115b4 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -181,8 +181,8 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
     }
 
     ImagingSectionEnter(&cookie);
-    for (yy = 0; yy < imOut->ysize; yy++) {
-        if (imIn->image8) {
+    if (imIn->image8) {
+        for (yy = 0; yy < imOut->ysize; yy++) {
             for (xx = 0; xx < xsize; xx++) {
                 xmin = xbounds[xx * 2 + 0];
                 xmax = xbounds[xx * 2 + 1];
@@ -192,8 +192,10 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     ss0 += ((UINT8) imIn->image8[yy][x]) * k[x - xmin];
                 imOut->image8[yy][xx] = clip8(ss0);
             }
-        } else if (imIn->type == IMAGING_TYPE_UINT8) {
-            if (imIn->bands == 2) {
+        }
+    } else if (imIn->type == IMAGING_TYPE_UINT8) {
+        if (imIn->bands == 2) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
                     xmax = xbounds[xx * 2 + 1];
@@ -206,7 +208,9 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     imOut->image[yy][xx*4 + 0] = clip8(ss0);
                     imOut->image[yy][xx*4 + 3] = clip8(ss1);
                 }
-            } else if (imIn->bands == 3) {
+            }
+        } else if (imIn->bands == 3) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
                     xmax = xbounds[xx * 2 + 1];
@@ -221,7 +225,9 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     imOut->image[yy][xx*4 + 1] = clip8(ss1);
                     imOut->image[yy][xx*4 + 2] = clip8(ss2);
                 }
-            } else {
+            }
+        } else {
+            for (yy = 0; yy < imOut->ysize; yy++) {
                 for (xx = 0; xx < xsize; xx++) {
                     xmin = xbounds[xx * 2 + 0];
                     xmax = xbounds[xx * 2 + 1];

From babaaf9bbe42b3803a58fb994edbb34731ae9399 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Fri, 6 May 2016 00:28:13 +0300
Subject: [PATCH 14/21] raise on special cases even if imIn->image8

---
 libImaging/Resample.c | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 30cbea115b4..915da6a7697 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -326,7 +326,9 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
     if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
         return (Imaging) ImagingError_ModeError();
 
-    if (imIn->image8) {
+    if (imIn->type == IMAGING_TYPE_SPECIAL) {
+        return (Imaging) ImagingError_ModeError();
+    } else if (imIn->image8) {
         ResampleHorizontal = ImagingResampleHorizontal_8bpc;
     } else {
         switch(imIn->type) {

From f44a7f8b11feb47175bad87ed498d9d8fad12202 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 9 May 2016 13:37:50 +0200
Subject: [PATCH 15/21] implement round fix windows build

---
 Tests/test_image_resample.py | 43 ++++++++++++++++++------------------
 libImaging/Resample.c        |  6 ++++-
 2 files changed, 26 insertions(+), 23 deletions(-)

diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 7d08f7ec4a6..efd2b4c2989 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -156,35 +156,34 @@ def make_case(self, mode, fill):
         im = Image.new(mode, (512, 9), fill)
         return (im.resize((9, 512), Image.LANCZOS), im.load()[0, 0])
 
-    def run_cases(self, cases):
-        for channel, color in cases:
-            px = channel.load()
-            for x in range(channel.size[0]):
-                for y in range(channel.size[1]):
-                    if px[x, y] != color:
-                        message = "{} != {} for pixel {}".format(
-                            px[x, y], color, (x, y))
-                        self.assertEqual(px[x, y], color, message)
+    def run_case(self, case):
+        channel, color = case
+        px = channel.load()
+        for x in range(channel.size[0]):
+            for y in range(channel.size[1]):
+                if px[x, y] != color:
+                    message = "{} != {} for pixel {}".format(
+                        px[x, y], color, (x, y))
+                    self.assertEqual(px[x, y], color, message)
 
     def test_8u(self):
         im, color = self.make_case('RGB', (0, 64, 255))
-        self.run_cases(zip(im.split(), color))
+        r, g, b = im.split()
+        self.run_case((r, color[0]))
+        self.run_case((g, color[1]))
+        self.run_case((b, color[2]))
 
     def test_32i(self):
-        self.run_cases([
-            self.make_case('I', 12),
-            self.make_case('I', 0x7fffffff),
-            self.make_case('I', -12),
-            self.make_case('I', -1 << 31),
-        ])
+        self.run_case(self.make_case('I', 12))
+        self.run_case(self.make_case('I', 0x7fffffff))
+        self.run_case(self.make_case('I', -12))
+        self.run_case(self.make_case('I', -1 << 31))
 
     def test_32f(self):
-        self.run_cases([
-            self.make_case('F', 1),
-            self.make_case('F', 3.40282306074e+38),
-            self.make_case('F', 1.175494e-38),
-            self.make_case('F', 1.192093e-07),
-        ])
+        self.run_case(self.make_case('F', 1))
+        self.run_case(self.make_case('F', 3.40282306074e+38))
+        self.run_case(self.make_case('F', 1.175494e-38))
+        self.run_case(self.make_case('F', 1.192093e-07))
 
 
 if __name__ == '__main__':
diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 915da6a7697..4e446af22a1 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -2,6 +2,10 @@
 
 #include <math.h>
 
+
+#define ROUND_UP(f) ((int) ((f) >= 0.0 ? (f) + 0.5F : (f) - 0.5F))
+
+
 struct filter {
     float (*filter)(float x);
     float support;
@@ -288,7 +292,7 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                     ss = 0.0;
                     for (x = xmin; x < xmax; x++)
                         ss += IMAGING_PIXEL_I(imIn, x, yy) * k[x - xmin];
-                    IMAGING_PIXEL_I(imOut, xx, yy) = lround(ss);
+                    IMAGING_PIXEL_I(imOut, xx, yy) = ROUND_UP(ss);
                 }
             }
             break;

From d3749ccc7cbf722bc0ecc2a94498565d2869b4c1 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 9 May 2016 13:48:14 +0200
Subject: [PATCH 16/21] possible increase precision and avoid compiler warnings

---
 libImaging/Resample.c | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 4e446af22a1..723b9729c14 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -7,11 +7,11 @@
 
 
 struct filter {
-    float (*filter)(float x);
-    float support;
+    double (*filter)(double x);
+    double support;
 };
 
-static inline float sinc_filter(float x)
+static inline double sinc_filter(double x)
 {
     if (x == 0.0)
         return 1.0;
@@ -19,7 +19,7 @@ static inline float sinc_filter(float x)
     return sin(x) / x;
 }
 
-static inline float lanczos_filter(float x)
+static inline double lanczos_filter(double x)
 {
     /* truncated sinc */
     if (-3.0 <= x && x < 3.0)
@@ -27,7 +27,7 @@ static inline float lanczos_filter(float x)
     return 0.0;
 }
 
-static inline float bilinear_filter(float x)
+static inline double bilinear_filter(double x)
 {
     if (x < 0.0)
         x = -x;
@@ -36,7 +36,7 @@ static inline float bilinear_filter(float x)
     return 0.0;
 }
 
-static inline float bicubic_filter(float x)
+static inline double bicubic_filter(double x)
 {
     /* https://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm */
 #define a -0.5
@@ -83,7 +83,7 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
     double *kk, *k;
 
     /* prepare for horizontal stretch */
-    filterscale = scale = (float) inSize / outSize;
+    filterscale = scale = (double) inSize / outSize;
     if (filterscale < 1.0) {
         filterscale = 1.0;
     }

From 7722d00daff339900ad106cbef074f153064a022 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 9 May 2016 14:33:14 +0200
Subject: [PATCH 17/21] simpler coefficients calculation in
 ImagingResampleHorizontal_8bpc

---
 libImaging/Resample.c | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 723b9729c14..dc1ba781c35 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -103,11 +103,11 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
         return 0;
 
     /* coefficient buffer */
-    kk = malloc(outSize * kmax * sizeof(double));
+    kk = calloc(outSize * kmax, sizeof(double));
     if ( ! kk)
         return 0;
 
-    xbounds = malloc(outSize * 2 * sizeof(int));
+    xbounds = calloc(outSize * 2, sizeof(int));
     if ( ! xbounds) {
         free(kk);
         return 0;
@@ -159,20 +159,15 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
         return (Imaging) ImagingError_MemoryError();
     }
     
-    kk = malloc(xsize * kmax * sizeof(int));
+    kk = calloc(xsize * kmax, sizeof(int));
     if ( ! kk) {
         free(xbounds);
         free(prekk);
         return (Imaging) ImagingError_MemoryError();
     }
 
-    for (xx = 0; xx < xsize; xx++) {
-        xmin = xbounds[xx * 2 + 0];
-        xmax = xbounds[xx * 2 + 1];
-        k = &kk[xx * kmax];
-        for (x = 0; x < xmax - xmin; x++) {
-            k[x] = (int) floor(0.5 + prekk[xx * kmax + x] * (1 << PRECISION_BITS));
-        }
+    for (x = 0; x < xsize * kmax; x++) {
+        kk[x] = (int) (0.5 + prekk[x] * (1 << PRECISION_BITS));
     }
 
     free(prekk);

From eec17c0c7aa3070bb63aa637c576e5b084818809 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 9 May 2016 14:39:08 +0200
Subject: [PATCH 18/21] change xmax meaning to avoid substructions

---
 libImaging/Resample.c | 45 ++++++++++++++++++++++---------------------
 1 file changed, 23 insertions(+), 22 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index dc1ba781c35..c248d31591e 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -123,13 +123,14 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
         xmax = (int) ceil(center + support);
         if (xmax > inSize)
             xmax = inSize;
+        xmax -= xmin;
         k = &kk[xx * kmax];
-        for (x = xmin; x < xmax; x++) {
-            double w = filterp->filter((x - center + 0.5) * ss);
-            k[x - xmin] = w;
+        for (x = 0; x < xmax; x++) {
+            double w = filterp->filter((x + xmin - center + 0.5) * ss);
+            k[x] = w;
             ww += w;
         }
-        for (x = 0; x < xmax - xmin; x++) {
+        for (x = 0; x < xmax; x++) {
             if (ww != 0.0)
                 k[x] /= ww;
         }
@@ -187,8 +188,8 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                 xmax = xbounds[xx * 2 + 1];
                 k = &kk[xx * kmax];
                 ss0 = 1 << (PRECISION_BITS -1);
-                for (x = xmin; x < xmax; x++)
-                    ss0 += ((UINT8) imIn->image8[yy][x]) * k[x - xmin];
+                for (x = 0; x < xmax; x++)
+                    ss0 += ((UINT8) imIn->image8[yy][x + xmin]) * k[x];
                 imOut->image8[yy][xx] = clip8(ss0);
             }
         }
@@ -200,9 +201,9 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     xmax = xbounds[xx * 2 + 1];
                     k = &kk[xx * kmax];
                     ss0 = ss1 = 1 << (PRECISION_BITS -1);
-                    for (x = xmin; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                        ss1 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
+                    for (x = 0; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss1 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 3]) * k[x];
                     }
                     imOut->image[yy][xx*4 + 0] = clip8(ss0);
                     imOut->image[yy][xx*4 + 3] = clip8(ss1);
@@ -215,10 +216,10 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     xmax = xbounds[xx * 2 + 1];
                     k = &kk[xx * kmax];
                     ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
-                    for (x = xmin; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                        ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                        ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
+                    for (x = 0; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss1 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 1]) * k[x];
+                        ss2 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 2]) * k[x];
                     }
                     imOut->image[yy][xx*4 + 0] = clip8(ss0);
                     imOut->image[yy][xx*4 + 1] = clip8(ss1);
@@ -232,11 +233,11 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
                     xmax = xbounds[xx * 2 + 1];
                     k = &kk[xx * kmax];
                     ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
-                    for (x = xmin; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
-                        ss1 += ((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
-                        ss2 += ((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
-                        ss3 += ((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
+                    for (x = 0; x < xmax; x++) {
+                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss1 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 1]) * k[x];
+                        ss2 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 2]) * k[x];
+                        ss3 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 3]) * k[x];
                     }
                     imOut->image[yy][xx*4 + 0] = clip8(ss0);
                     imOut->image[yy][xx*4 + 1] = clip8(ss1);
@@ -285,8 +286,8 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                     xmax = xbounds[xx * 2 + 1];
                     k = &kk[xx * kmax];
                     ss = 0.0;
-                    for (x = xmin; x < xmax; x++)
-                        ss += IMAGING_PIXEL_I(imIn, x, yy) * k[x - xmin];
+                    for (x = 0; x < xmax; x++)
+                        ss += IMAGING_PIXEL_I(imIn, x + xmin, yy) * k[x];
                     IMAGING_PIXEL_I(imOut, xx, yy) = ROUND_UP(ss);
                 }
             }
@@ -299,8 +300,8 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
                     xmax = xbounds[xx * 2 + 1];
                     k = &kk[xx * kmax];
                     ss = 0.0;
-                    for (x = xmin; x < xmax; x++)
-                        ss += IMAGING_PIXEL_F(imIn, x, yy) * k[x - xmin];
+                    for (x = 0; x < xmax; x++)
+                        ss += IMAGING_PIXEL_F(imIn, x + xmin, yy) * k[x];
                     IMAGING_PIXEL_F(imOut, xx, yy) = ss;
                 }
             }

From 555d2ace49079e95ffa6f95860caa4f952059497 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 9 May 2016 19:29:05 +0200
Subject: [PATCH 19/21] add test for L

---
 Tests/test_image_resample.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index efd2b4c2989..90dedae0400 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -172,6 +172,7 @@ def test_8u(self):
         self.run_case((r, color[0]))
         self.run_case((g, color[1]))
         self.run_case((b, color[2]))
+        self.run_case(self.make_case('L', 12))
 
     def test_32i(self):
         self.run_case(self.make_case('I', 12))

From 1c3def1d5a1fb48aa243f955786ce28f3f7bed62 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Mon, 16 May 2016 15:45:20 +0300
Subject: [PATCH 20/21] return malloc

---
 libImaging/Resample.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index c248d31591e..6cb622f230f 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -103,11 +103,11 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
         return 0;
 
     /* coefficient buffer */
-    kk = calloc(outSize * kmax, sizeof(double));
+    kk = malloc(outSize * kmax * sizeof(double));
     if ( ! kk)
         return 0;
 
-    xbounds = calloc(outSize * 2, sizeof(int));
+    xbounds = malloc(outSize * 2 * sizeof(int));
     if ( ! xbounds) {
         free(kk);
         return 0;
@@ -160,7 +160,7 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
         return (Imaging) ImagingError_MemoryError();
     }
     
-    kk = calloc(xsize * kmax, sizeof(int));
+    kk = malloc(xsize * kmax * sizeof(int));
     if ( ! kk) {
         free(xbounds);
         free(prekk);

From 5ffd9e53bc26c4ca1500c99447c2a6daeee60aa7 Mon Sep 17 00:00:00 2001
From: homm <homm86@gmail.com>
Date: Thu, 26 May 2016 02:28:35 +0300
Subject: [PATCH 21/21] use calloc and INT_MAX

---
 Tests/test_image_resample.py |  1 -
 libImaging/Resample.c        | 10 +++++-----
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 90dedae0400..0b78217ce56 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -151,7 +151,6 @@ def test_enlarge_lanczos(self):
 
 
 class CoreResampleConsistencyTest(PillowTestCase):
-
     def make_case(self, mode, fill):
         im = Image.new(mode, (512, 9), fill)
         return (im.resize((9, 512), Image.LANCZOS), im.load()[0, 0])
diff --git a/libImaging/Resample.c b/libImaging/Resample.c
index 6cb622f230f..8ad1c9d97a0 100644
--- a/libImaging/Resample.c
+++ b/libImaging/Resample.c
@@ -95,19 +95,19 @@ ImagingPrecompute(int inSize, int outSize, struct filter *filterp,
     kmax = (int) ceil(support) * 2 + 1;
 
     // check for overflow
-    if (outSize > SIZE_MAX / (kmax * sizeof(double)))
+    if (outSize > INT_MAX / (kmax * sizeof(double)))
         return 0;
 
     // sizeof(double) should be greater than 0 as well
-    if (outSize > SIZE_MAX / (2 * sizeof(double)))
+    if (outSize > INT_MAX / (2 * sizeof(double)))
         return 0;
 
     /* coefficient buffer */
-    kk = malloc(outSize * kmax * sizeof(double));
+    kk = calloc(outSize * kmax, sizeof(double));
     if ( ! kk)
         return 0;
 
-    xbounds = malloc(outSize * 2 * sizeof(int));
+    xbounds = calloc(outSize * 2, sizeof(int));
     if ( ! xbounds) {
         free(kk);
         return 0;
@@ -160,7 +160,7 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
         return (Imaging) ImagingError_MemoryError();
     }
     
-    kk = malloc(xsize * kmax * sizeof(int));
+    kk = calloc(xsize * kmax, sizeof(int));
     if ( ! kk) {
         free(xbounds);
         free(prekk);