static PyObject *
complex_richcompare(PyObject *v, PyObject *w, int op)
{
    PyObject *res;
    Py_complex i;
    int equal;

    if (op != Py_EQ && op != Py_NE) {
        goto Unimplemented;
    }

    assert(PyComplex_Check(v));
    TO_COMPLEX(v, i);

    if (PyLong_Check(w)) {
        /* Check for 0.0 imaginary part first to avoid the rich
         * comparison when possible.
         */
        if (i.imag == 0.0) {
            PyObject *j, *sub_res;
            j = PyFloat_FromDouble(i.real);
            if (j == NULL)
                return NULL;

            sub_res = PyObject_RichCompare(j, w, op);
            Py_DECREF(j);
            return sub_res;
        }
        else {
            equal = 0;
        }
    }
    else if (PyFloat_Check(w)) {
        equal = (i.real == PyFloat_AsDouble(w) && i.imag == 0.0);
    }
    else if (PyComplex_Check(w)) {
        Py_complex j;

        TO_COMPLEX(w, j);
        equal = (i.real == j.real && i.imag == j.imag);
    }
    else {
        goto Unimplemented;
    }

    if (equal == (op == Py_EQ))
         res = Py_True;
    else
         res = Py_False;

    Py_INCREF(res);
    return res;

Unimplemented:
    Py_RETURN_NOTIMPLEMENTED;
}

Py_complex
PyComplex_AsCComplex(PyObject *op)
{
    Py_complex cv;
    PyObject *newop = NULL;
    static PyObject *complex_str = NULL;

    assert(op);
    /* If op is already of type PyComplex_Type, return its value */
    if (PyComplex_Check(op)) {
        return ((PyComplexObject *)op)->cval;
    }
    /* If not, use op's __complex__  method, if it exists */

    /* return -1 on failure */
    cv.real = -1.;
    cv.imag = 0.;

    if (complex_str == NULL) {
        if (!(complex_str = PyUnicode_FromString("__complex__")))
            return cv;
    }

    {
        PyObject *complexfunc;
        complexfunc = _PyType_Lookup(op->ob_type, complex_str);
        /* complexfunc is a borrowed reference */
        if (complexfunc) {
            newop = PyObject_CallFunctionObjArgs(complexfunc, op, NULL);
            if (!newop)
                return cv;
        }
    }

    if (newop) {
        if (!PyComplex_Check(newop)) {
            PyErr_SetString(PyExc_TypeError,
                "__complex__ should return a complex object");
            Py_DECREF(newop);
            return cv;
        }
        cv = ((PyComplexObject *)newop)->cval;
        Py_DECREF(newop);
        return cv;
    }
    /* If neither of the above works, interpret op as a float giving the
       real part of the result, and fill in the imaginary part as 0. */
    else {
        /* PyFloat_AsDouble will return -1 on failure */
        cv.real = PyFloat_AsDouble(op);
        return cv;
    }
}

