Subversion Repositories forwardcom

/*******************************  integrate.as  *******************************

* Author:        Agner Fog

* date created:  2018-03-30

* Last modified: 2021-05-25

* Version:       1.11

* Project:       ForwardCom library math.li

* Description:   Numerical integration of a function f(x) over a given x-interval

*                Uses 4-point Gauss-Legendre integration method

*

* Parameters:

* v0 = start of x interval, double scalar

* v1 = end   of x interval, double scalar

* r0 = function pointer

* r1 = number of function points, int32

* r2 = method. currently unused

* r3 = optional extra parameter to the function, can be integer or pointer

* return value: double scalar in v0

* will return NAN if the maximum vector length is less than 256 bits or if the function returns NAN

*

* The function that is integrated is pointed to by r0. This function must take a double vector

* as input and return a double vector as results. It can have one optional extra parameter which

* can be an integer or a pointer. The function must conform to the default register use.

*

* The number of function points will be rounded up to a power of 2 or a multiple of the maximum

* vector length. It must be at least 4.

*

* Copyright 2018-2021 GNU General Public License http://www.gnu.org/licenses

NOTE: not debugged after conversion to version 1.11 !!

to do: use masked replace rather then repeat_block to broadcast constants

*****************************************************************************/

// define constants

%xval1   = 0.861136311594                        // x-values for 4-point Gauss-Legendre integration

%xval2   = 0.339981043585

%weight1 = 0.347854845137                        // weights  for 4-point Gauss-Legendre integration

%weight2 = 0.652145154863

%nan     = 0x7FF8100000000000                    // NAN error return

const section read ip                            // read-only data section

xvals:   double (1-xval1)/2, (1-xval2)/2, (1+xval2)/2, (1+xval1)/2  // x-points in interval from 0 - 1

weights: double weight1/2, weight2/2, weight2/2, weight1/2          // corresponding weights

const end

code section execute align = 4                   // code section

_integrate function public

push (r16, 18)                                   // save registers

push (v16, 19)                                   // save registers

int64  r16 = r0                                  // function pointer

int64  r18 = r3                                  // save optional extra parameter to function

if (int32 r1 < 4) {

   int32+ r1 = 4                                 // must have at least 4 points

}

double v2 = broadcast_max(0)                     // find maximum vector length

int64 r3 = get_num(v2)                           // maximum elements

if (int32 r1 < r3) {                             // number fits into a single vector

   int32  r1 = roundp2(r1, 1)                    // round up to nearest power of 2

}

else {                                           // more than one vector needed

   int32 r1 = r1 + r3 - 1                        // round up to nearest multiple of maximum vector length

   int32 r3 = -r3

   int32+ r1 &= r3

}

uint32+ r3 = r1 >> 2                             // number of steps, nsteps = points / 4

int64  r0 = 4*8                                  // length of one step

int64  r2 = r3 * r0                              // total length, bytes

double v16 = [xvals,   length = r0]              // x values for one step

double v17 = [weights, length = r0]              // weights for one step

int64 r4 = get_len(v16)

if (int64 r4 < r0) {                             // maximum vector length is insufficient for one step. a fix for this is not implemented

   int64 v0 = nan                                // return NAN to indicate error

   jump ERROREXIT

}

double v18 = v1 - v0                             // length of x interval

int64  v2  = gp2vec(r3)                          // number of steps

double v2  = int2float(v2, 0)                    // same, as float

double v18 /= v2                                 // x-step size = length / nsteps

double v18 = broad(v18, r2)                      // broadcast x-step size

double v16 *= v18                                // first 4 x-values relative to start

double v0  = broad(v0, r2)                       // broadcast x start value

double v16 += v0                                 // first 4 x-values

double v17 *= v18                                // multiply weights by step size

double v16 = repeat_block(v16, r2, 4*8)          // repeat x-values block nsteps times or until the maximum vector length

double v17 = repeat_block(v17, r2, 4*8)          // repeat weights  block nsteps times or until the maximum vector length

int64  v2  = make_sequence(r1, 0)                // 0, 1, 2, ..

int64  v2  >>= 2                                 // 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, ...

double v2 = int2float(v2, 0)                     // same, as double

double v2 *= v18                                 // multiply by x-step size

double v16 += v2                                 // make each x block = previous block + x-step size

int64  r3 = get_len(v2)                          // actual vector length

uint32 r17 = r2 / r3                             // number of iterations

int64 r2 = get_num(v2)                           // vector length

int64 r2--

double v2 = extract(v2, r2)                      // get last addend of last block, broadcast

//double v2 = rotate_up(r3, v2)

double v2 += v18                                 // start value for next x vector

//double v19 = broad(r3, v2)

double v19 = v2 + v18                            // (x-step size) * (steps per vector). add this to x vector to get next x vector

double v18 = replace(v16, 0)                     // vector of zeroes for summation

// calculation loop. calculate as many function values per iteration as the maximum vector length allows

while (int32+ r17 > 0) {

   double v0 = v16                               // x-values

   int64  r0 = r18                               // optional extra parameter

   call (r16)                                    // function(x)

   double v0  *= v17                             // multiply results by weights

   double v18 += v0                              // summation

   double v16 += v19                             // next vector of x-values

   int32+ r17--                                  // loop counter

}

// the sums are in v18. make horizontal sum

int64  r2 = get_len(v18)                         // vector length

while (uint64 r2 > 8) {                          // loop to calculate horizontal sum

   uint64 r2 >>= 1                               // the vector length is halved

   double v1 = shift_reduce(v18, r2)             // get upper half of vector

   double v18 = v1 + v18                         // Add upper half and lower half

}

double v0 = v18                                  // return value

ERROREXIT:

pop (v16, 19)                                    // restore registers

pop (r16, 18)                                    // restore registers

return

_integrate end

code end