Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/intfact.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % INTFACT computes integer factorisation
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % INTFACT tries to find two integers, P and Q, that satisfy
5 %
6 % y = P/Q * x
7 %
8 % >> [p,q] = intfact(y,x)
9 %
10 % The following call returns a parameter list object that contains the
11 % default parameter values:
12 %
13 % >> pl = intfact(utils, 'Params')
14 %
15 % The following call returns a string that contains the routine CVS version:
16 %
17 % >> version = intfact(utils,'Version')
18 %
19 % The following call returns a string that contains the routine category:
20 %
21 % >> category = intfact(utils,'Category')
22 %
23 % VERSION: $Id: intfact.m,v 1.9 2009/02/06 12:03:52 hewitson Exp $
24 %
25 % M Hewitson 26-05-08
26 %
27 %
28 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29
30 function varargout = intfact(varargin)
31
32 % Get input sample rates
33 fs2 = floor(1e10*varargin{1});
34 fs1 = floor(1e10*varargin{2});
35
36 g = gcd(fs2,fs1);
37
38 varargout{1} = fs2/g;
39 varargout{2} = fs1/g;
40
41 return
42
43 % Get min and max fs
44 x = max(fs1,fs2);
45 y = min(fs1,fs2);
46
47 % If we already have two integers then it's easy
48 if rem(fs2,1) == 0 && rem(fs1,1) == 0
49 P = fs2;
50 Q = fs1;
51 % Get the greatest common divisor
52 g = gcd(P,Q);
53 % factor that out
54 P = P / g;
55 Q = Q / g;
56 % If fs1 and fs2 are factors of each other, then it's easy
57 elseif rem(fs1,fs2) == 0
58 P = fs1/fs2;
59 Q = 1;
60 elseif rem(fs2,fs1) == 0
61 Q = fs2/fs1;
62 P = 1;
63 % Otherwise we search for two integers
64 else
65 % Start from N=0
66 N = 0;
67 s = x * 10^N;
68 % Now compute initial P and Q
69 if fs2 > fs1
70 P = 10^N;
71 Q = s / y;
72 else
73 Q = 10^N;
74 P = s / y;
75 end
76 % Now look for the integers
77 while rem(s,1) > 0 || rem(P,1) > 0 || rem(Q,1) > 0
78 s = x * 10^N;
79 % Now compute P and Q
80 if fs2 > fs1
81 P = 10^N;
82 Q = s / y;
83 else
84 Q = 10^N;
85 P = s / y;
86 end
87 N = N+1;
88 end
89
90 % Get the greatest common divisor
91 g = gcd(P,Q);
92
93 % factor that out
94 P = P / g;
95 Q = Q / g;
96 end
97
98 % Set outputs
99 if nargout == 2
100 varargout{1} = P;
101 varargout{2} = Q;
102 else
103 error('### Incorrect number of outputs');
104 end
105
106 end
107 % END