{******************************************************************************}
{                                                                              }
{ Project JEDI Code Library (JCL)                                              }
{                                                                              }
{ The contents of this file are subject to the Mozilla Public License Version  }
{ 1.1 (the "License"); you may not use this file except in compliance with the }
{ License. You may obtain a copy of the License at http://www.mozilla.org/MPL/ }
{                                                                              }
{ Software distributed under the License is distributed on an "AS IS" basis,   }
{ WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License for }
{ the specific language governing rights and limitations under the License.    }
{                                                                              }
{ The Original Code is JclComplex.pas.                                         }
{                                                                              }
{ The Initial Developer of the Original Code is documented in the accompanying }
{ help file JCL.chm. Portions created by these individuals are Copyright (C)   }
{ of these individuals.                                                        }
{                                                                              }
{******************************************************************************}
{                                                                              }
{ Class for working with complex numbers.                                      }
{                                                                              }
{ Unit owner: Alexei Koudinov                                                  }
{ Last modified: January 30, 2001                                              }
{                                                                              }
{******************************************************************************}

unit JclComplex;

{$I jcl.inc}

{$WEAKPACKAGEUNIT ON}

interface

uses
  SysUtils,
  JclBase, JclMath, JclResources, JclStrings, JclSysUtils;

const
  TComplex_VERSION = 5.01;

type
  TComplexKind = (crRectangular, crPolar);

  TCoords = record
    X: Float;     // rectangular real
    Y: Float;     // rectangular imaginary
    R: Float;     // polar 1
    Theta: Float; // polar 2
  end;

  TRectCoord = record
    X: Float;
    Y: Float;
  end;

  TJclComplex = class (TObject)
  private   {z = x + yi}
    // -----------------------------
    FCoord: TCoords;
    FFracLen: Byte;
    function MiscalcSingle(const X: Float): Float;
    procedure MiscalcComplex; // eliminates miscalculation
    procedure FillCoords(const CompexType: TComplexKind);
    function GetRectangularString: string;
    function GetPolarString: string;
    procedure SetRectangularString(StrToParse: string);
    procedure SetPolarString(StrToParse: string);
    procedure SetFracLen(const X: Byte);
    function GetRadius: Float;
    function GetAngle: Float;
    function NormalizeAngle(Value: Float): Float;
  protected
    function Assign(const Coord: TCoords; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CoreAdd(const First, Second: TRectCoord): TRectCoord;
    function CoreDiv(const First, Second: TRectCoord): TRectCoord;
    function CoreMul(const First, Second: TRectCoord): TRectCoord;
    function CoreSub(const First, Second: TRectCoord): TRectCoord;
    function CoreLn (const LnValue: TRectCoord): TRectCoord;
    function CoreExp(const ExpValue: TRectCoord): TRectCoord;
    function CorePwr(First, Second, Polar: TRectCoord): TRectCoord;
    function CoreIntPwr(First: TRectCoord; const Polar: TRectCoord; const Pwr: Integer): TRectCoord;
    function CoreRealPwr(First: TRectCoord; const Polar: TRectCoord; const Pwr: Float): TRectCoord;
    function CoreRoot(First: TRectCoord; const Polar: TRectCoord; const K, N: Word): TRectCoord;
    function CoreCos(const Value: TRectCoord): TRectCoord;
    function CoreSin(const Value: TRectCoord): TRectCoord;
    function CoreTan(const Value: TRectCoord): TRectCoord;
    function CoreCot(const Value: TRectCoord): TRectCoord;
    function CoreSec(const Value: TRectCoord): TRectCoord;
    function CoreCsc(const Value: TRectCoord): TRectCoord;
    function CoreCosH(const Value: TRectCoord): TRectCoord;
    function CoreSinH(const Value: TRectCoord): TRectCoord;
    function CoreTanH(const Value: TRectCoord): TRectCoord;
    function CoreCotH(const Value: TRectCoord): TRectCoord;
    function CoreSecH(const Value: TRectCoord): TRectCoord;
    function CoreCscH(const Value: TRectCoord): TRectCoord;
    function CoreI0(const Value: TRectCoord): TRectCoord;
    function CoreJ0(const Value: TRectCoord): TRectCoord;
    function CoreApproxLnGamma(const Value: TRectCoord): TRectCoord;
    function CoreLnGamma(Value: TRectCoord): TRectCoord;
    function CoreGamma(const Value: TRectCoord): TRectCoord;
  public
    property FracLength: Byte read FFracLen write SetFracLen default 8;

    //----------- constructors
    constructor Create; overload;
    constructor Create(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}); overload;

    //----------- getting different parts of the number
    property RealPart: Float read FCoord.X;
    property ImaginaryPart: Float read FCoord.Y;
    property Radius: Float read GetRadius;
    property Angle: Float read GetAngle;

    //----------- format output
    property AsString: string read GetRectangularString write SetRectangularString;
    property AsPolarString: string read GetPolarString write SetPolarString;

    //----------- complex numbers assignment routines
    function Assign(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function AssignZero: TJclComplex;
    function AssignOne: TJclComplex;
    function Duplicate: TJclComplex;

    //----------- arithmetics -- modify the object itself
    function CAdd(const AddValue: TJclComplex): TJclComplex; overload;
    function CAdd(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CDiv(const DivValue: TJclComplex): TJclComplex; overload;
    function CDiv(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CMul(const MulValue: TJclComplex): TJclComplex; overload;
    function CMul(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CSub(const SubValue: TJclComplex): TJclComplex; overload;
    function CSub(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNeg: TJclComplex;
    function CConjugate: TJclComplex;

    //----------- arithmetics -- creates new resulting object
    function CNewAdd(const AddValue: TJclComplex): TJclComplex; overload;
    function CNewAdd(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNewDiv(const DivValue: TJclComplex): TJclComplex; overload;
    function CNewDiv(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNewMul(const MulValue: TJclComplex): TJclComplex; overload;
    function CNewMul(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNewSub(const SubValue: TJclComplex): TJclComplex; overload;
    function CNewSub(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNewNeg: TJclComplex;
    function CNewConjugate: TJclComplex;

    //----------- natural log and exponential functions
    function CLn: TJclComplex;
    function CNewLn: TJclComplex;
    function CExp: TJclComplex;
    function CNewExp: TJclComplex;
    function CPwr(const PwrValue: TJclComplex): TJclComplex; overload;
    function CPwr(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CNewPwr(PwrValue: TJclComplex): TJclComplex; overload;
    function CNewPwr(const X, Y: Float; const ComplexType: TComplexKind
      {$IFDEF SUPPORTS_DEFAULTPARAMS} = crRectangular {$ENDIF}): TJclComplex; overload;
    function CIntPwr(const Pwr: Integer): TJclComplex; overload;
    function CNewIntPwr(const Pwr: Integer): TJclComplex; overload;
    function CRealPwr(const Pwr: Float): TJclComplex; overload;
    function CNewRealPwr(const Pwr: Float): TJclComplex; overload;
    function CRoot(const K, N: Word): TJclComplex; overload;
    function CNewRoot(const K, N: Word): TJclComplex; overload;
    function CSqrt: TJclComplex; overload;
    function CNewSqrt: TJclComplex; overload;

    //----------- trigonometric functions
    function CCos: TJclComplex;
    function CNewCos: TJclComplex;
    function CSin: TJclComplex;
    function CNewSin: TJclComplex;
    function CTan: TJclComplex;
    function CNewTan: TJclComplex;
    function CCot: TJclComplex;
    function CNewCot: TJclComplex;
    function CSec: TJclComplex;
    function CNewSec: TJclComplex;
    function CCsc: TJclComplex;
    function CNewCsc: TJclComplex;

    //----------- complex hyperbolic functions
    function CCosH: TJclComplex;
    function CNewCosH: TJclComplex;
    function CSinH: TJclComplex;
    function CNewSinH: TJclComplex;
    function CTanH: TJclComplex;
    function CNewTanH: TJclComplex;
    function CCotH: TJclComplex;
    function CNewCotH: TJclComplex;
    function CSecH: TJclComplex;
    function CNewSecH: TJclComplex;
    function CCscH: TJclComplex;
    function CNewCscH: TJclComplex;

    //----------- complex Bessel functions of order zero
    function CI0: TJclComplex;
    function CNewI0: TJclComplex;
    function CJ0: TJclComplex;
    function CNewJ0: TJclComplex;

    function CApproxLnGamma: TJclComplex;
    function CNewApproxLnGamma: TJclComplex;
    function CLnGamma: TJclComplex;
    function CNewLnGamma: TJclComplex;
    function CGamma: TJclComplex;
    function CNewGamma: TJclComplex;

    //----------- miscellaneous routines
    function AbsoluteValue: Float; overload;
    function AbsoluteValue(const Coord: TRectCoord): Float; overload;
    function AbsoluteValueSqr: Float; overload;
    function AbsoluteValueSqr(const Coord: TRectCoord): Float; overload;
    function FormatExtended(const X: Float): string;
  end;

var
  ComplexPrecision: Float = 1E-14;

const
  MaxTerm: Byte = 35;
  EpsilonSqr: Float = 1E-20;

implementation

const
  MaxFracLen = 18;
  RectOne: TRectCoord = (X: 1.0; Y: 0.0);
  RectZero: TRectCoord = (X: 0.0; Y: 0.0);
  RectInfinity: TRectCoord = (X: Infinity; Y: Infinity);

//------------------------------------------------------------------------------

function Coordinates(const cX, cY: Float; CoordType: TComplexKind): TCoords;
begin
  case CoordType of
    crRectangular:
      begin
        Result.X := cX;
        Result.Y := cY;
        Result.R := 0.0;
        Result.Theta := 0.0;
      end;
    crPolar:
      begin
        Result.X := 0.0;
        Result.Y := 0.0;
        Result.R := cX;
        Result.Theta := cY;
      end;
  end;
end;

//------------------------------------------------------------------------------

function RectCoord(X, Y: Float): TRectCoord; overload;
begin
  Result.X := X;
  Result.Y := Y;
end;

//------------------------------------------------------------------------------

function RectCoord(Value: TJclComplex): TRectCoord; overload;
begin
  Result.X := Value.FCoord.X;
  Result.Y := Value.FCoord.Y;
end;

//------------------------------------------------------------------------------

constructor TJclComplex.Create;
begin
  inherited Create;
  AssignZero;
  FFracLen := MaxFracLen;
end;

//------------------------------------------------------------------------------

constructor TJclComplex.Create(const X, Y: Float; const ComplexType: TComplexKind);
begin
  inherited Create;

  Assign(X, Y, ComplexType);
  FFracLen := MaxFracLen;
end;

//------------------------------------------------------------------------------

procedure TJclComplex.FillCoords(const CompexType: TComplexKind);
begin
  MiscalcComplex;
  case CompexType of
    crPolar:
      begin
        FCoord.X := FCoord.R * Cos(FCoord.Theta);
        FCoord.Y := FCoord.R * Sin(FCoord.Theta);
      end;
    crRectangular:
      begin
        if FCoord.X = 0.0 then
        begin
          FCoord.R := Abs(FCoord.Y);
          FCoord.Theta := PiOn2 * Sgn(FCoord.Y);
        end
        else
        begin
          FCoord.R := Absolute