Documentation : Scilab is Not Naive

Documentation : Scilab is Not Naive Commit Details

Date:2010-12-20 11:32:04 (7 years 8 months ago)
Author:Michael Baudin
Commit:33
Parents: 32
Message:Fixed Stewart algorithm, with Corrigendum from TOMS/1986.
Changes:
M/en_US/scripts/cdiv.tst
M/en_US/scripts/cdiv_smithStewart.sci

File differences

en_US/scripts/cdiv.tst
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
// Copyright (C) 2010 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Some attempts to improve Smith's algorithm
haveimproved = %t;
// cdiv_naive.sci : a naive method
exec("cdiv_naive.sci");
// cdiv_smith.sci : original Smith's method
exec("cdiv_smith.sci");
// cdiv_smithStewart.sci : original Smith's method, improved by Steward
exec("cdiv_smithStewart.sci");
// A Scilab port of the ANSI/ISO C algorithm
exec("cdiv_ansiisoC.sci");
// A Scilab port of the Li et al. improved Smith
exec("cdiv_smithLi.sci");
// A division based on polar form
exec("cdiv_polar.sci");
// An improved Smith method, with ideas from Stewart (prodminmax)
exec("cdiv_smith2.sci");
// An improved Smith method, with an exp-log product
exec("cdiv_smith3.sci");
if ( haveimproved ) then
// An improved Smith method
exec("cdiv_smith4.sci");
// An improved Smith, with ideas from Li et al.
exec("cdiv_smith5.sci");
// An improved Smith, with ideas from Li et al.
exec("cdiv_smith6.sci");
end
exec("assert_computedigits.sci");
exec("assert_datasetread.sci");
ieee(2);
function [msg,digits] = fmtdig ( x , y , expected , divfunc , headermsg , footermsg )
// Performs the complex division x/y, compare to expected, and produve a message
// containing the number of digits in the real and imaginary parts.
r = divfunc ( a + %i * b , c + %i * d );
digits = assert_computedigits(r,expectedR,basis);
msg = msprintf("%s=%2d%s\n", headermsg , digits , footermsg )
endfunction
function r = cdiv_scilab ( x , y )
// Performs the complex division with Scilab /
r = x/y;
endfunction
////////////////////////////////////////////////////////////////////////
//
// Check accuracy
format("e",10);
path=pwd();
basis = 2;
msg=[];
dg=[];
dataset = fullfile(path,"cdiv.dataset.csv");
table = assert_datasetread ( dataset , "#" , "," , %t );
ntests = size(table,"r");
atable = evstr(table(:,1));
btable = evstr(table(:,2));
ctable = evstr(table(:,3));
dtable = evstr(table(:,4));
etable = evstr(table(:,5));
ftable = evstr(table(:,6));
for k = 1 : ntests
a = atable(k);
b = btable(k);
c = ctable(k);
d = dtable(k);
e = etable(k);
f = ftable(k);
//
expectedR = e+%i*f;
x = a + %i * b;
y = c + %i * d;
msg(1) = msprintf("Test #%3d/%d, %70s, ", k,ntests,sci2exp([a b c d e f]) );
dg(k,1)=0;
j=2;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_scilab , "Sc" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_naive , "Na" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith , "Sm" , ", " );
j=j+1;
// [msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smithStewart , "St" , ", " );
// j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith2 , "S2" , ", " );
//j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith3 , "S3" , ", " );
//j=j+1;
if ( haveimproved ) then
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith4 , "S4" , ", " );
j=j+1;
end
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_ansiisoC , "AC", ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smithLi , "Li" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith5 , "S5" , "" );
j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith6 , "S6" , ", " );
//j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_polar , "PO" , "" );
//j=j+1;
mprintf("%s\n",strcat(msg,""));
end
mprintf("Average:\n");
j=2;
mprintf("Sc=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Na=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Sm=%d\n",mean(dg(:,j)));
j=j+1;
//mprintf("ST=%d\n",mean(dg(:,j)));
//j=j+1;
//mprintf("S2=%d\n",mean(dg(:,j)));
//j=j+1;
//mprintf("S3=%d\n",mean(dg(:,j)));
//j=j+1;
mprintf("S4=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("AC=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Li=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("S5=%d\n",mean(dg(:,j)));
j=j+1;
//mprintf("S6=%d\n",mean(dg(:,j)));
//j=j+1;
//mprintf("PO=%d\n",mean(dg(:,j)));
//j=j+1;
// Discovered Issues:
// (0+%i*1 )/(0+%i*0) Scilab = Nan + Nan i Octave = NaN + Inf i Matlab = NaN + Inf i
// (1+%i*-1)/(0+%i*0) Scilab = Nan + Nan i Octave = Inf - Inf i Matlab = Inf - Inf i
// (1+%i*1)/(0+%i*0) Scilab = NaN + NaN i Octave = Inf + Inf i Matlab = Inf + Inf i
// cdiv_ansiisoC(1+%i*1,0+%i*0) = NaN + Infinity i
// This is because %inf + %inf * %i produces NaN + Infinity i. Is this right ?
// (%inf+%i*1)/(1+%i*1) Scilab = NaN - Infinity i Octave = Inf - Inf i, Matlab = Inf - Inf i
// but smith, ansiisoc = NaN - Infinity i
//
// smith3 - improved Smith method, with an exp-log product - does not improve much.
// Li improves for some case, but not all.
// Only smith5 works in all the current test cases.
//
// Copyright (C) 2010 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Some attempts to improve Smith's algorithm
haveimproved = %t;
// cdiv_naive.sci : a naive method
exec("cdiv_naive.sci");
// cdiv_smith.sci : original Smith's method
exec("cdiv_smith.sci");
// cdiv_smithStewart.sci : original Smith's method, improved by Steward
exec("cdiv_smithStewart.sci");
// A Scilab port of the ANSI/ISO C algorithm
exec("cdiv_ansiisoC.sci");
// A Scilab port of the Li et al. improved Smith
exec("cdiv_smithLi.sci");
// A division based on polar form
exec("cdiv_polar.sci");
// An improved Smith method, with ideas from Stewart (prodminmax)
exec("cdiv_smith2.sci");
// An improved Smith method, with an exp-log product
exec("cdiv_smith3.sci");
if ( haveimproved ) then
// An improved Smith method
exec("cdiv_smith4.sci");
// An improved Smith, with ideas from Li et al.
exec("cdiv_smith5.sci");
// An improved Smith, with ideas from Li et al.
exec("cdiv_smith6.sci");
end
exec("assert_computedigits.sci");
exec("assert_datasetread.sci");
ieee(2);
function [msg,digits] = fmtdig ( x , y , expected , divfunc , headermsg , footermsg )
// Performs the complex division x/y, compare to expected, and produve a message
// containing the number of digits in the real and imaginary parts.
r = divfunc ( a + %i * b , c + %i * d );
digits = assert_computedigits(r,expectedR,basis);
msg = msprintf("%s=%2d%s\n", headermsg , digits , footermsg )
endfunction
function r = cdiv_scilab ( x , y )
// Performs the complex division with Scilab /
r = x/y;
endfunction
////////////////////////////////////////////////////////////////////////
//
// Check accuracy
format("e",10);
path=pwd();
basis = 2;
msg=[];
dg=[];
dataset = fullfile(path,"cdiv.dataset.csv");
table = assert_datasetread ( dataset , "#" , "," , %t );
ntests = size(table,"r");
atable = evstr(table(:,1));
btable = evstr(table(:,2));
ctable = evstr(table(:,3));
dtable = evstr(table(:,4));
etable = evstr(table(:,5));
ftable = evstr(table(:,6));
for k = 1 : ntests
a = atable(k);
b = btable(k);
c = ctable(k);
d = dtable(k);
e = etable(k);
f = ftable(k);
//
expectedR = e+%i*f;
x = a + %i * b;
y = c + %i * d;
msg(1) = msprintf("Test #%3d/%d, %70s, ", k,ntests,sci2exp([a b c d e f]) );
dg(k,1)=0;
j=2;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_scilab , "Sc" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_naive , "Na" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith , "Sm" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smithStewart , "St" , ", " );
j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith2 , "S2" , ", " );
//j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith3 , "S3" , ", " );
//j=j+1;
if ( haveimproved ) then
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith4 , "S4" , ", " );
j=j+1;
end
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_ansiisoC , "AC", ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smithLi , "Li" , ", " );
j=j+1;
[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith5 , "S5" , "" );
j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_smith6 , "S6" , ", " );
//j=j+1;
//[msg(j),dg(k,j)] = fmtdig ( x , y , expectedR , cdiv_polar , "PO" , "" );
//j=j+1;
mprintf("%s\n",strcat(msg,""));
end
mprintf("Average:\n");
j=2;
mprintf("Sc=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Na=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Sm=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("ST=%d\n",mean(dg(:,j)));
j=j+1;
//mprintf("S2=%d\n",mean(dg(:,j)));
//j=j+1;
//mprintf("S3=%d\n",mean(dg(:,j)));
//j=j+1;
mprintf("S4=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("AC=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("Li=%d\n",mean(dg(:,j)));
j=j+1;
mprintf("S5=%d\n",mean(dg(:,j)));
j=j+1;
//mprintf("S6=%d\n",mean(dg(:,j)));
//j=j+1;
//mprintf("PO=%d\n",mean(dg(:,j)));
//j=j+1;
// Discovered Issues:
// (0+%i*1 )/(0+%i*0) Scilab = Nan + Nan i Octave = NaN + Inf i Matlab = NaN + Inf i
// (1+%i*-1)/(0+%i*0) Scilab = Nan + Nan i Octave = Inf - Inf i Matlab = Inf - Inf i
// (1+%i*1)/(0+%i*0) Scilab = NaN + NaN i Octave = Inf + Inf i Matlab = Inf + Inf i
// cdiv_ansiisoC(1+%i*1,0+%i*0) = NaN + Infinity i
// This is because %inf + %inf * %i produces NaN + Infinity i. Is this right ?
// (%inf+%i*1)/(1+%i*1) Scilab = NaN - Infinity i Octave = Inf - Inf i, Matlab = Inf - Inf i
// but smith, ansiisoc = NaN - Infinity i
//
// smith3 - improved Smith method, with an exp-log product - does not improve much.
// Li improves for some case, but not all.
// Only smith5 works in all the current test cases.
//
en_US/scripts/cdiv_smithStewart.sci
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
// Copyright (C) 2010 - Michael Baudin
// A Scilab port of Stewart's algorithm.
// References
//
// "A Note on Complex Division"
// G. W. Stewart, University of Maryland
// ACM Transactions on Mathematical Software, Vol. 11
// No 3, September 1985, Pages 238-241
//
// TODO : fix this using :
// G. W. Stewart
// Corrigendum: ``{A} Note on Complex Division''",
// TOMS,
// volume 12
// number 3
// pages = 285--285
// month = September, 1986
function r = cdiv_smithStewart ( x,y )
error("Does not work")
c = real(x)
d = imag(x)
e = real(y)
f = imag(y)
flip = %f
if ( abs(f) >= abs(e) ) then
[f,e]=switch(e,f)
[c,d]=switch(d,c)
flip=%t
end
s=1/e
t=1/(e+f*(f*s))
if ( abs(f) >= abs(s) ) then
[f,s]=switch(s,f)
end
if ( abs(d) >= abs(s) ) then
a=t*(c+s*(d*f))
elseif ( abs(d) >= abs(f) ) then
a=t*(c+d*(s*f))
else
a=t*(c+f*(s*d))
end
if ( abs(c) >= abs(s) ) then
b=t*(d+s*(c*f))
elseif ( abs(c) >= abs(f) ) then
b=t*(d+c*(s*f))
else
b=t*(d+f*(s*c))
end
if ( flip ) then
b=-b
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(a,b)
endfunction
function [c,d]=switch(a,b)
c=a
d=b
endfunction
function r = cdiv_smithStewart_Other ( x,y )
c = real(x)
d = imag(x)
e = real(y)
f = imag(y)
flip = %f
if ( abs(f) >= abs(e) ) then
// Switch e and f
tmp=e
e=f
f=tmp
// Switch c and d
tmp=c
c=d
d=tmp
flip=%t
end
s=1/e
t=1/(e+f*(f*s))
if ( abs(f) >= abs(s) ) then
// Switch f and s
tmp=f
f=s
s=tmp
end
if ( abs(d) >= abs(s) ) then
a=t*(c+s*(d*f))
elseif ( abs(d) >= abs(f) ) then
a=t*(c+d*(s*f))
else
a=t*(c+f*(s*d))
end
if ( abs(c) >= abs(s) ) then
b=t*(d+s*(c*f))
elseif ( abs(c) >= abs(f) ) then
b=t*(d+c*(s*f))
else
b=t*(d+f*(s*c))
end
if ( flip ) then
b=-b
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(a,b)
endfunction
// Copyright (C) 2010 - Michael Baudin
// A Scilab port of Stewart's algorithm.
// References
//
// "A Note on Complex Division"
// G. W. Stewart, University of Maryland
// ACM Transactions on Mathematical Software, Vol. 11
// No 3, September 1985, Pages 238-241
//
// G. W. Stewart
// Corrigendum: ``{A} Note on Complex Division''",
// TOMS, Volume 12, Number 3, Pages = 285--285, September, 1986
function r = cdiv_smithStewart ( x,y )
c = real(x)
d = imag(x)
e = real(y)
f = imag(y)
flip = %f
if ( abs(f) >= abs(e) ) then
[f,e]=switch(e,f)
[c,d]=switch(d,c)
flip=%t
end
s=1/e
t=1/(e+f*(f*s))
if ( abs(f) >= abs(s) ) then
[f,s]=switch(s,f)
end
if ( abs(d) >= abs(s) ) then
a=t*(c+s*(d*f))
elseif ( abs(d) >= abs(f) ) then
a=t*(c+d*(s*f))
else
a=t*(c+f*(s*d))
end
if ( abs(c) >= abs(s) ) then
b=t*(d-s*(c*f))
elseif ( abs(c) >= abs(f) ) then
b=t*(d-c*(s*f))
else
b=t*(d-f*(s*c))
end
if ( flip ) then
b=-b
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(a,b)
endfunction
function [c,d]=switch(a,b)
c=a
d=b
endfunction
function r = cdiv_smithStewart_Other ( x,y )
c = real(x)
d = imag(x)
e = real(y)
f = imag(y)
flip = %f
if ( abs(f) >= abs(e) ) then
// Switch e and f
tmp=e
e=f
f=tmp
// Switch c and d
tmp=c
c=d
d=tmp
flip=%t
end
s=1/e
t=1/(e+f*(f*s))
if ( abs(f) >= abs(s) ) then
// Switch f and s
tmp=f
f=s
s=tmp
end
if ( abs(d) >= abs(s) ) then
a=t*(c+s*(d*f))
elseif ( abs(d) >= abs(f) ) then
a=t*(c+d*(s*f))
else
a=t*(c+f*(s*d))
end
if ( abs(c) >= abs(s) ) then
b=t*(d+s*(c*f))
elseif ( abs(c) >= abs(f) ) then
b=t*(d+c*(s*f))
else
b=t*(d+f*(s*c))
end
if ( flip ) then
b=-b
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(a,b)
endfunction

Archive Download the corresponding diff file

Revision: 33