Manage a rational number without loss of precision in Matlab?

Question

I want to use this rational number in calculations without losing image accuracy in Matlab:

f = 359.0 + 16241/16250.0 

I think that saving, for example, f = uint64(359.0 + 16241/16250.0) loses the accuracy visible as 360 in Matlab.

I think the best way to handle this is to never store the value, but store its factors, such as

 % f = a + b/c a = 359 b = 16241 c = 16250 

and then the calculation is performed on the variables a, b and c and giving the result as an image.

Is this a good way to maintain accuracy?

Answer 1

As you believe, if you absolutely do not want to lose accuracy when storing a rational number, the best solution is probably to save the number in terms of its integer components.

Instead of three components ( f = a + b/c ), you can reduce the reprint by two components: f = n/d . Thus, each rational number will be defined (and stored) as a two-component integer vector [nd] . For example, the number f in your example corresponds to n=5849991 and d=16250 .

To simplify the processing of rational numbers stored in this way, you can define a helper function that converts from the representation [nd] to n/d before applying the desired operation:

 useInteger = @(x, nd, fun) fun(x,double(nd(1))/double(nd(2))); 

Then

 >> x = sqrt(pi); >> nd = int64([5849991 16250]); >> useInteger(x, nd, @plus) ans = 361.7719 >> useInteger(x, nd, @times) ans = 638.0824 

If you want to achieve arbitrary high precision calculations, you should use variable precision arithmetic ( vpa ) with string arguments. With this approach, you can specify how many digits you want:

 >> vpa('sqrt(pi)*5849991/16250', 50) ans = 638.08240465923757600307902117159072301901656248436 

Answer 2

Perhaps create a Rational class and define the necessary operations ( plus , minus , times , etc.). Start with something like this:

Rational.m

 classdef Rational properties n; d; end methods function obj = Rational(n,d) GCD = gcd(n,d); obj.n = n./GCD; obj.d = d./GCD; end function d = dec(obj) d = double(obj.n)/double(obj.d); end % X .* Y function R = times(X,Y) chkxy(X,Y); if isnumeric(X), N = X .* Yn; D = Yd; elseif isnumeric(Y), N = Xn .* Y; D = Xd; else N = Xn .* Yn; D = Xd .* Yd; end R = Rational(N,D); end % X * Y function R = mtimes(X,Y) R = times(X,Y); end % X ./ Y function R = rdivide(X,Y) if isnumeric(Y), y = Rational(1,Y); else y = Rational(Yd,Yn); end R = times(X,y); end % X / Y function R = mrdivide(X,Y) R = rdivide(X,Y); end % X + Y function R = plus(X,Y) chkxy(X,Y); if isnumeric(X), N = X.*Yd + Yn; D = Yd; elseif isnumeric(Y), N = Y.*Xd + Xn; D = Xd; else D = lcm(Xd,Yd); N = sum([Xn Yn].*(D./[Xd Yd])); end R = Rational(N,D); end % X - Y function R = minus(X,Y) R = plus(X,-Y); end % -X function R = uminus(X) R = Rational(-Xn,Xd); end function chkxy(X,Y) if (~isa(X, 'Rational') && ~isnumeric(X)) || ... (~isa(Y, 'Rational') && ~isnumeric(Y)), error('X and Y must be Rational or numeric.'); end end end end 

<strong> Examples

Build Objects:

 >> clear all % reset class definition >> r1 = Rational(int64(1),int64(2)) r1 = Rational with properties: n: 1 d: 2 >> r2 = Rational(int64(3),int64(4)) r2 = Rational with properties: n: 3 d: 4 

Add and Subtract:

 >> r1+r2 ans = Rational with properties: n: 5 d: 4 >> r1-r2 ans = Rational with properties: n: -1 d: 4 

Multiply and divide:

 >> r1*r2 ans = Rational with properties: n: 3 d: 8 >> r1/r2 ans = Rational with properties: n: 2 d: 3 

Get decimal value:

 >> r12 = r1/r2; % 2/3 ((1/2)/(3/4)) >> f = r12.dec f = 0.6667 

Answer 3

LuisMendo answer extension

I got this as an error for your py suggestion

 >>> a = 638.08240465923757600307902117159072301901656248436059 >>> a 638.0824046592376 % do not know if Python is computing here with exact number >>> b = 638.0824 >>> ave = abs(b+a)/2 >>> diff = abs(ba) >>> ave = abs(b+a)/2 >>> diff/ave 7.30193709165014e-09 

which is greater than the estimated error storing error above.

I run in WolframAlpha

 x = sqrt(pi) x*5849991/16250 

and get

 509.11609919757198016211937362635174599076143654820109 

I'm not sure if this is what you meant in your comment on your answer.

Answer 4

Extension for chappjc answer.

I have now

 [B,T,F] = tfrwv(data1, 1:length(data1), length(data1)); % here F double fs = Rational(uint64(5849991), uint64(16250)); t = 1/fs; imagesc(T*t, F*fs, B); 

I run it

 Error using .* Integers can only be combined with integers of the same class, or scalar doubles. Error in .* (line 23) N = X .* Yn; D = Yd; Error in * (line 34) R = times(X,Y); 

How can you multiply double in this class with Rational?