Module: Numo::GSL::Poly

Defined in:
ext/numo/gsl/poly/gsl_poly.c

Defined Under Namespace

Classes: ComplexWorkspace

Class Method Summary collapse

Class Method Details

.complex_solve(a) ⇒ Numo::DComplex

This function computes the roots of the general polynomial $P(x) = a_0 + a_1 x + a_2 x^2 + … + a_n-1 x^n-1$ P(x) = a_0 + a_1 x + a_2 x^2 + … + a_[n-1] x^[n-1] using balanced-QR reduction of the companion matrix. The parameter n specifies the length of the coefficient array. The coefficient of the highest order term must be non-zero. The function requires a workspace w of the appropriate size. The n-1 roots are returned in the packed complex array z of length 2(n-1), alternating real and imaginary parts.

The function returns GSL_SUCCESS if all the roots are found. If the QR reduction does not converge, the error handler is invoked with an error code of GSL_EFAILED. Note that due to finite precision, roots of higher multiplicity are returned as a cluster of simple roots with reduced accuracy. The solution of polynomials with higher-order roots requires specialized algorithms that take the multiplicity structure into account (see e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp 218–236).

Parameters:

  • a (Numo::DFloat)

Returns:

  • (Numo::DComplex) —

    z result



564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 564

static VALUE
poly_s_complex_solve(VALUE mod, VALUE v1)
{
    size_t shape[0];
    ndfunc_arg_in_t ain[1] = {{cDF,1}};
    ndfunc_arg_out_t aout[1] = {{cDC,1,shape}};
    ndfunc_t ndf = {iter_poly_s_complex_solve, NO_LOOP, 1,1, ain,aout};
    gsl_poly_complex_workspace *w;
    void     *opts[1];
    VALUE     vz, vws;
    size_t    n;
    narray_t *na;

    v1 = rb_funcall(cDF, rb_intern("cast"), 1, v1);
    GetNArray(v1,na);
    if (na->ndim == 0) {
        rb_raise(nary_eDimensionError,"ndim(=%d) should >= %d", na->ndim, 0);
    }
    if (na->shape[na->ndim-1] < 2) {
        rb_raise(nary_eShapeError, "last axis size must be >= 2");
    }
    n = na->shape[na->ndim-1];
    shape[0] = n-1;

    vws = poly_complex_workspace_s_new(cComplexWorkspace, SIZET2NUM(n));
    TypedData_Get_Struct(vws, gsl_poly_complex_workspace, &poly_complex_workspace_data_type, w);
    opts[0] = w;

    vz = na_ndloop3(&ndf, opts, 1, v1);
    RB_GC_GUARD(vws);
    RB_GC_GUARD(v1);
    return vz;
}

.complex_solve_cubic(a, b, c) ⇒ [Numo::DComplex, Numo::DComplex, Numo::DComplex, Numo::Int]

This function finds the complex roots of the cubic equation,

z^3 + a z^2 + b z + c = 0

The number of complex roots is returned (always three) and the locations of the roots are stored in z0, z1 and z2. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components.

Parameters:

  • a (Numo::DFloat)
  • b (Numo::DFloat)
  • c (Numo::DFloat)

Returns:

  • ([Numo::DComplex, Numo::DComplex, Numo::DComplex, Numo::Int]) —

    array of [z0, z1, z2, return]



503
504
505
506
507
508
509
510
511
512
513
514
515
516
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 503

static VALUE
poly_s_complex_solve_cubic(VALUE mod,VALUE v0,VALUE v1,VALUE v2)
{
    
    

    ndfunc_arg_in_t ain[3] = {{cDF,0},{cDF,0},{cDF,0}};
    ndfunc_arg_out_t aout[4] = {{cDC,0},{cDC,0},{cDC,0},{cI,0}};
    ndfunc_t ndf = {iter_poly_s_complex_solve_cubic,NO_LOOP|NDF_INPLACE|NDF_EXTRACT,
                    3,4,ain,aout};
    
    
    return na_ndloop(&ndf,3,v0,v1,v2); 
}

.complex_solve_quadratic(a, b, c) ⇒ [Numo::DComplex, Numo::DComplex, Numo::Int]

This function finds the complex roots of the quadratic equation,

a z^2 + b z + c = 0

The number of complex roots is returned (either one or two) and the locations of the roots are stored in z0 and z1. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components. If only one real root is found (i.e. if a=0) then it is stored in z0.

Parameters:

  • a (Numo::DFloat)
  • b (Numo::DFloat)
  • c (Numo::DFloat)

Returns:

  • ([Numo::DComplex, Numo::DComplex, Numo::Int]) —

    array of [z0, z1, return]



366
367
368
369
370
371
372
373
374
375
376
377
378
379
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 366

static VALUE
poly_s_complex_solve_quadratic(VALUE mod,VALUE v0,VALUE v1,VALUE v2)
{
    
    

    ndfunc_arg_in_t ain[3] = {{cDF,0},{cDF,0},{cDF,0}};
    ndfunc_arg_out_t aout[3] = {{cDC,0},{cDC,0},{cI,0}};
    ndfunc_t ndf = {iter_poly_s_complex_solve_quadratic,NO_LOOP|NDF_INPLACE|NDF_EXTRACT,
                    3,3,ain,aout};
    
    
    return na_ndloop(&ndf,3,v0,v1,v2); 
}

.eval(c, x) ⇒ Numo::DFloat or DComplex

