Egyiptomi szorzás

Ismert történelem: Nílus áradása után újra ki kell jelölni a határokat.
Papok egy kiválasztott csoportja, a "kötélfeszítők" foglalkoztak evvel.
Kevés írásos emlék maradt fent, az egyik a Rhind matematikai papirusz
(Rhind skót gyüjtő volt), i.e.1650-ből, amit Ahmes írnok jegyzett le.
A tekercs tartalmaz szorzási és osztási algoritmust is.

Az egyiptomiak nem használtak helyiértékeket és nem imerték a nullát.
Ennek eredményeképpen a műveletek, pl. a szorzás nehéz volt és csak
kevesen tudták végrehajtani.

A szorzás: ismételt összeadás.


            1a = a

        (n+1)a = na + a


In C++:

#include <iostream>
#include <cstdlib>

int mul0( int n, int a)
{
  if ( 1 == n ) return a;
  return mul0(n-1,a) + a;
}

int main(int argc, char *argv[])
{
  if ( argc == 3)
  {
    int n = atoi(argv[1]);
    int a = atoi(argv[2]);

    std::cout << mu10(n,a) << std::endl;
  }
  return 0;
}

Természetesen a>0 és n>0, az egyiptomiak tudomásunk szerint nem használtak
negatív számokat.

Az egyiptomi szorzás (Knuth: orosz paraszt algoritmus) Ahmes példája nyomán:

        n = 41  a = 59

         1 #      59  a
         2       118 (a+a)
         4       236 (a+a) + (a+a)
         8 #     472
        16       944
        32 #    1888

        41 x 59 = (1x59)+(8x59)+(32x59)

Figyeljük meg a #-jellel jelölt számokat: ez a 41 bináris reprezentációja:

        00000000 000101001 (16 biten)

Az algoritmus ellenőrzi, hogy "n" páros vagy páratlan, ezért feltételezhetjük,
bár biztosan nem tudjuk, hogy az egyiptomiak ismerték a páros és páratlan szám
fogalmát. A görögök ismerték.

    XXXXXXXXXX  ->   XXXXX XXXXX  páros
    XXXXXXXXXXX ->  XXXXX X XXXXX páratlan




#include <iostream>
#include <cstdlib>

bool odd(int n)  { return  n & 0x1; } // utolso bit
int  half(int n) { return n >> 1; }   // maradek nelkul

int mul1( int n, int a)
{
  if ( 1 == n ) return a;

  int result = mul1(half(n),a+a);
  if ( odd(n) ) result = result + a;

  return result;
}

int main(int argc, char *argv[])
{
  if ( argc == 3)
  {
    int n = atoi(argv[1]);
    int a = atoi(argv[2]);

    std::cout << mul1(n,a) << std::endl;
  }
  return 0;
}


Ennek az algoritmusnak a költsége = log(n) + (v(n)-1)
ahol v(n) az 1-esek száma (population count) n bináris ábrázolásában.


- Megj1:  3 x 15 könnyebben számítható, mint 15 x 3, ez könnyű optimizálás

- Megj2:  A függvényhívások költségesek



Dolgozzunk kicsit többet, számítsuk ki "na" helyett "r + na" -t.
Gyakran már hardver szinten erre a "multiply-accumulate"-re optimalizálnak.



#include <iostream>
#include <cstdlib>

bool odd(int n)  { return  n & 0x1; } // utolso bit
int  half(int n) { return n >> 1; }   // maradek nelkul

int mulacc0( int r, int n, int a)
{
  if ( 1 == n ) return r + a;

  if ( odd(n) )
    return mulacc0(r+a, half(n), a+a);
  else
    return mulacc0(r, half(n), a+a);
}

int main(int argc, char *argv[])
{
  if ( argc == 3)
  {
    int n = atoi(argv[1]);
    int a = atoi(argv[2]);

    std::cout << mulacc0(0,n,a) << std::endl;
  }
  return 0;
}


Kis átalakítással:


int mulacc1( int r, int n, int a)
{
  if ( 1 == n ) return r + a;

  if ( odd(n) )
    r = r + a;
  return mulacc1(r, half(n), a+a);    // tail recursive
}


Továbbá,
  - n ritkán 1
  - ha n páros, akkor nem kell ellenőrizni, hogy 1-e



int mulacc2( int r, int n, int a)
{
  if ( odd(n) )
  {
    r = r + a;
    if ( 1 == n ) return r;
  }
  return mulacc2(r, half(n), a+a);    // vég-rekurzív (tail-recursive)
}



Def: egy algoritmus szigorúan vég-rekurzív, ha a vég-rekurzív hívása
pont ugyanazokkal a formális paraméterekkel történik.


int mulacc2( int r, int n, int a)
{
  if ( odd(n) )
  {
    r = r + a;
    if ( 1 == n ) return r;
  }
  n = half(n);
  a = a + a;
  return mulacc2(r, n, a);    // szigorúan vég-rekurzív
}


Ez már könnyen alakítható ciklussá:


int mulacc3( int r, int n, int a)
{
  while (true )
  {
    if ( odd(n) )
    {
      r = r + a;
      if ( 1 == n ) return r;
    }
    n = half(n);
    a = a + a;
  }
}







Mérjünk!




#include <iostream>
#include <cstdlib>

int eqcounter = 0;
int pluscounter = 0;
int minuscounter = 0;

struct MyInt
{
  MyInt( int i ) : value(i) {}
  int value;

};

MyInt operator+(MyInt a, MyInt b)
{
  ++pluscounter;
  return MyInt(a.value + b.value);
}

MyInt operator-(MyInt a, int b)
{
  ++minuscounter;
  return MyInt(a.value - b);
}

bool operator==(int a, MyInt b)
{
  ++eqcounter;
  return a == b.value;
}

std::ostream& operator<<( std::ostream& os, MyInt n)
{
  os << n.value;
  return os;
}

bool  odd(MyInt n)  { return  n.value & 0x1; }         // utolso bit
MyInt half(MyInt n) { return  MyInt(n.value >> 1); }   // maradek nelkul

MyInt mulacc3( MyInt r, MyInt n, MyInt a)
{
  while (true )
  {
    if ( odd(n) )
    {
      r = r + a;
      if ( 1 == n ) return r;
    }
    n = half(n);
    a = a + a;
  }
}

MyInt multiply( MyInt n, MyInt a)
{
  while ( !odd(n) )
  {
    a = a + a;
    n = half(n);
  }
  if ( 1 == n ) return a;
  return mulacc3( a, n-1, a);
}

int main(int argc, char *argv[])
{
  if ( argc == 3)
  {
    MyInt n = atoi(argv[1]);
    MyInt a = atoi(argv[2]);

    std::cout << multiply(n,a) << std::endl;
    std::cout << "#plus = " << pluscounter << std::endl;
    std::cout << "#minus = " << minuscounter << std::endl;
    std::cout << "#eq = " << eqcounter << std::endl;
  }
  return 0;
}



$ ./a.out 7 8
56
#plus = 4
#minus = 1
#eq = 3




# elso verzio
$ ./a.out 76543 21098
1614904214
#plus = 76542
#minus = 76542
#eq = 76543



# kozbulso verzio
$ ./a.out 76543 21098
1614904214
#plus = 27
#minus = 0
#eq = 17



#utolso verzio
$ ./a.out 76543 21098
1614904214
#plus = 27
#minus = 1
#eq = 12
$