Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

[ec6234]: mp_TESTING2.m Maximize Restore History

Download this file

mp_TESTING2.m    92 lines (82 with data), 4.0 kB

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
% Thanks for Carlos Lopez for this file and many other enhancements to the mp toolbox.
%quick and dirty tests
precision=1000;%notice that tol~=(.5)^precision;
a42=mp(rand(4,2),precision);%It will be nice to use mprand, but it doesn't work
b42=mp(rand(4,2),precision);
%2) direct/inverse pairs
c_mp=norm(asin(sin(a42))-a42);disp(['2.1 sin/asin' setstr(9) num2str(c_mp)])
c_mp=norm(acos(cos(a42))-a42);disp(['2.2 cos/acos' setstr(9) num2str(c_mp)])
c_mp=norm(atan(tan(a42))-a42);disp(['2.3 tan/atan' setstr(9) num2str(c_mp)])
c_mp=norm(asinh(sinh(a42))-a42);disp(['2.4 sinh/asinh' setstr(9) num2str(c_mp)])
c_mp=norm(acosh(cosh(a42))-a42);disp(['2.5 cosh/acosh' setstr(9) num2str(c_mp)])
c_mp=norm(atanh(tanh(a42))-a42);disp(['2.6 tanh/atanh' setstr(9) num2str(c_mp)])
c_mp=norm(exp(log(a42))-a42);disp( ['2.7 exp/log ' setstr(9) num2str(c_mp)])
c_mp=norm(sqrt(a42.^2)-a42);disp( ['2.8 sqrt/^2 ' setstr(9) num2str(c_mp)])
%4) Catastrophic cancellation
z=sqrt(117)-sqrt(116);zExact=1/(sqrt(117)+sqrt(116));
disp('4.1 Improvement in the catastrophic cancellation of sqrt(117)-sqrt(116)')
disp([setstr(9) 'Standard double precision:' num2str(abs(z/zExact-1)) ])
for precision=300:100:1000
zmp=sqrt(mp(117,precision))-sqrt(mp(116,precision));zmpExact=1/(sqrt(mp(117,precision))+sqrt(mp(116,precision)));
disp([setstr(9) 'Multiple precision = ' num2str(precision) ':' num2str(abs((zmp/zmpExact-1)))])
end
disp('4.2 Accuracy of tan(pi/4)-1')
for precision=300:100:1000
zmp=mp('pi',precision)/4;
disp([setstr(9) 'Multiple precision = ' num2str(precision) ':' num2str(((tan(zmp)-1)))])
end
disp('4.3 Accuracy of 4*atan2(-1,1)+pi')
for precision=300:100:1000
zmp=mp('pi',2*precision);%Request a higher precision pi in order to compare
one=mp(1,precision);
disp([setstr(9) 'Multiple precision = ' num2str(precision) ':' ...
num2str(((atan2(-one,one)*4+zmp)))])
end
%5) Test matrix; we will only test that the accuracy
%of an exact inverse times the matrix is closer to the identity in mp rather than in double
n=14;
disp(['5.1 Hilbert matrix and its exact inverse.'])
disp([' norm in double: ' num2str([norm(hilb(n)*invhilb(n)-eye(n))])])
for precision=300:100:1000
mp_set_defaults(precision);
%5.1 Hilbert
one=mp(1);
J = 1:n;
J = J(ones(n,1),:);
I = J';
E = mp(ones(n,n));
mpHilb = E./(I+J-one);
%Now the inverse; this is a slightly modified version of invhilb.m to produce mp values
p = n;
H = mp(zeros(n,n));
for k = 1:n
i=mp(k);
if k > 1, p = ((n-i+one)*p*(n+i-one))/(k-one)^2; end
r = p*p;
H(k,k) = r/(2*i-one);
for j = i+1:n
% r = -((n-j+one)*r*(n+j-one))/(j-one)^2;
num = -((n-j+one)*r*(n+j-one));
den= (j-one)^2;
r = num/den;
H(k,j) = r/(i+j-one);
H(j,k) = r/(i+j-one);
end
end
disp(['precision=' num2str(precision) ':' setstr(9) num2str(norm(H*mpHilb-eye(n)))]);
end
%check also some auxiliary functions (unrelated with precision!)
X = mp([2 8 4;7 3 9]);
%see "help min" for full documentation
if all(min(X,[],1) == [2 3 4]), disp('Test 6.1 passed'), else, disp('Test 6.1 failed'),end
if all(min(X,[],2) == [2;3]), disp('Test 6.2 passed'), else, disp('Test 6.2 failed'),end
if all(min(X,5)==[2 5 4;5 3 5]), disp('Test 6.3 passed'), else, disp('Test 6.3 failed'),end
if all(max(X,[],1)==[7 8 9]), disp('Test 6.4 passed'), else, disp('Test 6.4 failed'),end
if all(max(X,[],2)==[8;9]), disp('Test 6.5 passed'), else, disp('Test 6.5 failed'),end
if all(max(X,5)==[5 8 5;7 5 9]),disp('Test 6.6 passed'), else, disp('Test 6.6 failed'),end
X = mp([0 1 2;3 4 5]);
if all(sum(X,1) == [3 5 7]), disp('Test 6.7 passed'), else, disp('Test 6.7 failed'),end
if all(sum(X,2) == [ 3;12]), disp('Test 6.8 passed'), else, disp('Test 6.8 failed'),end
X = mp([3 7 5;0 4 2]);
if all(sort(X,1) == [0 4 2;3 7 5] ), disp('Test 6.9 passed'), else, disp('Test 6.9 failed'),end
if all(sort(X,2) == [3 5 7;0 2 4] ), disp('Test 6.10 passed'), else, disp('Test 6.10 failed'),end