Path: unixg.ubc.ca!vanbc.wimsey.com!scipio.cyberstore.ca!math.ohio-state.edu!howland.reston.ans.net!europa.eng.gtefsd.com!paladin.american.edu!auvm!AX.APC.ORG!barbuto
Comments: Gated by NETNEWS@AUVM.AMERICAN.EDU
Newsgroups: bit.listserv.frac-l
X-Envelope-to: FRAC-L@GITVM1.BITNET
Message-ID: <199410301921.LAA21119@ax.ibase.br>
Date: Sun, 30 Oct 1994 11:21:42 -0800
Sender: "\"FRACTAL\" discussion list" <FRAC-L@GITVM1.BITNET>
From: Fausto Arinos de Almeida Barbuto <barbuto@AX.APC.ORG>
Subject: Routine in Assembler
Lines: 142
 
From barbuto  Sun Oct 30 11:12:02 1994
>From barbuto  Sun Oct 30 17:12:02 1994
Received: by ax.ibase.br (8.6.9/Revision: 1.13 )
        id LAA20788; Sun, 30 Oct 1994 11:12:02 -0800
Date: Sun, 30 Oct 1994 11:12:02 -0800
From: Fausto Arinos de Almeida Barbuto <barbuto>
Message-Id: <199410301912.LAA20788@ax.ibase.br>
To: barbuto
Subject: Routine in Assembler
Status: RO
 
 
  October 31, 1994
 
  FRAC-Lers:
 
     Fractal programmers use a wide variety of programming languages,  mainly
  C/C++,  Basic and Pascal.  C has the intrinsic function "pow" to  calculate
  powers of the type x^y (x > 0),  but Pascal users often use Exp(y*Ln(x)) to
  perform these calculations.
 
     One of my workmates, Alberto H. V. Sarmento, has developed an  Assembler
  function to be used with Pascal  which has a  much better  performance than
  the usual Exp(y*Ln(x)). I was able to port it to C; some minor changes were
  enough to make it workable (though I have had a lot of work to do it). Both
  Pascal and C/C++ versions are included at the bottom of this message.
 
     The benchmarks' results are shown below  (Sarmento's function is  called
  POWER, and the machine employed is a DX2-66):
 
  --------------------------------Execution time (seconds)-------------------
  Number of x^y    Pascal & exp(..)   Pascal & POWER    C & pow     C & POWER
  operations
 
   10000           0.16               0.11              0.164835    0.109890
   12000           0.22               0.17              0.219780    0.109890(?)
   20000           0.39               0.28              0.384615    0.219780
   50000           1.11               0.66              1.043956    0.604396
  100000           1.98               1.33              2.417582    1.428571
  ---------------------------------------------------------------------------
 
     This function can be used not only as a fractal programming aid, but al-
  so for any other purpose, PROVIDED x > 0 in x^y. And here are the sources:
 
---PASCAL---Begin---CUT HERE---
Function POWER (b,e:double):double;
var
  status : word;
  resp   : extended;
 
begin
    asm
                    { st          st(1)       st(2)      st(3)    }
    fld e           { e                                           }
    fld b           { b           e                               }
    fyl2x           { e.log(b)=a                                  }
    fld st          { a           a                               }
    frndint         { round(a)=r  a                               }
    fxch            { a           r                               }
    fld st(1)       { r           a           r                   }
    fsubp st(1),st  { a-r         r                               }
    ftst
    fstsw status
    mov   ax, status
    and   ax, $4500
    cmp   ax, $0100                                  { a-r<0.0 ? }
    je    @@1                                        { yes, jump }
    f2xm1           { 2^(a-r)-1   r                              }
    fld1            { 1           2^(a-r)-1   r                  }
    faddp st(1),st  { 2^(a-r)     r                              }
    fscale          { 2^(a-r)*2^r r                              }
    jmp @@2
    @@1:
    fabs
    f2xm1           { 2^(a-r)-1   r                              }
    fld1            { 1           2^(a-r)-1   r                  }
    faddp st(1),st  { 2^(a-r)     r                              }
    fld1
    fdivrp st(1),st
    fscale          { 2^(a-r)*2^r r                              }
    @@2:
    fxch  st(1)     { r           2^(a-r)*2^r                    }
    fstp  resp
    fstp  resp
  end;
  power := resp;
end;
---PASCAL---End---CUT HERE---
 
---C/C++---Begin---CUT HERE---
double POWER(double b, double e)
{
  unsigned int status;
  double resp, mypower;
 
#pragma inline +eff
    asm extrn _fabs:near
                        //* st          st(1)       st(2)      st(3)    *//
    asm fld e           //* e                                           *//
    asm fld b           //* b           e                               *//
    asm fyl2x           //* e.log(b)=a                                  *//
    asm fld st          //* a           a                               *//
    asm frndint         //* round(a)=r  a                               *//
    asm fxch            //* a           r                               *//
    asm fld st(1)       //* r           a           r                   *//
    asm fsubp st(1),st  //* a-r         r                               *//
    asm ftst
    asm fstsw status
    asm mov   ax, status
    asm and   ax, 4500h
    asm cmp   ax, 0100h                                  //* a-r<0.0 ? *//
    asm je    Label1                                     //* yes, jump *//
    asm f2xm1           //* 2^(a-r)-1   r                              *//
    asm fld1            //* 1           2^(a-r)-1   r                  *//
    asm faddp st(1),st  //* 2^(a-r)     r                              *//
    asm fscale          //* 2^(a-r)*2^r r                              *//
    asm jmp Label2
    Label1:
    asm FABS
    asm f2xm1           //* 2^(a-r)-1   r                              *//
    asm fld1            //* 1           2^(a-r)-1   r                  *//
    asm faddp st(1),st  //* 2^(a-r)     r                              *//
    asm fld1
    asm fdivrp st(1),st
    asm fscale          //* 2^(a-r)*2^r r                              *//
    Label2:
    asm fxch  st(1)     //* r           2^(a-r)*2^r                    *//
    asm fstp  resp
    asm fstp  resp
#pragma inline -eff
  mypower = resp;
  return (mypower);
}
---C/C++---End---CUT HERE---
 
  Hope you may enjoy (and use) this...
 
  ------- Fausto Arinos de Almeida Barbuto --------
  BJ06%C53000.PETROBRAS.ANRJ.BR@fpsp.fapesp.ansp.br
  INTERNET:                 barbuto@ax.ibase.org.br
  Rio de  Janeiro, RJ, Federal  Republic  of Brazil
  FAX:(55)(21)5986796/5986441 Phone:(55)(21)5986712
  ----"Friends come and go; enemies accumulate"----