-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathroots.c
More file actions
101 lines (93 loc) · 1.86 KB
/
roots.c
File metadata and controls
101 lines (93 loc) · 1.86 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
92
93
94
95
96
97
98
99
100
101
#include "roots.h"
/*
* This function calculates 4 complex roots for the passed value.
*
* dst - storage for the found roots,
* the roots are located for the passed i-th value each as follows:
* dst[i][0] - 1st
* dst[i][1] - 2nd
* dst[i][2] - 3d
* dst[i][3] - 4th
*
* src - array of the double values to calculate the roots for
*
* len - the 'src' length
*
* */
int get4Roots(struct Z dst[][4], const double *src, int len)
{
if (len <= 0)
return -1;
struct Z *r = *dst; /* roots */
for(; len--; src++)
{
double val = *src;
if NEIGHBOR(val, 0, EPS_EQ)
{
for(int i = 0; i < 4; i++, r++)
{
r->Re = r->Im = 0;
}
}
else
{
double x = 10.0;
if (val < 0) val = -val;
double nx;
/*
* 4th root
for (;;)
{
nx = (3*x + val / (x*x*x)) / 4;
if NEIGHBOR(x, nx, EPS_STOP)
break;
x = nx;
}
*/
/*
* 2nd root 2 times
*/
for (;;)
{
nx = (x + val / x) / 2;
if NEIGHBOR(x, nx, EPS_STOP)
break;
x = nx;
}
for (val = x;;)
{
nx = (x + val / x) / 2;
if NEIGHBOR(x, nx, EPS_STOP)
break;
x = nx;
}
if (*src > 0)
{
r->Re = x; r->Im = 0;
r++;
r->Re = 0; r->Im = x;
r++;
r->Re = -x; r->Im = 0;
r++;
r->Re = 0; r->Im = -x;
r++;
}
else
{
r->Re = x * SIN_PI_4;
r->Im = x * SIN_PI_4;
r++;
r->Re = -x * SIN_PI_4;
r->Im = x * SIN_PI_4;
r++;
r->Re = -x * SIN_PI_4;
r->Im = -x * SIN_PI_4;
r++;
r->Re = x * SIN_PI_4;
r->Im = -x * SIN_PI_4;
r++;
}
}
}
return 0;
}