cpu/complex.h

#ifndef PSCF_CPU_COMPLEX_H
#define PSCF_CPU_COMPLEX_H
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <pscf/math/arithmetic.h>  // built-in and standard lib types
 
#include <fftw3.h>
#include <complex>
#include <iostream>
 
namespace Pscf{

 
   // Real and imaginary components of fftw_complex numbers

   inline 
   double real(fftw_complex const & a)
   {  return a[0]; }

   inline 
   double imag(fftw_complex const & a)
   {  return a[1]; }
 
   // Absolute magnitude

   inline 
   double abs(fftw_complex const & a)
   {  return sqrt(a[0] * a[0] + a[1] * a[1]); }

   inline 
   double absSq(fftw_complex const & a)
   {  return (a[0] * a[0] + a[1] * a[1]); }
 
   // Complex Conjugation

   inline
   void conj(fftw_complex& z, fftw_complex const & a)
   {
      z[0] = a[0];
      z[1] = -a[1];
   }

   inline
   void conj(fftw_complex& a)
   {
      a[0] = a[0];
      a[1] = -a[1];
   }
 
   // Assignment

   inline
   void assign(fftw_complex& z, double const & a, double const & b)
   {
      z[0] = a;
      z[1] = b;
   }

   inline
   void assign(fftw_complex& z, double const & a)
   {
      z[0] = a;
      z[1] = 0.0;
   }

   inline
   void assign(fftw_complex& z, fftw_complex const & a)
   {
      z[0] = a[0];
      z[1] = a[1];
   }

   inline
   void assign(fftw_complex & z, std::complex<double> const& a)
   {
      z[0] = a.real();
      z[1] = a.imag();
   }

   inline
   void assign(std::complex<double> & z, fftw_complex const& a)
   {  z = std::complex<double>(a[0], a[1]); }
 
   // Addition

   inline
   void add(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] + b[0];
      z[1] = a[1] + b[1];
   }

   inline
   void add(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0] + b;
      z[1] = a[1];
   }

   inline
   void addEq(fftw_complex& a, fftw_complex const& b)
   {
      a[0] += b[0];
      a[1] += b[1];
   }

   inline
   void addEq(fftw_complex& a, double const& b)
   {
      a[0] += b;
   }
 
   // Subtraction

   inline
   void sub(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] - b[0];
      z[1] = a[1] - b[1];
   }

   inline
   void sub(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0] - b;
      z[1] = a[1];
   }

   inline
   void subEq(fftw_complex & a, fftw_complex const& b)
   {
      a[0] -= b[0];
      a[1] -= b[1];
   }

   inline
   void subEq(fftw_complex & a, double const& b)
   {
      a[0] -= b;
   }

   inline 
   double absSqDiff(fftw_complex const& a, fftw_complex const& b)
   {
      fftw_complex z;
      sub(z, a, b);
      return absSq(z);
   }
 
   // Multiplication

   inline
   void mul(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] * b[0] - a[1] * b[1];
      z[1] = a[1] * b[0] + a[0] * b[1];
   }

   inline
   void mul(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0]*b;
      z[1] = a[1]*b;
   }

   inline
   void mulEq(fftw_complex & a, fftw_complex const& b)
   {
      double a0;
      a0   = a[0] * b[0] - a[1] * b[1];
      a[1] = a[1] * b[0] + a[0] * b[1];
      a[0] = a0;
   }

   inline
   void mulEq(fftw_complex & a, double const& b)
   {
      a[0] *= b;
      a[1] *= b;
   }

   inline
   void square(fftw_complex& z, fftw_complex const& a)
   {
      z[0] = a[0] * a[0] - a[1] * a[1];
      z[1] = 2.0 * a[1] * a[0];
   }
 
   // Division

   inline
   void div(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      double bSq = b[0] * b[0] + b[1] * b[1];
      z[0] = (a[0] * b[0] + a[1] * b[1])/bSq;
      z[1] = (a[1] * b[0] - a[0] * b[1])/bSq;
   }

   inline
   void div(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0]/b;
      z[1] = a[1]/b;
   }

   inline
   void divEq(fftw_complex & a, fftw_complex const & b)
   {
      double bSq = b[0] * b[0] + b[1] * b[1];
      double a0 = (a[0] * b[0] + a[1] * b[1])/bSq;
      a[1] = (a[1] * b[0] - a[0] * b[1])/bSq;
      a[0] = a0;
   }

   inline
   void divEq(fftw_complex & a, double const & b)
   {
      a[0] /= b;
      a[1] /= b;
   }
 
   // Inversion

   inline
   void inverse(fftw_complex& z, fftw_complex const & a)
   {
      double aSq = a[0] * a[0] + a[1] * a[1];
      z[0] =   a[0]/aSq;
      z[1] = - a[1]/aSq;
   }
 
   // Exponentiation and logarithm

   inline
   void assignExp(fftw_complex & z, fftw_complex const & a)
   {
      std::complex<double> arg = std::complex<double>(a[0], a[1]); 
      std::complex<double> result = std::exp(arg);
      z[0] = result.real();
      z[1] = result.imag();
   }

   inline
   void assignLog(fftw_complex & z, fftw_complex const & a)
   {  
      std::complex<double> arg = std::complex<double>(a[0], a[1]); 
      std::complex<double> result = std::log(arg);
      z[0] = result.real();
      z[1] = result.imag();
   }
 
   // Stream IO operators

   std::istream& operator >> (std::istream& is, fftw_complex & z);

   std::ostream& operator << (std::ostream& os, fftw_complex const & z);
 
} // namespace Pscf
#endif