diff --git a/Tests/test_imaging_stretch.py b/Tests/test_imaging_stretch.py
new file mode 100644
index 000000000..d2fbe1c79
--- /dev/null
+++ b/Tests/test_imaging_stretch.py
@@ -0,0 +1,40 @@
+"""
+Tests for ImagingCore.stretch functionality.
+"""
+
+from helper import unittest, PillowTestCase
+
+from PIL import Image
+
+
+im = Image.open("Tests/images/hopper.ppm").copy()
+
+
+class TestImagingStretch(PillowTestCase):
+
+    def test_modes(self):
+        self.assertRaises(ValueError, im.convert("1").im.stretch,
+                          (15, 12), Image.ANTIALIAS)
+        self.assertRaises(ValueError, im.convert("P").im.stretch,
+                          (15, 12), Image.ANTIALIAS)
+        for mode in ["L", "I", "F", "RGB", "RGBA", "CMYK", "YCbCr"]:
+            s = im.convert(mode).im
+            r = s.stretch((15, 12), Image.ANTIALIAS)
+            self.assertEqual(r.mode, mode)
+            self.assertEqual(r.size, (15, 12))
+            self.assertEqual(r.bands, s.bands)
+
+    def test_reduce_filters(self):
+        # There is no Image.NEAREST because im.stretch implementation
+        # is not NEAREST for reduction. It should be removed
+        # or renamed to supersampling.
+        for f in [Image.BILINEAR, Image.BICUBIC, Image.ANTIALIAS]:
+            r = im.im.stretch((15, 12), f)
+            self.assertEqual(r.mode, "RGB")
+            self.assertEqual(r.size, (15, 12))
+
+    def test_enlarge_filters(self):
+        for f in [Image.BILINEAR, Image.BICUBIC, Image.ANTIALIAS]:
+            r = im.im.stretch((764, 414), f)
+            self.assertEqual(r.mode, "RGB")
+            self.assertEqual(r.size, (764, 414))
diff --git a/libImaging/Antialias.c b/libImaging/Antialias.c
index d413fbb6a..be49bc827 100644
--- a/libImaging/Antialias.c
+++ b/libImaging/Antialias.c
@@ -64,16 +64,14 @@ static struct filter BILINEAR = { bilinear_filter, 1.0 };
 
 static inline float bicubic_filter(float x)
 {
-    /* FIXME: double-check this algorithm */
-    /* FIXME: for best results, "a" should be -0.5 to -1.0, but we'll
-       set it to zero for now, to match the 1.1 magnifying filter */
-#define a 0.0
+    /* http://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm */
+#define a -0.5
     if (x < 0.0)
         x = -x;
     if (x < 1.0)
-        return (((a + 2.0) * x) - (a + 3.0)) * x*x + 1;
+        return ((a + 2.0) * x - (a + 3.0)) * x*x + 1;
     if (x < 2.0)
-        return (((a * x) - 5*a) * x + 8) * x - 4*a;
+        return (((x - 5) * x + 8) * x - 4) * a;
     return 0.0;
 #undef a
 }
@@ -95,7 +93,10 @@ ImagingStretch(Imaging imOut, Imaging imIn, int filter)
 
     /* check modes */
     if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
-	return (Imaging) ImagingError_ModeError();
+        return (Imaging) ImagingError_ModeError();
+
+    if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
+        return (Imaging) ImagingError_ModeError();
 
     /* check filter */
     switch (filter) {
@@ -124,14 +125,13 @@ ImagingStretch(Imaging imOut, Imaging imIn, int filter)
         /* prepare for vertical stretch */
         filterscale = scale = (float) imIn->ysize / imOut->ysize;
     } else
-	return (Imaging) ImagingError_Mismatch();
+        return (Imaging) ImagingError_Mismatch();
 
     /* determine support size (length of resampling filter) */
     support = filterp->support;
 
     if (filterscale < 1.0) {
         filterscale = 1.0;
-        support = 0.5;
     }
 
     support = support * filterscale;
@@ -154,7 +154,7 @@ ImagingStretch(Imaging imOut, Imaging imIn, int filter)
                 ymin = 0.0;
             ymax = ceil(center + support);
             if (ymax > (float) imIn->ysize)
-		ymax = (float) imIn->ysize;
+                ymax = (float) imIn->ysize;
             for (y = (int) ymin; y < (int) ymax; y++) {
                 float w = filterp->filter((y - center + 0.5) * ss) * ss;
                 k[y - (int) ymin] = w;
@@ -230,7 +230,7 @@ ImagingStretch(Imaging imOut, Imaging imIn, int filter)
                 xmin = 0.0;
             xmax = ceil(center + support);
             if (xmax > (float) imIn->xsize)
-		xmax = (float) imIn->xsize;
+                xmax = (float) imIn->xsize;
             for (x = (int) xmin; x < (int) xmax; x++) {
                 float w = filterp->filter((x - center + 0.5) * ss) * ss;
                 k[x - (int) xmin] = w;