Documentation : Scilab is Not Naive

Documentation : Scilab is Not Naive Commit Details

Date:2010-12-21 22:50:06 (7 years 7 months ago)
Author:Michael Baudin
Commit:34
Parents: 33
Message:Added licence header. Added C version of Priest.
Changes:
A/en_US/scripts/cdiv_priest.c
A/en_US/scripts/cdiv_priest-driver.c
M/en_US/scripts/cdiv_smith3.sci
M/en_US/scripts/cdiv_smithLi.sci
M/en_US/scripts/assert_computedigits.sci
M/en_US/scripts/cdiv_polar.sci
M/en_US/scripts/cdiv_smithStewart.sci
M/en_US/scripts/complexdivision.sce
M/en_US/scripts/cdiv_smith.sci
M/en_US/scripts/cdiv_smith2.sci

File differences

en_US/scripts/assert_computedigits.sci
11
2
3
4
5
6
7
28
39
410
// Copyright (C) 2009 - 2010 - 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
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
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
// 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
// Copyright (C) 2010 - 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
// 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
en_US/scripts/cdiv_priest-driver.c
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
// Copyright (C) 2010 - 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
// References
//
// "Efficient scaling for complex division"
// Douglas M. Priest
// Sun Microsystems
// ACM Transactions on Mathematical software, Vol. 30, No. 4, December 2004
// A driver for Priest's division
/*
To compile :
$ gcc -c cdiv_priest.c
$ gcc -c cdiv_priest-driver.c
$ gcc -o priest.exe cdiv_priest.o cdiv_priest-driver.o
$ ./priest.exe
(1.000000e+00+I*2.000000e+00)/(3.000000e+00+I*4.000000e+00)=4.400000e-01+I*8.000000e-02
*/
#include <stdio.h>
#include <stdlib.h>
// Set v[0] + i v[1] := (z[0] + i z[1]) / (w[0] + i w[1]) assuming:
void cdiv(double *v, const double *z, const double *w);
int main(void)
{
double v[2];
double z[2];
double w[2];
// Perform (1+i*2)/(3+i*4)
z[0]=1;
z[1]=2;
w[0]=3;
w[1]=4;
cdiv(v, z, w);
printf("(%e+I*%e)/(%e+I*%e)=%e+I*%e\n", z[0],z[1],w[0],w[1],v[0],v[1]);
return EXIT_SUCCESS;
}
en_US/scripts/cdiv_smithLi.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
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
// Copyright (C) 2010 - Michael Baudin
// An improved Smith method, from Li et al.
// Li,, Xiaoye S. and Demmel,, James W. and Bailey,, David H. and Henry,, Greg and Hida,, Yozo and Iskandar,, Jimmy and Kahan,, William and Kang,, Suh Y. and Kapur,, Anil and Martin,, Michael C. and Thompson,, Brandon J. and Tung,, Teresa and Yoo,, Daniel J.},
// Design, implementation and testing of extended and mixed precision BLAS},
// ACM Trans. Math. Softw., vol. 28, numb. 2, 2002
// Tested in the normal range.
// TODO : test in the scaled cases.
function q = cdiv_smithLi ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
AB = max(abs([a b]))
CD = max(abs([c d]))
// The choice of Li et al. for B:
B = 2
S = 1
OV = number_properties("huge")
UN = number_properties("tiny")
Be = B/%eps^2
// Scaling
if ( AB > OV/16 ) then
// Scale down a, b
a = a/16
b = b/16
S = S*16
end
if ( CD > OV/16 ) then
// Scale down c, d
c = c/16
d = d/16
S = S/16
end
if ( AB < UN*B/%eps ) then
// Scale up a, b
a = a * Be
b = b * Be
S = S/Be
end
if ( CD < UN*B/%eps ) then
// Scale up c, d
c = c * Be
d = d * Be
// The paper has a mistake there
S = S*Be
end
// Now a, b, c, d in [UN*B/%eps,OV/16]
if ( abs(c) > abs(d) ) then
r = d/c
t = 1/(c+d*r)
qre = (a+b*r)*t
qim = (b-a*r)*t
else
r = c/d
t = 1/(d+c*r)
qre = (b+a*r)*t
qim = (-a+b*r)*t
end
// Scale back
qre = qre * S
qim = qim * S
q = complex(qre,qim)
endfunction
// Copyright (C) 2010 - 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
// An improved Smith method, from Li et al.
// Li,, Xiaoye S. and Demmel,, James W. and Bailey,, David H. and Henry,, Greg and Hida,, Yozo and Iskandar,, Jimmy and Kahan,, William and Kang,, Suh Y. and Kapur,, Anil and Martin,, Michael C. and Thompson,, Brandon J. and Tung,, Teresa and Yoo,, Daniel J.},
// Design, implementation and testing of extended and mixed precision BLAS},
// ACM Trans. Math. Softw., vol. 28, numb. 2, 2002
// Tested in the normal range.
// TODO : test in the scaled cases.
function q = cdiv_smithLi ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
AB = max(abs([a b]))
CD = max(abs([c d]))
// The choice of Li et al. for B:
B = 2
S = 1
OV = number_properties("huge")
UN = number_properties("tiny")
Be = B/%eps^2
// Scaling
if ( AB > OV/16 ) then
// Scale down a, b
a = a/16
b = b/16
S = S*16
end
if ( CD > OV/16 ) then
// Scale down c, d
c = c/16
d = d/16
S = S/16
end
if ( AB < UN*B/%eps ) then
// Scale up a, b
a = a * Be
b = b * Be
S = S/Be
end
if ( CD < UN*B/%eps ) then
// Scale up c, d
c = c * Be
d = d * Be
// The paper has a mistake there
S = S*Be
end
// Now a, b, c, d in [UN*B/%eps,OV/16]
if ( abs(c) > abs(d) ) then
r = d/c
t = 1/(c+d*r)
qre = (a+b*r)*t
qim = (b-a*r)*t
else
r = c/d
t = 1/(d+c*r)
qre = (b+a*r)*t
qim = (-a+b*r)*t
end
// Scale back
qre = qre * S
qim = qim * S
q = complex(qre,qim)
endfunction
en_US/scripts/cdiv_smith.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
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
// Copyright (C) 2010 - Michael Baudin
//
// smith --
// Returns the complex division r = x / y by Smith's method
//
function r = cdiv_smith ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
r = d/c
den = c + d * r
e = (a + b * r) / den
f = (b - a * r) / den
else
r = c/d
den = c * r + d
e = (a * r + b) / den
f = (b * r - a) / den
end
r = e + %i * f
endfunction
// Copyright (C) 2010 - 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
//
// smith --
// Returns the complex division r = x / y by Smith's method
//
function r = cdiv_smith ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
r = d/c
den = c + d * r
e = (a + b * r) / den
f = (b - a * r) / den
else
r = c/d
den = c * r + d
e = (a * r + b) / den
f = (b * r - a) / den
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(e,f)
endfunction
en_US/scripts/complexdivision.sce
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
// Copyright (C) 2009 - Michael Baudin
exec("cdiv_naive.sci");
exec("cdiv_smith.sci");
exec("cdiv_smith2.sci");
ieee(2);
// An easy case.
r = cdiv_naive ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + 2.0 * %i)/(3.0 + 4.0 * %i)
// A more difficult test case, with overflows
r = cdiv_naive ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_naive ( 1.0 + %i * 1.0 , 1.e-307 + %i * 1.e-307 )
(1.0 + %i * 1.0)/(1.e-307 + %i * 1.e-307)
//
// compare --
// Compare the naive complex division and Scilab's division operator
// a,b,c,d : real numbers of the division (a + ib)/(c + id)
// e,f : expected real and imaginary part of the division
//
function compare (x, y, z)
a = real(x)
b = real(x)
c = real(y)
d = real(y)
e = real(z)
f = real(z)
printf("****************\n");
printf(" %10.5e + I * %10.5e \n" , a , b );
printf(" / %10.5e + I * %10.5e \n" , c , d );
printf(" = %10.5e + I * %10.5e \n" , e , f );
// Use naive computation
r = cdiv_naive ( x, y);
printf("Naive --> %10.5e + I * %10.5e\n" , real(r) , imag(r) );
complexerrors ( r , z );
// Use Scilab's division operator.
div = x/y
printf("Scilab --> %10.5e + I * %10.5e\n" , real(div) , imag(div) );
complexerrors ( div , z );
endfunction
// Compute various errors: in module,
// w.r.t. complex magnitude, w.r.t. to real part, w.r.t. to imaginary part
// x : computed complex
// y : expected complex
function complexerrors ( x , y )
err1 = myerror ( x , y );
err2 = myerror ( real(x) , real(y) );
err3 = myerror ( imag(x) , imag(y) );
printf(" e1 = %10.5e, e2 = %10.5e, e3 = %10.5e\n", err1, err2, err3);
endfunction
// Compute a relative error when possible, and an absolute error if not.
// x : computed complex
// y : expected complex
function err = myerror ( x , y )
if ( y == 0 ) then
err = abs(x - y);
else
err = abs(x - y)/abs(y);
end
endfunction
ieee(2);
// Check that naive implementation does not have a bug
compare (1+%i*2,3+%i*4, 11/25+%i* 2/25);
// Check that naive implementation is not robust with respect to overflow
compare (1+%i* 1, 1+%i* 1e307, 1e-307+%i* -1e-307);
// Check that naive implementation is not robust with respect to underflow
compare (1+%i* 1, 1e-307+%i* 1e-307, 1e307+%i* 0.);
//
// rootsdemo.sf --
// A root-finder algorithm based on the eigenvalues of the
// companion matrix of the polynomial.
//
function result = rootsdemo(p)
A=companion(p);
result = spec(A)
endfunction
// An easy case.
r = cdiv_smith ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + %i * 2.0 )/(3.0 + %i * 4.0)
// A more difficult test case, with overflows
r = cdiv_smith ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_smith ( 1.0 + %i * 1.0 , 1.e-307 + %i * 1.e-307 )
(1.0 + %i * 1.0)/(1.e-307 + %i * 1.e-307)
// A difficult case for Smith
r = cdiv_naive ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
r = cdiv_smith ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
(1.e307 + %i * 1.e-307)/(1.e204 + %i * 1.e-204)
//////////////////////////////////////////////////////
// Test smith2, with ideas from Stewart
// An easy case.
r = cdiv_smith2 ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + %i * 2.0)/(3.0 + %i * 4.0)
// A more difficult test case, with overflows
r = cdiv_smith2 ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_smith2 ( 1.0 + %i * 1.0 , 1.e-308 + %i * 1.e-308 ) // Wrong !!! Fix it !
(1.0 + %i * 1.0)/(1.e-308 + %i * 1.e-308)
// A difficult case for Smith
r = cdiv_smith2 ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
(1.e307 + %i * 1.e-307)/(1.e204 + %i * 1.e-204)
//////////////////////////////////////////////////////
// Find a counter-example for Smith's algorithm
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// There are 22484 such couples.
N = 0
for m = 5:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) ) then
N = N + 1;
//mprintf ( "------------\n")
//mprintf ( "n=%d, m=%d\n", n, m)
//mprintf ( "n-m=%d, n-3m=%d\n", n-m, n-3*m)
//pause
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and n and m are multiples of 20
// There are 59 such couples.
N = 0;
for m = 5:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & modulo(n,20)==0 & modulo(m,20)==0 ) then
N = N + 1;
mprintf ( "n=%d, m=%d, n-m=%d, n-3m=%d\n", n, m,n-m, n-3*m)
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and Smith's method fails.
// There are 2752 such couples.
N = 0
for m = 163:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & (n+m>324) ) then
N = N + 1;
verbose = %t;
if ( verbose & modulo(N,10) == 0 ) then
mprintf ( "N=%d, n=%d, m=%d,n-m=%d, n-3m=%dn n+m =%d\n", N , n, m,n-m, n-3*m , n+m)
end
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and Smith's method fails.
// and n,m are multiples of 20.
// There are 5 such couples.
N = 0
for m = 163:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & (n+m>324) & modulo(n,10)==0 & modulo(m,10)==0 ) then
N = N + 1;
verbose = %t;
if ( verbose ) then
mprintf ( "N=%d, n=%d, m=%d, n-m=%d, n-3m=%dn n+m =%d\n", N, n, m,n-m, n-3*m , n+m)
end
end
end
end
n = [210,220,230,240,250,260,270,280,290,300,240,250,260,270,280,290,300,270,280,290,300,300]
m = [170,170,170,170,170,170,170,170,170,170,180,180,180,180,180,180,180,190,190,190,190,200]
// Copyright (C) 2009 - 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
exec("cdiv_naive.sci");
exec("cdiv_smith.sci");
exec("cdiv_smith2.sci");
ieee(2);
// An easy case.
r = cdiv_naive ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + 2.0 * %i)/(3.0 + 4.0 * %i)
// A more difficult test case, with overflows
r = cdiv_naive ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_naive ( 1.0 + %i * 1.0 , 1.e-307 + %i * 1.e-307 )
(1.0 + %i * 1.0)/(1.e-307 + %i * 1.e-307)
//
// compare --
// Compare the naive complex division and Scilab's division operator
// a,b,c,d : real numbers of the division (a + ib)/(c + id)
// e,f : expected real and imaginary part of the division
//
function compare (x, y, z)
a = real(x)
b = real(x)
c = real(y)
d = real(y)
e = real(z)
f = real(z)
printf("****************\n");
printf(" %10.5e + I * %10.5e \n" , a , b );
printf(" / %10.5e + I * %10.5e \n" , c , d );
printf(" = %10.5e + I * %10.5e \n" , e , f );
// Use naive computation
r = cdiv_naive ( x, y);
printf("Naive --> %10.5e + I * %10.5e\n" , real(r) , imag(r) );
complexerrors ( r , z );
// Use Scilab's division operator.
div = x/y
printf("Scilab --> %10.5e + I * %10.5e\n" , real(div) , imag(div) );
complexerrors ( div , z );
endfunction
// Compute various errors: in module,
// w.r.t. complex magnitude, w.r.t. to real part, w.r.t. to imaginary part
// x : computed complex
// y : expected complex
function complexerrors ( x , y )
err1 = myerror ( x , y );
err2 = myerror ( real(x) , real(y) );
err3 = myerror ( imag(x) , imag(y) );
printf(" e1 = %10.5e, e2 = %10.5e, e3 = %10.5e\n", err1, err2, err3);
endfunction
// Compute a relative error when possible, and an absolute error if not.
// x : computed complex
// y : expected complex
function err = myerror ( x , y )
if ( y == 0 ) then
err = abs(x - y);
else
err = abs(x - y)/abs(y);
end
endfunction
ieee(2);
// Check that naive implementation does not have a bug
compare (1+%i*2,3+%i*4, 11/25+%i* 2/25);
// Check that naive implementation is not robust with respect to overflow
compare (1+%i* 1, 1+%i* 1e307, 1e-307+%i* -1e-307);
// Check that naive implementation is not robust with respect to underflow
compare (1+%i* 1, 1e-307+%i* 1e-307, 1e307+%i* 0.);
//
// rootsdemo.sf --
// A root-finder algorithm based on the eigenvalues of the
// companion matrix of the polynomial.
//
function result = rootsdemo(p)
A=companion(p);
result = spec(A)
endfunction
// An easy case.
r = cdiv_smith ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + %i * 2.0 )/(3.0 + %i * 4.0)
// A more difficult test case, with overflows
r = cdiv_smith ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_smith ( 1.0 + %i * 1.0 , 1.e-307 + %i * 1.e-307 )
(1.0 + %i * 1.0)/(1.e-307 + %i * 1.e-307)
// A difficult case for Smith
r = cdiv_naive ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
r = cdiv_smith ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
(1.e307 + %i * 1.e-307)/(1.e204 + %i * 1.e-204)
//////////////////////////////////////////////////////
// Test smith2, with ideas from Stewart
// An easy case.
r = cdiv_smith2 ( 1.0 + %i * 2.0 , 3.0 + %i * 4.0 )
(1.0 + %i * 2.0)/(3.0 + %i * 4.0)
// A more difficult test case, with overflows
r = cdiv_smith2 ( 1.0 + %i * 1.0 , 1.0 + %i * 1.e307 )
(1.0 + %i * 1.0)/(1.0 + %i * 1.e307)
// A more difficult test case, with underflows
r = cdiv_smith2 ( 1.0 + %i * 1.0 , 1.e-308 + %i * 1.e-308 ) // Wrong !!! Fix it !
(1.0 + %i * 1.0)/(1.e-308 + %i * 1.e-308)
// A difficult case for Smith
r = cdiv_smith2 ( 1.e307 + %i * 1.e-307 , 1.e204 + %i * 1.e-204 )
(1.e307 + %i * 1.e-307)/(1.e204 + %i * 1.e-204)
//////////////////////////////////////////////////////
// Find a counter-example for Smith's algorithm
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// There are 22484 such couples.
N = 0
for m = 5:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) ) then
N = N + 1;
//mprintf ( "------------\n")
//mprintf ( "n=%d, m=%d\n", n, m)
//mprintf ( "n-m=%d, n-3m=%d\n", n-m, n-3*m)
//pause
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and n and m are multiples of 20
// There are 59 such couples.
N = 0;
for m = 5:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & modulo(n,20)==0 & modulo(m,20)==0 ) then
N = N + 1;
mprintf ( "n=%d, m=%d, n-m=%d, n-3m=%d\n", n, m,n-m, n-3*m)
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and Smith's method fails.
// There are 2752 such couples.
N = 0
for m = 163:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & (n+m>324) ) then
N = N + 1;
verbose = %t;
if ( verbose & modulo(N,10) == 0 ) then
mprintf ( "N=%d, n=%d, m=%d,n-m=%d, n-3m=%dn n+m =%d\n", N , n, m,n-m, n-3*m , n+m)
end
end
end
end
// Print all m,n so that the following complex division is exact
// with doubles floating point numbers.
// (10^n + i 10^-n)/(10^m + i 10^-m) = 10^(n-m) + i 10^(n-3m)
// and Smith's method fails.
// and n,m are multiples of 20.
// There are 5 such couples.
N = 0
for m = 163:308
for n = 13:308
if ( (m + 8 < n) & (n-m>=0) & (n-m <308) & (n-3*m >= -307) & (n-3*m<=0) & (n+m>324) & modulo(n,10)==0 & modulo(m,10)==0 ) then
N = N + 1;
verbose = %t;
if ( verbose ) then
mprintf ( "N=%d, n=%d, m=%d, n-m=%d, n-3m=%dn n+m =%d\n", N, n, m,n-m, n-3*m , n+m)
end
end
end
end
n = [210,220,230,240,250,260,270,280,290,300,240,250,260,270,280,290,300,270,280,290,300,300]
m = [170,170,170,170,170,170,170,170,170,170,180,180,180,180,180,180,180,190,190,190,190,200]
en_US/scripts/cdiv_smith2.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
114
115
116
117
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
// Copyright (C) 2010 - Michael Baudin
// An improved Smith method, with ideas from Stewart
function r = cdiv_smith2 ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
z = 1 / c
den = c + d * (d * z)
e = (a + prodminmax ( [ b d z ] ) ) / den
f = (b - prodminmax ( [ a d z ] ) ) / den
else
z = 1 / d
den = c * (c * z) + d
e = ( prodminmax ( [ a c z ] ) + b) / den
f = ( prodminmax ( [ b c z ] ) - a) / den
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(e,f)
endfunction
// A stable floating point product of a vector.
// Based on Stewart's ideas
function p = prodminmax ( v )
n = size(v,"*")
v = matrix(v,n,1)
for i = 1 : n - 1
// Search the maximum absolute value
[ma,kmax]=max(abs(v))
// Search the minimum absolute value in other values
o = [v(1:kmax-1);v(kmax+1:$)]
[mi,komin]=min(abs(o))
// komin is the index in o
// Compute the index kmin in v.
if ( komin<=kmax-1 ) then
kmin = komin
else
kmin = komin+1
end
v(kmin) = v(kmax) * v(kmin)
v(kmax) = []
end
p = v(1)
endfunction
// A stable floating point product of a vector.
// Based on Stewart's ideas
// Does not delete entries in the vector v.
// More clever, but longer: forces to generate the matrices (1:i).
function p = prodminmax2 ( v )
n = size(v,"*")
v = matrix(v,n,1)
for i = n-1 : -1 : 1
// Search the maximum absolute value in v(1:i+1)
[w,k]=max(abs(v(1:i+1)))
// Switch the max and v(i+1)
v([i+1,k]) = v([k,i+1])
// Search the minimum absolute value in v(1:i)
[w,k]=min(abs(v(1:i)))
// Put the product into the min
v(k) = v(k) * v(i+1)
end
p = v(1)
endfunction
if ( %f ) then
p = prodminmax ( [1 3 1 2] ) // 6
p = prodminmax ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodminmax ( [-1 1 1] ) // -1
p = prodminmax ( [-1 -1 1] ) // 1
p = prodminmax ( [-1 -1 -1] ) // -1
p = prodminmax ( [-1 2 3] ) // -6
p = prodminmax ( [-1 -2 3] ) // 6
p = prodminmax ( [-1 2 -3] ) // 6
p = prodminmax ( [1 2 -3] ) // -6
p = prodminmax ( [1.e200 1.e200 1.e-100] ) // 1.+300
for k = 1 : 10
v=round(10*rand(10,1)-5);
v(find(v==0))=[];
p1 = prodminmax ( v );
p2 = prod ( v );
disp([k p1 p2])
end
//
p = prodminmax2 ( [1 3 1 2] ) // 6
p = prodminmax2 ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodminmax2 ( [-1 1 1] ) // -1
p = prodminmax2 ( [-1 -1 1] ) // 1
p = prodminmax2 ( [-1 -1 -1] ) // -1
p = prodminmax2 ( [-1 2 3] ) // -6
p = prodminmax2 ( [-1 -2 3] ) // 6
p = prodminmax2 ( [-1 2 -3] ) // 6
p = prodminmax2 ( [1 2 -3] ) // -6
p = prodminmax2 ( [1.e200 1.e200 1.e-100] ) // 1.+300
for k = 1 : 10
v=round(10*rand(10,1)-5);
v(find(v==0))=[];
p1 = prodminmax2 ( v );
p2 = prod ( v );
disp([k p1 p2])
end
end
if ( %f ) then
//
// Test perfs
exec("benchfun.sci")
n=5000;
v=round(10*rand(n,1)-5);
benchfun ( "prodminmax " , prodminmax , list(v) , 1 , 10 );
benchfun ( "prodminmax2" , prodminmax2 , list(v) , 1 , 10 );
benchfun ( "prod " , prod , list(v) , 1 , 10 );
end
// Copyright (C) 2010 - 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
// An improved Smith method, with ideas from Stewart
function r = cdiv_smith2 ( x,y )
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
z = 1 / c
den = c + d * (d * z)
e = (a + prodminmax ( [ b d z ] ) ) / den
f = (b - prodminmax ( [ a d z ] ) ) / den
else
z = 1 / d
den = c * (c * z) + d
e = ( prodminmax ( [ a c z ] ) + b) / den
f = ( prodminmax ( [ b c z ] ) - a) / den
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(e,f)
endfunction
// A stable floating point product of a vector.
// Based on Stewart's ideas
function p = prodminmax ( v )
n = size(v,"*")
v = matrix(v,n,1)
for i = 1 : n - 1
// Search the maximum absolute value
[ma,kmax]=max(abs(v))
// Search the minimum absolute value in other values
o = [v(1:kmax-1);v(kmax+1:$)]
[mi,komin]=min(abs(o))
// komin is the index in o
// Compute the index kmin in v.
if ( komin<=kmax-1 ) then
kmin = komin
else
kmin = komin+1
end
v(kmin) = v(kmax) * v(kmin)
v(kmax) = []
end
p = v(1)
endfunction
// A stable floating point product of a vector.
// Based on Stewart's ideas
// Does not delete entries in the vector v.
// More clever, but longer: forces to generate the matrices (1:i).
function p = prodminmax2 ( v )
n = size(v,"*")
v = matrix(v,n,1)
for i = n-1 : -1 : 1
// Search the maximum absolute value in v(1:i+1)
[w,k]=max(abs(v(1:i+1)))
// Switch the max and v(i+1)
v([i+1,k]) = v([k,i+1])
// Search the minimum absolute value in v(1:i)
[w,k]=min(abs(v(1:i)))
// Put the product into the min
v(k) = v(k) * v(i+1)
end
p = v(1)
endfunction
if ( %f ) then
p = prodminmax ( [1 3 1 2] ) // 6
p = prodminmax ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodminmax ( [-1 1 1] ) // -1
p = prodminmax ( [-1 -1 1] ) // 1
p = prodminmax ( [-1 -1 -1] ) // -1
p = prodminmax ( [-1 2 3] ) // -6
p = prodminmax ( [-1 -2 3] ) // 6
p = prodminmax ( [-1 2 -3] ) // 6
p = prodminmax ( [1 2 -3] ) // -6
p = prodminmax ( [1.e200 1.e200 1.e-100] ) // 1.+300
for k = 1 : 10
v=round(10*rand(10,1)-5);
v(find(v==0))=[];
p1 = prodminmax ( v );
p2 = prod ( v );
disp([k p1 p2])
end
//
p = prodminmax2 ( [1 3 1 2] ) // 6
p = prodminmax2 ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodminmax2 ( [-1 1 1] ) // -1
p = prodminmax2 ( [-1 -1 1] ) // 1
p = prodminmax2 ( [-1 -1 -1] ) // -1
p = prodminmax2 ( [-1 2 3] ) // -6
p = prodminmax2 ( [-1 -2 3] ) // 6
p = prodminmax2 ( [-1 2 -3] ) // 6
p = prodminmax2 ( [1 2 -3] ) // -6
p = prodminmax2 ( [1.e200 1.e200 1.e-100] ) // 1.+300
for k = 1 : 10
v=round(10*rand(10,1)-5);
v(find(v==0))=[];
p1 = prodminmax2 ( v );
p2 = prod ( v );
disp([k p1 p2])
end
end
if ( %f ) then
//
// Test perfs
exec("benchfun.sci")
n=5000;
v=round(10*rand(n,1)-5);
benchfun ( "prodminmax " , prodminmax , list(v) , 1 , 10 );
benchfun ( "prodminmax2" , prodminmax2 , list(v) , 1 , 10 );
benchfun ( "prod " , prod , list(v) , 1 , 10 );
end
en_US/scripts/cdiv_smith3.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
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
// Copyright (C) 2010 - Michael Baudin
// An improved Smith method, with an exp-log product
function r = cdiv_smith3 (x,y)
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
z = 1 / c
den = c + d * (d * z)
e = (a + prodexplog ( [ b d z ] ) ) / den
f = (b - prodexplog ( [ a d z ] ) ) / den
else
z = 1 / d
den = c * (c * z) + d
e = ( prodexplog ( [ a c z ] ) + b) / den
f = ( prodexplog ( [ b c z ] ) - a) / den
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(e,f)
endfunction
// A product of a vector.
// Based on log(exp) trick : cost much more !
function p = prodexplog ( v )
pe = exp(sum(log(v)))
p = abs(pe) * sign(real(pe))
endfunction
if ( %f ) then
p = prodexplog ( [1 3 1 2] ) // 6
p = prodexplog ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodexplog ( [-1 1 1] ) // -1
p = prodexplog ( [1 -1 1] ) // -1
p = prodexplog ( [1.e200 1.e200 1.e-100] ) // 1.+300
end
// Copyright (C) 2010 - 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
// An improved Smith method, with an exp-log product
function r = cdiv_smith3 (x,y)
a = real(x)
b = imag(x)
c = real(y)
d = imag(y)
if ( abs(d) <= abs(c) ) then
z = 1 / c
den = c + d * (d * z)
e = (a + prodexplog ( [ b d z ] ) ) / den
f = (b - prodexplog ( [ a d z ] ) ) / den
else
z = 1 / d
den = c * (c * z) + d
e = ( prodexplog ( [ a c z ] ) + b) / den
f = ( prodexplog ( [ b c z ] ) - a) / den
end
// Avoid using x + %i * y, which may create
// %inf+%i*%inf, which actually multiplies 0 and %inf,
// generating an unwanted %nan.
r = complex(e,f)
endfunction
// A product of a vector.
// Based on log(exp) trick : cost much more !
function p = prodexplog ( v )
pe = exp(sum(log(v)))
p = abs(pe) * sign(real(pe))
endfunction
if ( %f ) then
p = prodexplog ( [1 3 1 2] ) // 6
p = prodexplog ( [-1 2 -3 4 1 -2 3 4] ) // -576
p = prodexplog ( [-1 1 1] ) // -1
p = prodexplog ( [1 -1 1] ) // -1
p = prodexplog ( [1.e200 1.e200 1.e-100] ) // 1.+300
end
en_US/scripts/cdiv_priest.c
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
// Copyright (C) 2010 - 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
// References
//
// "Efficient scaling for complex division"
// Douglas M. Priest
// Sun Microsystems
// ACM Transactions on Mathematical software, Vol. 30, No. 4, December 2004
/*
* Set v[0] + i v[1] := (z[0] + i z[1]) / (w[0] + i w[1]) assuming:
*
* 1. z[0], z[1], w[0], w[1] are all finite, and w[0], w[1] are not
* both zero
*
* 2. "int" refers to a 32-bit integer type, "long long int" refers
* to a 64-bit integer type, and "double" refers to an IEEE 754-
* compliant 64-bit floating point type
*/
void cdiv(double *v, const double *z, const double *w)
{
union {
long long int i;
double d;
} aa, bb, cc, dd, ss;
double a, b, c, d, t;
int ha, hb, hc, hd, hz, hw, hs;
/* name the components of z and w */
a = z[0];
b = z[1];
c = w[0];
d = w[1];
/* extract high-order 32 bits to estimate |z| and |w| */
aa.d = a;
bb.d = b;
ha = (aa.i >> 32) & 0x7fffffff;
hb = (bb.i >> 32) & 0x7fffffff;
hz = (ha > hb)? ha : hb;
cc.d = c;
dd.d = d;
hc = (cc.i >> 32) & 0x7fffffff;
hd = (dd.i >> 32) & 0x7fffffff;
hw = (hc > hd)? hc : hd;
/* compute the scale factor */
if (hz < 0x07200000 && hw >= 0x32800000 && hw < 0x47100000) {
/* |z| < 2^-909 and 2^-215 <= |w| < 2^114 */
hs = (((0x47100000 - hw) >> 1) & 0xfff00000) + 0x3ff00000;
}
else
{
hs = (((hw >> 2) - hw) + 0x6fd7ffff) & 0xfff00000;
}
ss.i = (long long int) hs << 32;
/* scale c and d, and compute the quotient */
c *= ss.d;
d *= ss.d;
t = 1.0 / (c * c + d * d);
c *= ss.d;
d *= ss.d;
v[0] = (a * c + b * d) * t;
v[1] = (b * c - a * d) * t;
}
en_US/scripts/cdiv_polar.sci
11
2
3
4
5
6
7
28
39
410
// Copyright (C) 2010 - 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
// A complex division based on polar form.

Archive Download the corresponding diff file

Revision: 34