static int
complex_coerce(PyObject **pv, PyObject **pw)
{
	Py_complex cval;
	cval.imag = 0.;
	if (PyInt_Check(*pw)) {
		cval.real = (double)PyInt_AsLong(*pw);
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	else if (PyLong_Check(*pw)) {
		cval.real = PyLong_AsDouble(*pw);
		if (cval.real == -1.0 && PyErr_Occurred())
			return -1;
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	else if (PyFloat_Check(*pw)) {
		cval.real = PyFloat_AsDouble(*pw);
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	else if (PyComplex_Check(*pw)) {
		Py_INCREF(*pv);
		Py_INCREF(*pw);
		return 0;
	}
	return 1; /* Can't do it */
}

static MPC_Object *
GMPy_MPC_From_PyComplex(PyObject *obj, mpfr_prec_t rprec, mpfr_prec_t iprec,
                        CTXT_Object *context)
{
    MPC_Object *result;

    assert(PyComplex_Check(obj));

    CHECK_CONTEXT(context);

    if (rprec == 0)
        rprec = GET_REAL_PREC(context);
    else if (rprec == 1)
        rprec = DBL_MANT_DIG;

    if (iprec == 0)
        iprec = GET_IMAG_PREC(context);
    else if (iprec == 1)
        rprec = DBL_MANT_DIG;

    if ((result = GMPy_MPC_New(rprec, iprec, context))) {
        result->rc = mpc_set_d_d(result->c, PyComplex_RealAsDouble(obj),
                                 PyComplex_ImagAsDouble(obj), GET_MPC_ROUND(context));
        if (rprec != 1 || iprec != 1) {
            GMPY_MPC_CHECK_RANGE(result, context);
        }
        GMPY_MPC_SUBNORMALIZE(result, context);
        GMPY_MPC_EXCEPTIONS(result, context);
    }
    return result;
}

int extractComplexList(PyObject *list_in, double2** p_array_out, int argnum) {

  //Checking input type
  if (!PyList_Check(list_in)) {
     char errmsg[37];
     sprintf(errmsg,"Argument %d must be a list of complex\n", argnum);fflush(0);
     printf("Argument %d must be a list of complex\n", argnum);fflush(0);
     fprintf(stderr, "%s", errmsg);
     PyErr_SetString(PyExc_TypeError, errmsg); 
     return -1;
  }
  if (!PyComplex_Check(PyList_GetItem(list_in,0))) {
     char errmsg[37];
     sprintf(errmsg,"Argument %d must be a list of complex\n", argnum);fflush(0);
     printf("Argument %d must be a list of complex\n", argnum);fflush(0);
     fprintf(stderr, "%s", errmsg);
     PyErr_SetString(PyExc_TypeError, errmsg); 
     return -1;
  }

  int list_size = PyList_Size(list_in);
  PyObject* iter = PyObject_GetIter(list_in);
  *p_array_out = (double2*)malloc(sizeof(double2)*list_size);
  double2* array_out = *p_array_out; 
  int q;
  PyObject* item;
  for (q=0;q<PyList_Size(list_in);q++) {
     item = PyIter_Next(iter);
     array_out[q].x = PyComplex_RealAsDouble(item);
     array_out[q].y = PyComplex_ImagAsDouble(item);
  }
  return 0;
}

Py_complex
PyComplex_AsCComplex(PyObject *op)
{
    Py_complex cv;
    PyObject *newop = NULL;

    assert(op);
    /* If op is already of type PyComplex_Type, return its value */
    if (PyComplex_Check(op)) {
        return ((PyComplexObject *)op)->cval;
    }
    /* If not, use op's __complex__  method, if it exists */

    /* return -1 on failure */
    cv.real = -1.;
    cv.imag = 0.;

    newop = try_complex_special_method(op);

    if (newop) {
        cv = ((PyComplexObject *)newop)->cval;
        Py_DECREF(newop);
        return cv;
    }
    else if (PyErr_Occurred()) {
        return cv;
    }
    /* If neither of the above works, interpret op as a float giving the
       real part of the result, and fill in the imaginary part as 0. */
    else {
        /* PyFloat_AsDouble will return -1 on failure */
        cv.real = PyFloat_AsDouble(op);
        return cv;
    }
}

 static std::complex<double> py2c(PyObject *ob) {
  if (PyComplex_Check(ob)) {
    auto r = PyComplex_AsCComplex(ob);
    return {r.real, r.imag};
  }
  return PyFloat_AsDouble(ob);
 }

static PyObject *
try_complex_special_method(PyObject *op)
{
    PyObject *f;
    _Py_IDENTIFIER(__complex__);

    f = _PyObject_LookupSpecial(op, &PyId___complex__);
    if (f) {
        PyObject *res = _PyObject_CallNoArg(f);
        Py_DECREF(f);
        if (!res || PyComplex_CheckExact(res)) {
            return res;
        }
        if (!PyComplex_Check(res)) {
            PyErr_Format(PyExc_TypeError,
                "__complex__ returned non-complex (type %.200s)",
                res->ob_type->tp_name);
            Py_DECREF(res);
            return NULL;
        }
        /* Issue #29894: warn if 'res' not of exact type complex. */
        if (PyErr_WarnFormat(PyExc_DeprecationWarning, 1,
                "__complex__ returned non-complex (type %.200s).  "
                "The ability to return an instance of a strict subclass of complex "
                "is deprecated, and may be removed in a future version of Python.",
                res->ob_type->tp_name)) {
            Py_DECREF(res);
            return NULL;
        }
        return res;
    }
    return NULL;
}

NPY_NO_EXPORT PyArray_Descr *
_array_find_python_scalar_type(PyObject *op)
{
    if (PyFloat_Check(op)) {
        return PyArray_DescrFromType(NPY_DOUBLE);
    }
    else if (PyComplex_Check(op)) {
        return PyArray_DescrFromType(NPY_CDOUBLE);
    }
    else if (PyInt_Check(op)) {
        /* bools are a subclass of int */
        if (PyBool_Check(op)) {
            return PyArray_DescrFromType(NPY_BOOL);
        }
        else {
            return  PyArray_DescrFromType(NPY_LONG);
        }
    }
    else if (PyLong_Check(op)) {
        /* if integer can fit into a longlong then return that*/
        if ((PyLong_AsLongLong(op) == -1) && PyErr_Occurred()) {
            PyErr_Clear();
            return PyArray_DescrFromType(NPY_OBJECT);
        }
        return PyArray_DescrFromType(NPY_LONGLONG);
    }
    return NULL;
}

static int pyobject_to_ligotimegps(PyObject *obj, LIGOTimeGPS *gps)
{
	if(pylal_LIGOTimeGPS_Check(obj)) {
		*gps = ((pylal_LIGOTimeGPS *) obj)->gps;
	} else if(PyInt_Check(obj)) {
		XLALGPSSet(gps, PyInt_AsLong(obj), 0);
	} else if(PyLong_Check(obj)) {
		XLALGPSSet(gps, PyLong_AsLongLong(obj), 0);
	} else if(PyFloat_Check(obj)) {
		XLALGPSSetREAL8(gps, PyFloat_AsDouble(obj));
	} else if(PyComplex_Check(obj)) {
		if(PyComplex_ImagAsDouble(obj) != 0.0) {
			XLALGPSSet(gps, 0, 0);
			PyErr_SetObject(PyExc_ValueError, obj);
			return 0;
		}
		XLALGPSSetREAL8(gps, PyComplex_RealAsDouble(obj));
	} else {
		PyObject *s_attr = PyObject_GetAttrString(obj, "seconds");
		PyObject *n_attr = PyObject_GetAttrString(obj, "nanoseconds");
		XLALGPSSet(gps, PyInt_AsLong(s_attr), PyInt_AsLong(n_attr));
		Py_XDECREF(s_attr);
		Py_XDECREF(n_attr);
		if(PyErr_Occurred()) {
			PyErr_SetObject(PyExc_TypeError, obj);
			return 0;
		}
	}
	return 1;
}

static mxArray *numeric2mx(PyObject *pSrc)
{
  mxArray *lDst = NULL;

  pyassert(PyArray_API, "Unable to perform this function without NumPy installed");
  if (PyArray_Check(pSrc)) {
    lDst = makeMxFromNumeric((const PyArrayObject *)pSrc);
  } else if (PySequence_Check(pSrc)) {
    lDst = makeMxFromSeq(pSrc);
  } else if (PyObject_HasAttrString(pSrc, "__array__")) {
    PyObject *arp;
    arp = PyObject_CallMethod(pSrc, "__array__", NULL);
    lDst = makeMxFromNumeric((const PyArrayObject *)arp);
    Py_DECREF(arp);             // FIXME check this is correct;
  }
    else if (PyInt_Check(pSrc) || PyLong_Check(pSrc) ||
             PyFloat_Check(pSrc) || PyComplex_Check(pSrc)){
    PyObject *t;
    t = Py_BuildValue("(O)", pSrc);
//     t = PyTuple_New(1);
//     PyTuple_SetItem(t, 0, pSrc);
    lDst = makeMxFromSeq(t);
    Py_DECREF(t); // FIXME check this
  } else {

  }
  return lDst;

 error_return:
  return NULL;
  }

static int int_from_pyobj(int* v,PyObject *obj,const char *errmess) {
  PyObject* tmp = NULL;
  if (PyInt_Check(obj)) {
    *v = (int)PyInt_AS_LONG(obj);
    return 1;
  }
  tmp = PyNumber_Int(obj);
  if (tmp) {
    *v = PyInt_AS_LONG(tmp);
    Py_DECREF(tmp);
    return 1;
  }
  if (PyComplex_Check(obj))
    tmp = PyObject_GetAttrString(obj,"real");
  else if (PyString_Check(obj) || PyUnicode_Check(obj))
    /*pass*/;
  else if (PySequence_Check(obj))
    tmp = PySequence_GetItem(obj,0);
  if (tmp) {
    PyErr_Clear();
    if (int_from_pyobj(v,tmp,errmess)) {Py_DECREF(tmp); return 1;}
    Py_DECREF(tmp);
  }
  {
    PyObject* err = PyErr_Occurred();
    if (err==NULL) err = uts_scsmfo_error;
    PyErr_SetString(err,errmess);
  }
  return 0;
}

inline ComplexVector* ComplexVector_from_python(PyObject* py) {
  PyObject* seq = PySequence_Fast(py, "Argument must be a sequence of complex numbers.");
  if (seq == NULL)
    return 0;
  int size = PySequence_Fast_GET_SIZE(seq);
  ComplexVector* cpp = new ComplexVector(size);
  try {
    for (int i = 0; i < size; ++i) {
      PyObject* value = PySequence_Fast_GET_ITEM(seq, i);
      if (!PyComplex_Check(value)) {
        delete cpp;
        Py_DECREF(seq);
        PyErr_SetString(PyExc_TypeError, "Argument must be a sequence of complex numbers.");
        return 0;
      }
      Py_complex temp = PyComplex_AsCComplex(value);
      (*cpp)[i] = ComplexPixel(temp.real, temp.imag);
    }
  } catch (std::exception e) {
    delete cpp;
    Py_DECREF(seq);
    PyErr_SetString(PyExc_RuntimeError, e.what());
    return 0;
  }
  Py_DECREF(seq);
  return cpp;
}

 BOOST_PYTHON_DECL void expect_complex(PyObject* p)
 {
     if (!PyComplex_Check(p))
     {
         PyErr_SetString(PyExc_TypeError, "expected a complex number");
         boost::python::throw_argument_error();
     }
 }

//---------------------------------------------------------------------------
template<> double FromPyObject<double>(PyObject *O)
{
	if(PyComplex_Check(O))
    if(PyComplex_ImagAsDouble(O) != 0)
      PyErr_SetString(PyExc_TypeError, "complex number has an imaginary part");
		else
      return PyComplex_RealAsDouble(O);
  return PyFloat_AsDouble(O);
}

double
PyComplex_ImagAsDouble(PyObject *op)
{
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval.imag;
	}
	else {
		return 0.0;
	}
}

double
PyComplex_RealAsDouble(PyObject *op)
{
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval.real;
	}
	else {
		return PyFloat_AsDouble(op);
	}
}

static PyObject *
complex_richcompare(PyObject *v, PyObject *w, int op)
{
	int c;
	Py_complex i, j;
	PyObject *res;

	c = PyNumber_CoerceEx(&v, &w);
	if (c < 0)
		return NULL;
	if (c > 0) {
		Py_INCREF(Py_NotImplemented);
		return Py_NotImplemented;
	}
	/* Make sure both arguments are complex. */
	if (!(PyComplex_Check(v) && PyComplex_Check(w))) {
		Py_DECREF(v);
		Py_DECREF(w);
		Py_INCREF(Py_NotImplemented);
		return Py_NotImplemented;
	}

	i = ((PyComplexObject *)v)->cval;
	j = ((PyComplexObject *)w)->cval;
	Py_DECREF(v);
	Py_DECREF(w);

	if (op != Py_EQ && op != Py_NE) {
		PyErr_SetString(PyExc_TypeError,
			"no ordering relation is defined for complex numbers");
		return NULL;
	}

	if ((i.real == j.real && i.imag == j.imag) == (op == Py_EQ))
		res = Py_True;
	else
		res = Py_False;

	Py_INCREF(res);
	return res;
}

Py_complex
PyComplex_AsCComplex(PyObject *op)
{
	Py_complex cv;
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval;
	}
	else {
		cv.real = PyFloat_AsDouble(op);
		cv.imag = 0.;
		return cv;
	}   
}

inline ComplexPixel pixel_from_python<ComplexPixel>::convert(PyObject* obj) {
  if (!PyComplex_Check(obj)) {
    if (!is_RGBPixelObject(obj)) {
      if (!PyFloat_Check(obj)) {
        if (!PyInt_Check(obj)) {
          throw std::runtime_error("Pixel value is not convertible to a ComplexPixel");
        }
        return ComplexPixel((double)PyInt_AsLong(obj), 0.0);
      }
      return ComplexPixel(PyFloat_AsDouble(obj), 0.0);
    }
    return ComplexPixel(((RGBPixelObject*)obj)->m_x->luminance(), 0.0);
  }
  Py_complex temp = PyComplex_AsCComplex(obj);
  return ComplexPixel(temp.real, temp.imag);
}