Algorithms

const GROUPS_COUNTER = 10;
const MIN_ASCII = 65;
const MAX_ASCII = 69;

/*
 * generate children pairs.
 */
function generateGroups($n) {
    $matrix = [];
    $exist = [];
    echo "************************\n";

    for ($i = 0; $i < $n; $i++) {
        $a = chr(rand(MIN_ASCII, MAX_ASCII));
        $b = chr(rand(MIN_ASCII, MAX_ASCII));

        if ($a == $b || isset($exist[$a][$b]) || isset($exist[$b][$a])) {
            $i--;
            continue;
        }

        $matrix[$i][0] = $a;
        $matrix[$i][1] = $b;
        $exist[$a][$b] = true;

        echo $a . " " . $b . "\n";
    }

    echo "************************\n";

    return $matrix;
}

/*
 * find all children which are not like by parent.
 */
function findGroups($groups) {
    $dislikeList = [];
    $checked = [];

    foreach ($groups as $groupParent) {
        if (isset($checked[$groupParent[0]])) {
            continue;
        }

        $groupChildList = [];

        foreach ($groups as $groupChild) {
            if ($groupParent[0] == $groupChild[0]) {
                array_push($groupChildList, $groupChild[1]);
            }
        }

        echo $groupParent[0] . " => " . implode(", ", $groupChildList) . "\n";
        $dislikeList[$groupParent[0]] = $groupChildList;
        $checked[$groupParent[0]] = true;
    }

    echo "************************\n";

    return $dislikeList;
}

/*
 * phpTarjanEntry() iterates through the graph array rows, executing phpTarjan().
 */
function phpTarjanEntry($groups) {
    global $cycles;
    global $marked;
    global $markedStack;
    global $pointStack;

    $cycles = array();
    $marked = array();
    $markedStack = array();
    $pointStack = array();

    $marked = array();

    foreach ($groups as $key => $value) {
        $marked[$key] = FALSE;
    }

    foreach ($groups as $key => $value) {
        phpTarjan($key, $key);

        while (!empty($markedStack)) {
            $marked[array_pop($markedStack)] = FALSE;
        }
    }

    $cycles = array_keys($cycles);

    return $cycles;
}

/*
 * Recursive function to detect strongly connected components (cycles, loops).
 */
function phpTarjan($s, $v) {
    global $G;
    global $cycles;
    global $marked;
    global $markedStack;
    global $pointStack;

    $f = FALSE;
    $pointStack[] = $v;
    $marked[$v] = TRUE;
    $markedStack[] = $v;

    foreach ($G[$v] as $w) {
        if ($w < $s) {
            $G[$w] = array();
        } else if ($w == $s) {
            $cycles[implode('|', $pointStack)] = TRUE;
            $f = TRUE;
        } else if (isset($marked[$w]) && $marked[$w] === FALSE) {
            $g = phpTarjan($s, $w);

            if (!empty($f) OR ! empty($g)) {
                $f = TRUE;
            }
        }
    }

    if ($f === TRUE) {
        while (end($markedStack) != $v) {
            $marked[array_pop($markedStack)] = FALSE;
        }

        array_pop($markedStack);
        $marked[$v] = FALSE;
    }

    array_pop($pointStack);
    return $f;
}

$mixGroups = generateGroups(GROUPS_COUNTER);
$res = findGroups($mixGroups);
$G = $res;
$cycles = phpTarjanEntry($res);

echo '<p>Cycles found: ' . count($cycles);
print_r($cycles);

from numpy import *
import sys
import numpy as np
x=1
epsilon=sys.float_info.epsilon/2.0
# machine epsilon
h = np.linspace(1,epsilon,1000)
e1= (((sin(x)- sin(x-h))/h-cos(x))/cos(x))                      #realtive error in first order
e2= (((sin(x-2*h)-4*sin(x-h)+3*sin(x))/(2*h)-cos(x))/cos(x))    #realtive error in second order
e3= (((2*sin(x+h)+3*sin(x)-6*sin(x-h)+sin(x-2*h))/(6*h)-cos(x))/cos(x))#realtive error in third order
e4 = (((-sin(x-3*h)+6*sin(x-2*h)-18*sin(x-h)+10*sin(x)+3*sin(x+h))/(12*h)-cos(x))/cos(x))#realtive error in fourth order
print(h)
import matplotlib.pyplot as plt
plt.plot(h,e1,'-b',lw=3)
plt.plot(h,e2,'-g',lw=3)
plt.plot(h,e3,'-m',lw=3)
plt.plot(h,e4,'-r',lw=3)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('h')
plt.ylabel('Relative error')
plt.legend(['first order','second order','third order','fourth order'])
plt.title('Effect of order of accuracy on the numerical differentiation')
plt.show()

template <typename T>
constexpr MCSRmatrix<T>::MCSRmatrix( std::initializer_list<std::initializer_list<T>> rows)
{
      this->dim  = rows.size();
      auto _rows = *(rows.begin());
            
      aa_.resize(dim+1);
      ja_.resize(dim+1);
      
      if(dim != _rows.size())
      {
          throw InvalidSizeException("Error in costructor! MCSR format require square matrix!");  
      }
      
      itype w = 0 ;
      ja_.at(w) = dim+2 ;
      for(auto ii = rows.begin(), i=1; ii != rows.end() ; ++ii, i++)
      {
          for(auto ij = ii->begin(), j=1, elemCount = 0 ; ij != ii->end() ; ++ij, j++ )   
          {
              if(i==j)
                 aa_[i-1] = *ij ;
              else if( i != j && *ij != 0 )
              {   
                 ja_.push_back(j); 
                 aa_.push_back(*ij); 
                 elemCount++ ;
              }
              ja_[i] = ja_[i-1] + elemCount;           
          }
      }     
}

template <typename T>
std::size_t constexpr MCSRmatrix<T>::findIndex(const itype row ,  const itype col) const noexcept 
{    
     assert( row > 0 && row <= dim && col > 0 && col <= dim ); 
     if(row == col)
     {
       return row-1;
     }
     int i = -1;
     
      for(int i = ja_.at(row-1)-1 ; i < ja_.at(row)-1 ; i++ )
      {
            if( ja_.at(i) == col ) 
            { 
                  
              return i   ;
            }
      }
      return -1;
}

template <typename T>
const T MCSRmatrix<T>::operator()(const itype r , const itype c) const noexcept 
{     
     auto i = findIndex(r,c); 
     if( i != -1 && i < aa_.size() )
     {      
            return aa_.at(i) ;
     }
     else
     {
            return 0.0 ;
     }       
}

template<typename T>
auto constexpr MCSRmatrix<T>::insertAt(itype row ,itype col , T elem) noexcept 
{
      if(elem == 0)
      {
          std::cerr << "zero element to insert ! Exit 1" << std::endl; 
      }
      
      int index = findIndex(row,col);
      if( index >= dim+1  || row == col)
      {     
          aa_.at(index) = elem ;
      }
      else if(index == -1) 
      {

          for(auto i=row; i< dim+1 ; i++) 
          {     
                ja_.at(i) += 1 ;
          }
          
          if(ja_.at(row) >= aa_.size() ) //std::cout << aa_.size() << ' ' <<ja_.at(row) << std::endl;
          {
              aa_.push_back(elem);
              ja_.push_back(col);
          }
          else if ( ja_.at(row) < aa_.size() )
          {
                ja_.insert(ja_.begin() + dim  + row , col );
                aa_.insert(aa_.begin() + dim  + row , elem);
          }
      }
}