---------------------------------------------------------------------------

                 S O U R C E S    Y   S O U R C E S I T O S
         
  ---------------------------------------------------------------------------


   Nos gustaría en cada número, incluir algunas rutinas, bien en ensamblador o
  en cualquier otro lenguaje, pues creemos que pueden ser interesantes. Así
  que ya sabéis, si tenéis alguna, ¡enviarla ya!.

   En este primer número incluimos una rutina bastante interesante para ex-
  traer raíces cuadradas, tanto en Assembler como en C.


    ; Esta es una rutina super-rápida para extraer la raiz cuadrada de
    ; un número
    ; (En principio de 16 bits)
    ;
    ; © OTTO PETER, Innsbruck


   ROOT:

   ;       SETTINGS IN:      D0.w   -   number
   ;
   ;               OUT:      D0.w   -   result
   ;                                -   sign- and zeroflag are valid

   movem.w      d1-d3,-(sp)
   move.w  #$4000,d2

   root_loop:
   move.w       d1,d3
   add.w        d2,d3

   lsr.w        #1,d1
   cmp.w        d3,d0
   ble.s        boing
   sub.w        d3,d0
   or.w         d2,d1

   boing:
   lsr.w        #2,d2
   bne.s        root_loop

   cmp.w        d1,d0
   blt.s        no_round_up
   addq.w       #1,d1

   no_round_up:
   move.w       d1,d0
   movem.w      (sp)+,d1-d3
   rts


   ;*************************************************
   ;       Rutina para números en 32 bits
   ;*************************************************
   exp_offset:	equ $7f


   domain_error:        ;Error if number is negative
   movem.l      (a7)+,d1-d4
   rts

   sq_0:
   clr.l        d0
   movem.l      (a7)+,d1-d4
   rts

   sqrt:
   movem.l      d1-d4,-(a7)     ;Save d1-d4
   move.l       d0,d4
   bmi.s        domain_error    ;Error if number is negative
   swap	d4                      ;MSW of number
   and.l        #$7f80,d0       ;isolate exponent
   beq.s        sq_0            ;exponent is 0 => root is 0
   and.l        #$007fffff,d0           ;isolate mantisse
   sub.w        #exp_offset*$80,d4      ;exponent in `2^i`-format
   bclr         #7,d4                           ;exponent even?
   beq.s        even_exp
   add.l        d0,d0                   ;mantisse * 2
   add.l        #$01000000-$00800000,d0 ;set hidden-bit

   even_exp:
   asr.w        #1,d4                   ;exponent / 2
   add.w        #exp_offset*$80,d4      ;exponent in offset-format
   swap         d4

   lsl.l        #7,d0

   move.l       #$40000000,d2           ;xroot after 1. interation
   move.l       #$10000000,d3           ;m2 = 2 << (maxbit-1)

   loop1_0:
   move.l       d0,d1                   ;xx2 = x

   loop1_1:
   sub.l        d2,d1                   ;xx2 -= xroot
   lsr.l        #1,d2                   ;xroot >>= 1
   sub.l        d3,d1                   ;x2 -= m2
   bmi.s        dont_set1
   move.l       d1,d0                   ;x = xx2
   or.l         d3,d2                   ;xroot += m2
   lsr.l        #2,d3                   ;m2 >>= 2
   bne.s        loop1_1
   bra.s        d0_d1_same
   
   dont_set1:
   lsr.l        d2,d3                   ;m2 >>= 2
   bne.s        loop1_0                 ;repeat loop 15 times (bit 22..8)

   move.l       d0,d1

   d0_d1_same:
   sub.l        d2,d1                   ;xx2 -= xroot
   ror.l        #1,d2                   ;xroot >>= 1 (with carry)
   swap         d2                      ;turn to new aligment
   subq.l       #1,d1                   ;carry of 0-0x4000: x2 -= m2 (part 1)
   bmi.s        dont_set7

   or.l         #-$40000000,d1                  ;0-0x4000: x2 -= m2 (part 2)
   move.l       d1,d0
   or.w         #$4000,d2

   dont_set7:
   swap         d0                      ;turn x to new aligment

   move.w       #$1000,d3               ;m2 - bit 16..31 = 0

   loop2_0:
   move.l       d0,d1                   ;xx2 = x

   loop2_1:
   sub.l        d2,d1                   ;xx2 = xroot
   lsr.l        #1,d2                   ;xroot >>= 1
   sub.l        d3,d1
   bmi.s        dont_set2

   move.l       d1,d0                   ;x = xx2
   or.l         d3,d2                   ;xroot += m2
   lsr.w        #2,d3                   ;m2 >>= 2
   bne.s        loop2_1

   bra.s        finish

   dont_set2:
   lsr.w        #2,d3                   ;m2 >>= 2
   bne.s        loop2_0                 ;repeat loop 7 times (n=6..0)

   finish:
   sub.l        d2,d0                   ;round root?
   bls.s        no_inc
   addq.l       #1,d2                   ;round up!

   no_inc:
   bclr         #23,d2                  ;clear hidden bit
   or.l         d4,d2                   ;combine exponent and mantisse
   move.l       d2,d0                   ;result
   movem.l      (a7)+,d1-d4
   rts


    ----------------     C     -------------------------------


   #define N_BITS 32
   #define MAX_BIT ((N_BITS + 1) / 2 - 1)

   unsigned long int sqrt_1(x)
     unsigned long int x;
     {
     register int i;
     register unsigned long int m, r, root;

     root = 0;  m = 1 << MAX_BIT
     for (i = MAX_BIT;  i >= 0;  i--)
       {
       r = root + m;
       if (r * r <= X) root = r;
       m >>= 1;                           
       }
     return root;
     }


   unsigned long int sqrt_2(x)
     unsigned long int x;
     {
     register long int i;
     register unsigned long int m, r2, root;
     unsigned long int root2;

     root2 = root = 0;  m = 1 << MAX_BIT;
     for (i = MAX_BIT;  i >= 0;  i--)
       {
       r2 = root2 + (root << i + 1)  +  (1 << i + i);
       if (r2 <= X)
         {
         root2 = r2;  root += m;
         }
       m >>= 1;
       }
     return root;
     }


   unsigned long int sqrt_3(x)
     unsigned long int x;
     {
     register unsigned long int m2, r2, xroot, root2;
     int i;

     root2 = xroot = 0;  m2 = 1 << MAX_BIT * 2;
  for (i = MAX_BIT;  i >= 0;  i--)
       {
       r2 = root2 + xroot + m2;
       xroot >>= 1;
       if (r2 <= X)
         {
      root2 = r2;  xroot += m2;
         }
            m2 >>= 2;
       }
     return xroot;
     }


   unsigned long int sqrt_4(x)
  unsigned long int x;
     {
     register unsigned long int xroot, m2, x2;

  xroot = 0;  m2 = 1 << MAX_BIT * 2;
     do
       {
       x2 = xroot + m2;
       xroot >>= 1;
       if (x2 <= X)
         {
         x -= x2;  xroot += m2;
         }
       }
     while (m2 >>= 2);
    return xroot;
     }


   unsigned long int sqrt_5(x)
     unsigned long int x;
    {
     register unsigned long int xroot, m2, x2;

     xroot = 0;  m2 = 1 << MAX_BIT * 2;
     do
       {
       x2 = xroot + m2;
       xroot >>= 1;
       if (x2 <= X)
         {
         x -= x2;  xroot += m2;
         }
       }
  while (m2 >>= 2);
     if (xroot < x) return xroot + 1;
     return xroot;
  }