← Back to branch summary

~ubuntu-branches/ubuntu/utopic/mricron/utopic

« back to all changes in this revision

Viewing changes to fpmath/upolutil.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
fpmath/upolutil.pas
 
1
{ ******************************************************************
 
2
  Utility functions to handle roots of polynomials
 
3
  ****************************************************************** }
 
4
 
 
5
unit upolutil;
 
6
 
 
7
interface
 
8
 
 
9
uses
 
10
  utypes, uminmax;
 
11
 
 
12
function SetRealRoots(Deg : Integer; Z : PCompVector; Tol : Float) : Integer;
 
13
{ ------------------------------------------------------------------
 
14
  Set the imaginary part of a root to zero if it is less than a
 
15
  fraction Tol of its real part. This root is therefore considered
 
16
  real. The function returns the total number of real roots.
 
17
  ------------------------------------------------------------------ }
 
18
 
 
19
procedure SortRoots(Deg : Integer; Z : PCompVector);
 
20
{ ------------------------------------------------------------------
 
21
  Sort roots so that:
 
22
 
 
23
  (1) The Nr real roots are stored in elements [1..Nr] of vector Z,
 
24
      in increasing order.
 
25
 
 
26
  (2) The complex roots are stored in elements [(Nr + 1)..Deg] of
 
27
      vector Z and are unordered.
 
28
  ------------------------------------------------------------------ }
 
29
 
 
30
implementation
 
31
 
 
32
function SetRealRoots(Deg : Integer; Z : PCompVector; Tol : Float) : Integer;
 
33
var
 
34
  I, N : Integer;
 
35
begin
 
36
  for I := 1 to Deg do
 
37
    if (Z^[I].Y <> 0.0) and (Abs(Z^[I].Y) < Tol * Abs(Z^[I].X)) then
 
38
      Z^[I].Y := 0.0;
 
39
 
 
40
  { Count real roots }
 
41
  N := 0;
 
42
  for I := 1 to Deg do
 
43
    if Z^[I].Y = 0.0 then
 
44
      Inc(N);
 
45
 
 
46
  SetRealRoots := N;
 
47
end;
 
48
 
 
49
procedure SortRoots(Deg : Integer; Z : PCompVector);
 
50
var
 
51
  I, J, K, Nr, Nc : Integer;
 
52
  R, X, Y         : PVector;
 
53
 
 
54
  procedure Sort(X : PVector; N : Integer);
 
55
  { Sort vector X (insertion sort) }
 
56
  var
 
57
    I, J, K : Integer;
 
58
    A       : Float;
 
59
  begin
 
60
    for I := 1 to Pred(N) do
 
61
      begin
 
62
        K := I;
 
63
        A := X^[I];
 
64
        for J := Succ(I) to N do
 
65
          if X^[J] < A then
 
66
            begin
 
67
              K := J;
 
68
              A := X^[J];
 
69
            end;
 
70
        FSwap(X^[I], X^[K]);
 
71
      end;
 
72
  end;
 
73
 
 
74
begin
 
75
  { Count real and complex roots }
 
76
  Nr := 0; Nc := 0;
 
77
  for I := 1 to Deg do
 
78
    if Z^[I].Y = 0.0 then Inc(Nr) else Inc(Nc);
 
79
 
 
80
  DimVector(R, Nr);
 
81
  DimVector(X, Nc);
 
82
  DimVector(Y, Nc);
 
83
 
 
84
  { Store real roots in R and complex roots in (X,Y) }
 
85
  J := 0; K := 0;
 
86
  for I := 1 to Deg do
 
87
    if Z^[I].Y = 0.0 then
 
88
      begin
 
89
        Inc(J);
 
90
        R^[J] := Z^[I].X;
 
91
      end
 
92
    else
 
93
      begin
 
94
        Inc(K);
 
95
        X^[K] := Z^[I].X;
 
96
        Y^[K] := Z^[I].Y;
 
97
      end;
 
98
 
 
99
  { Sort vector R (insertion sort) }
 
100
  if Nr > 0 then Sort(R, Nr);
 
101
 
 
102
  { Transfer real roots into elements 1..Nr }
 
103
  for I := 1 to Nr do
 
104
    begin
 
105
      Z^[I].X := R^[I];
 
106
      Z^[I].Y := 0.0;
 
107
    end;
 
108
 
 
109
  { Transfer complex roots into elements (Nr+1)..Deg }
 
110
  for I := 1 to Nc do
 
111
    begin
 
112
      J := I + Nr;
 
113
      Z^[J].X := X^[I];
 
114
      Z^[J].Y := Y^[I];
 
115
    end;
 
116
 
 
117
  DelVector(R, Nr);
 
118
  DelVector(X, Nc);
 
119
  DelVector(Y, Nc);
 
120
end;
 
121
 
 
122
end.