This function calls gsl_poly_eval or gsl_poly_complex_eval or gsl_complex_poly_complex_eval according to whether argument is complex or not.

This function evaluates a polynomial with real coefficients for the real variable x.

Parameters:

  • c (Numo::DFloat or DComplex)
  • x (Numo::DFloat or DComplex)

Returns:

  • (Numo::DFloat or DComplex) 


185
186
187
188
189
190
191
192
193
194
195
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 185

static VALUE
poly_s_eval(VALUE mod, VALUE v0, VALUE v1)
{
    if (is_complex_nary(v0)) {
        return poly_eval_cc(v0,v1);
    } else if (is_complex_nary(v1)) {
        return poly_eval_fc(v0,v1);
    } else {
        return poly_eval_ff(v0,v1);
    }
}

.eval_derivs(c, x, lenres) ⇒ Numo::DFloat

This function evaluates a polynomial and its derivatives storing the results in the array res of size lenres. The output array contains the values of d^k P/d x^k for the specified value of x starting with k = 0.

Parameters:

  • c (Numo::DFloat)
  • x (Numo::DFloat)
  • lenres (Integer)

Returns:

  • (Numo::DFloat) 


231
232
233
234
235
236
237
238
239
240
241
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 231

static VALUE
poly_s_eval_derivs(VALUE mod, VALUE v0, VALUE v1, VALUE v2)
{
    size_t shape[1];
    ndfunc_arg_in_t ain[2] = {{cDF,1},{cDF,0}};
    ndfunc_arg_out_t aout[1] = {{cDF,1,shape}};
    ndfunc_t ndf = {iter_poly_s_eval_derivs,NO_LOOP|NDF_INPLACE|NDF_EXTRACT,2,1,ain,aout};

    shape[0] = NUM2SIZET(v2);
    return na_ndloop(&ndf,2,v0,v1);
}

.solve_cubic(a, b, c) ⇒ [Numo::DFloat, Numo::DFloat, Numo::DFloat, Numo::Int]

This function finds the real roots of the cubic equation,

x^3 + a x^2 + b x + c = 0

with a leading coefficient of unity. The number of real roots (either one or three) is returned, and their locations are stored in x0, x1 and x2. If one real root is found then only x0 is modified. When three real roots are found they are stored in x0, x1 and x2 in ascending order. The case of coincident roots is not considered special. For example, the equation (x-1)^3=0 will have three roots with exactly equal values. As in the quadratic case, finite precision may cause equal or closely-spaced real roots to move off the real axis into the complex plane, leading to a discrete change in the number of real roots.

Parameters:

  • a (Numo::DFloat)
  • b (Numo::DFloat)
  • c (Numo::DFloat)

Returns:

  • ([Numo::DFloat, Numo::DFloat, Numo::DFloat, Numo::Int]) —

    array of [x0, x1, x2, return]



437
438
439
440
441
442
443
444
445
446
447
448
449
450
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 437

static VALUE
poly_s_solve_cubic(VALUE mod,VALUE v0,VALUE v1,VALUE v2)
{
    
    

    ndfunc_arg_in_t ain[3] = {{cDF,0},{cDF,0},{cDF,0}};
    ndfunc_arg_out_t aout[4] = {{cDF,0},{cDF,0},{cDF,0},{cI,0}};
    ndfunc_t ndf = {iter_poly_s_solve_cubic,NO_LOOP|NDF_INPLACE|NDF_EXTRACT,
                    3,4,ain,aout};
    
    
    return na_ndloop(&ndf,3,v0,v1,v2); 
}

.solve_quadratic(a, b, c) ⇒ [Numo::DFloat, Numo::DFloat, Numo::Int]

This function finds the real roots of the quadratic equation,

a x^2 + b x + c = 0

The number of real roots (either zero, one or two) is returned, and their locations are stored in x0 and x1. If no real roots are found then x0 and x1 are not modified. If one real root is found (i.e. if a=0) then it is stored in x0. When two real roots are found they are stored in x0 and x1 in ascending order. The case of coincident roots is not considered special. For example (x-1)^2=0 will have two roots, which happen to have exactly equal values.

The number of roots found depends on the sign of the discriminant b^2 - 4 a c. This will be subject to rounding and cancellation errors when computed in double precision, and will also be subject to errors if the coefficients of the polynomial are inexact. These errors may cause a discrete change in the number of roots. However, for polynomials with small integer coefficients the discriminant can always be computed exactly.

Parameters:

  • a (Numo::DFloat)
  • b (Numo::DFloat)
  • c (Numo::DFloat)

Returns:

  • ([Numo::DFloat, Numo::DFloat, Numo::Int]) —

    array of [x0, x1, return]



302
303
304
305
306
307
308
309
310
311
312
313
314
315
 
# File 'ext/numo/gsl/poly/gsl_poly.c', line 302

static VALUE
poly_s_solve_quadratic(VALUE mod,VALUE v0,VALUE v1,VALUE v2)
{
    
    

    ndfunc_arg_in_t ain[3] = {{cDF,0},{cDF,0},{cDF,0}};
    ndfunc_arg_out_t aout[3] = {{cDF,0},{cDF,0},{cI,0}};
    ndfunc_t ndf = {iter_poly_s_solve_quadratic,NO_LOOP|NDF_INPLACE|NDF_EXTRACT,
                    3,3,ain,aout};
    
    
    return na_ndloop(&ndf,3,v0,v1,v2); 
}