91 for(k = 1; k < 4; k++ ) { p[
k] = p[
k]/p[0]; }
96 b = 0.5*( x*( t/1.5 - p[2] ) + p[3] );
103 d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
107 if( b > 0. ) { b = -
d; }
111 d = std::sqrt(0.75)*(b -
c);
117 if( ( b > 0. && x <= 0. ) || ( b < 0. && x > 0. ) )
134 if( b == 0. ) { d = std::atan(1.0)/1.5; }
135 else { d = std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
137 if( b < 0. ) { b = std::sqrt(t)*2.0; }
138 else { b = -2.0*std::sqrt(t); }
141 t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*
c;
146 if( std::fabs(c) > std::fabs(t) ) { r[1][3] =
c; }
152 if( std::fabs(d) > std::fabs(t) ) { r[1][2] =
d; }
160 for(k = 1; k < 4; k++ ) { r[2][
k] = 0.; }
178 for( k = 1; k < 5; k++) { p[
k] = p[
k]/p[0]; }
185 b = p[3] + b*( c - p[2] );
187 c = p[4] + e*( e*a - p[3] );
191 p[2] = (p[1]*p[1]-
c)*0.25;
198 for( k = 1; k < 4; k++ )
200 if( r[2][k] == 0. && r[1][k] > 0 )
205 if ( a >= 0. && b >= 0.) { p[1] = std::sqrt(d); }
206 else if( a <= 0. && b <= 0.) { p[1] = std::sqrt(d); }
207 else { p[1] = -std::sqrt(d); }
209 b = 0.5*( a + b/p[1] );
214 for( i = 1; i < 3; i++ )
216 for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
222 for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] -
e; }
234 if( d > 0. ) { p[1] = std::sqrt(d); }
238 if( p[1] > 0.) { b = std::sqrt(p[2])*2.0 + p[1]; }
239 else { b = -std::sqrt(p[2])*2.0 + p[1]; }
241 if( b != 0.) { p[1] = 0; }
244 for(k = 1; k < 5; k++ )
256 for( k = 1; k < 3; k++ )
258 for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][
k]; }
264 for( k = 1; k < 5; k++ ) { r[1][
k] = r[1][
k] -
e; }
280 for( k = 0; k < 4; k++ ) { reRoot[
k] =
DBL_MAX; }
284 for( k = 1; k < 5; k++) { p[
k] = p[
k]/p[0]; }
296 p[3] = 4*a2*a0 - a1*a1 - a3*a3*
a0;
300 for( k = 1; k < 4; k++ )
310 for( k = 1; k < 4; k++ )
312 if ( reRoot[k] < y1 ) { y1 = reRoot[
k]; }
315 R2 = 0.25*a3*a3 - a2 +
y1;
316 b = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
317 c = 0.75*a3*a3 - 2*a2;
330 r[1][1] = -0.25*a3 + 0.5*R + 0.5*
D;
331 r[1][2] = -0.25*a3 + 0.5*R - 0.5*
D;
338 r[1][1] = -0.25*a3 + 0.5*
R;
339 r[1][2] = -0.25*a3 + 0.5*
R;
346 r[1][3] = -0.25*a3 - 0.5*R + 0.5*
E;
347 r[1][4] = -0.25*a3 - 0.5*R - 0.5*
E;
354 r[1][3] = -0.25*a3 - 0.5*
R;
355 r[1][4] = -0.25*a3 - 0.5*
R;
366 r[1][1] = -0.25*a3 + 0.5*real(CD);
367 r[1][2] = -0.25*a3 - 0.5*real(CD);
368 r[2][1] = 0.5*R + 0.5*imag(CD);
369 r[2][2] = 0.5*R - 0.5*imag(CD);
373 r[1][3] = -0.25*a3 + 0.5*real(CE);
374 r[1][4] = -0.25*a3 - 0.5*real(CE);
375 r[2][3] = -0.5*R + 0.5*imag(CE);
376 r[2][4] = -0.5*R - 0.5*imag(CE);
382 D2 = c + std::sqrt(d);
383 E2 = c - std::sqrt(d);
388 r[1][1] = -0.25*a3 + 0.5*R + 0.5*
D;
389 r[1][2] = -0.25*a3 + 0.5*R - 0.5*
D;
396 r[1][1] = -0.25*a3 + 0.5*
R;
397 r[1][2] = -0.25*a3 + 0.5*
R;
404 r[1][3] = -0.25*a3 - 0.5*R + 0.5*
E;
405 r[1][4] = -0.25*a3 - 0.5*R - 0.5*
E;
412 r[1][3] = -0.25*a3 - 0.5*
R;
413 r[1][4] = -0.25*a3 - 0.5*
R;
424 r[1][1] = -0.25*a3 + 0.5*real(CD);
425 r[1][2] = -0.25*a3 - 0.5*real(CD);
426 r[2][1] = 0.5*R + 0.5*imag(CD);
427 r[2][2] = 0.5*R - 0.5*imag(CD);
432 r[1][3] = -0.25*a3 + 0.5*real(CE);
433 r[1][4] = -0.25*a3 - 0.5*real(CE);
434 r[2][3] = -0.5*R + 0.5*imag(CE);
435 r[2][4] = -0.5*R - 0.5*imag(CE);