// ************************************************************************** //
//                                                                            //
//    eses                   eses                                             //
//   eses                     eses                                            //
//  eses    eseses  esesese    eses   Embedded Systems Group                  //
//  ese    ese  ese ese         ese                                           //
//  ese    eseseses eseseses    ese   Department of Computer Science          //
//  eses   eses          ese   eses                                           //
//   eses   eseses  eseseses  eses    University of Kaiserslautern            //
//    eses                   eses                                             //
//                                                                            //
// ************************************************************************** //
// Compute the integer square root according to Heron's iteration. This is    //
// also the Newton iteration for computing the root of of f(x) = x^2-a, which //
// is in general x_{i+1} := x_i - (f(x_i) / f'(x_i)), and in our              //
// case x_{i+1} := (x_i + (a/x_i))/2. The iteration converges quadratically,  //
// so that the number of correct digits doubles in each step.                 //
// It can be proved that the sequence is monotonically decreasing for real    //
// numbers. For integer numbers, it is also decreasing until the root is      //
// found. Then, it may either be stable or oscillate between r and r+1.       //
// Since we wish to find the largest number root where root*root<=a holds, we //
// then select the smaller one of the oscillating values.                     //
// ************************************************************************** //

macro N = 2000000000000000000000000000001;

module Heron(nat{N} ?a, root,event !rdy) {
    nat{N} ra,new_root;
    ra = a;
    new_root = a;
    if(a>0)
        do {
            next(new_root) = (new_root + (ra / new_root)) / 2;
            next(root) = new_root;
            pause;
        } while(new_root < root);
    emit(rdy);
    assert(root*root <= ra & ra <= (root+1)*(root+1));
}
drivenby {
    a = 2000000000000000000000000000000;
    await(rdy);
}