CplxMath.pas

   1 unit CplxMath;
   2 {* This unit contains functins for working with complex numbers. To use with
   3    KOL library and its kolmath.pas unit instead of standard math.pas, define
   4    synmbol KOL in project options, or uncomment its definition below. }
   5
   6 interface
   7
   8 //{$DEFINE KOL}
   9
  10 {$IFNDEF KOL}
  11   {$IFDEF KOL_MCK}
  12     {$DEFINE KOL}
  13   {$ENDIF}
  14 {$ENDIF}
  15
  16 uses {$IFDEF KOL} kolmath, kol {$ELSE} math, sysutils {$ENDIF};
  17
  18 type
  19   Complex = record Re, Im: double end;
  20   {* }
  21
  22   function CfromReIm( Re, Im: Double ): Complex;
  23   {* Re + i * Im }
  24
  25   function Cadd( const X, Y: Complex ): Complex;
  26   {* X + Y }
  27
  28   function Cneg( const X: Complex ): Complex;
  29   {* -X }
  30
  31   function Csub( const X, Y: Complex ): Complex;
  32   {* X - Y }
  33
  34   function Cmul( const X, Y: Complex ): Complex;
  35   {* X * Y }
  36
  37   function CmulD( const X: Complex; D: Double ): Complex;
  38   {* X * D }
  39
  40   function CmulI( const X: Complex ): Complex;
  41   {* i * X }
  42
  43   function Cdiv( const X, Y: Complex ): Complex;
  44   {* X / Y }
  45
  46   function Cmod( const X: Complex ): Double;
  47   {* Q( X.Re^2 + X.Im^2 ) }
  48
  49   function Carg( const X: Complex ): Double;
  50   {* arctg( X.Im / X.Re ) }
  51
  52   function CfromModArg( R, Arg: Double ): Complex;
  53   {* R * ( cos Arg + i * sin Arg ) }
  54
  55   function Cpow( const X: Complex; Pow: Double ): Complex;
  56   {* X ^ Pow }
  57
  58   function Csqrt( const X: Complex ): Complex;
  59   {* Q( X ) }
  60
  61   function Cexp( const X: Complex ): Complex;
  62   {* exp( X ) }
  63
  64   function Cln( const X: Complex ): Complex;
  65   {* ln( X ) }
  66
  67   function Ccos( const X: Complex ): Complex;
  68   {* cos( X ) }
  69
  70   function Csin( const X: Complex ): Complex;
  71   {* sin( X ) }
  72
  73   function C2Str( const X: Complex ): String;
  74   {* }
  75
  76   function C2StrEx( const X: Complex ): String;
  77   {* experimental }
  78
  79 implementation
  80
  81   function CfromReIm( Re, Im: Double ): Complex;
  82   begin
  83     Result.Re := Re;
  84     Result.Im := Im;
  85   end;
  86
  87   function Cadd( const X, Y: Complex ): Complex;
  88   begin
  89     Result.Re := X.Re + Y.Re;
  90     Result.Im := X.Im + Y.Im;
  91   end;
  92
  93   function Cneg( const X: Complex ): Complex;
  94   begin
  95     Result.Re := -X.Re;
  96     Result.Im := -X.Im;
  97   end;
  98
  99   function Csub( const X, Y: Complex ): Complex;
 100   begin
 101     Result := Cadd( X, Cneg( Y ) );
 102   end;
 103
 104   function Cmul( const X, Y: Complex ): Complex;
 105   begin
 106     Result.Re := X.Re * Y.Re - X.Im * Y.Im;
 107     Result.Im := X.Re * Y.Im + X.Im * Y.Re;
 108   end;
 109
 110   function CmulD( const X: Complex; D: Double ): Complex;
 111   begin
 112     Result.Re := X.Re * D;
 113     Result.Im := X.Im * D;
 114   end;
 115
 116   function CmulI( const X: Complex ): Complex;
 117   begin
 118     Result.Re := -X.Im;
 119     Result.Im := X.Re;
 120   end;
 121
 122   function Cdiv( const X, Y: Complex ): Complex;
 123   var Z: Double;
 124   begin
 125     Z := 1.0 / ( Y.Re * Y.Re + Y.Im * Y.Im );
 126     Result.Re := (X.Re * Y.Re + X.Im * Y.Im ) * Z;
 127     Result.Im := (X.Im * Y.Re - X.Re * Y.Im ) * Z;
 128   end;
 129
 130   function Cmod( const X: Complex ): Double;
 131   begin
 132     Result := sqrt( X.Re * X.Re + X.Im * X.Im );
 133   end;
 134
 135   function Carg( const X: Complex ): Double;
 136   begin
 137     Result := ArcTan2( X.Im, X.Re );
 138   end;
 139
 140   function CfromModArg( R, Arg: Double ): Complex;
 141   begin
 142     Result.Re := R * cos( Arg );
 143     Result.Im := R * sin( Arg );
 144   end;
 145
 146   function Cpow( const X: Complex; Pow: Double ): Complex;
 147   var R, A: Double;
 148   begin
 149     R := power( Cmod( X ), Pow );
 150     A := Pow * Carg( X );
 151     Result := CfromModArg( R, A );
 152   end;
 153
 154   function Csqrt( const X: Complex ): Complex;
 155   begin
 156     Result := Cpow( X, 0.5 );
 157   end;
 158
 159   function Cexp( const X: Complex ): Complex;
 160   var Z: Double;
 161   begin
 162     Z := exp( X.Re );
 163     Result.Re := Z * cos( X.Im );
 164     Result.Im := Z * sin( X.Im );
 165   end;
 166
 167   function Cln( const X: Complex ): Complex;
 168   begin
 169     Result := CfromModArg( ln( Cmod( X ) ), Carg( X ) );
 170   end;
 171
 172   function Ccos( const X: Complex ): Complex;
 173   begin
 174     Result := CmulI( X );
 175     Result := CmulD( Cadd( Cexp( Result ), Cexp( Cneg( Result ) ) ),
 176                      0.5 );
 177   end;
 178
 179   function Csin( const X: Complex ): Complex;
 180   begin
 181     Result := CmulI( X );
 182     Result := CmulD( Csub( Cexp(Result), Cexp( Cneg(Result) ) ),
 183                      0.5 );
 184   end;
 185
 186   {$IFDEF KOL}
 187   function Abs( X: Double ): Double;
 188   begin
 189     Result := EAbs( X );
 190   end;
 191   {$ENDIF}
 192
 193   {$IFNDEF KOL}
 194   function Double2Str( D: Double ): String;
 195   begin
 196     Result := DoubleToStr( D );
 197   end;
 198   {$ENDIF}
 199
 200   function C2Str( const X: Complex ): String;
 201   begin
 202     if Abs( X.Im ) < 1e-307 then
 203     begin
 204       Result := Double2Str( X.Re );
 205     end
 206       else
 207     begin
 208       Result := '';
 209       if Abs( X.Re ) > 1e-307 then
 210       begin
 211         Result := Double2Str( X.Re );
 212         if X.Im > 0.0 then
 213           Result := Result + ' + ';
 214       end;
 215       if X.Im < 0.0 then
 216         Result := Result + '- i * ' + Double2Str( -X.Im )
 217       else
 218         Result := Result + 'i * ' + Double2Str( X.Im );
 219     end;
 220   end;
 221
 222   function C2StrEx( const X: Complex ): String;
 223   begin
 224     if Abs( X.Im ) < 1e-307 then
 225     begin
 226       Result := Double2StrEx( X.Re );
 227     end
 228       else
 229     begin
 230       Result := '';
 231       if Abs( X.Re ) > 1e-307 then
 232       begin
 233         Result := Double2StrEx( X.Re );
 234         if X.Im > 0.0 then
 235           Result := Result + ' + ';
 236       end;
 237       if X.Im < 0.0 then
 238         Result := Result + '- i * ' + Double2StrEx( -X.Im )
 239       else
 240         Result := Result + 'i * ' + Double2StrEx( X.Im );
 241     end;
 242   end;
 243
 244 end.