sys/i386/gnu/fpemul/polynomial.s

   1 /*
   2  *  polynomial.S
   3  *
   4  * Fixed point arithmetic polynomial evaluation.
   5  *
   6  * Call from C as:
   7  *   void polynomial(unsigned accum[], unsigned x[], unsigned terms[][2],
   8  *                   int n)
   9  *
  10  * Computes:
  11  * terms[0] + (terms[1] + (terms[2] + ... + (terms[n-1]*x)*x)*x)*x) ... )*x
  12  * The result is returned in accum.
  13  *
  14  *
  15  * Copyright (C) 1992,1993,1994
  16  *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
  17  *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
  18  * All rights reserved.
  19  *
  20  * This copyright notice covers the redistribution and use of the
  21  * FPU emulator developed by W. Metzenthen. It covers only its use
  22  * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
  23  * use is not permitted under this copyright.
  24  *
  25  * Redistribution and use in source and binary forms, with or without
  26  * modification, are permitted provided that the following conditions
  27  * are met:
  28  * 1. Redistributions of source code must retain the above copyright
  29  *    notice, this list of conditions and the following disclaimer.
  30  * 2. Redistributions in binary form must include information specifying
  31  *    that source code for the emulator is freely available and include
  32  *    either:
  33  *      a) an offer to provide the source code for a nominal distribution
  34  *         fee, or
  35  *      b) list at least two alternative methods whereby the source
  36  *         can be obtained, e.g. a publically accessible bulletin board
  37  *         and an anonymous ftp site from which the software can be
  38  *         downloaded.
  39  * 3. All advertising materials specifically mentioning features or use of
  40  *    this emulator must acknowledge that it was developed by W. Metzenthen.
  41  * 4. The name of W. Metzenthen may not be used to endorse or promote
  42  *    products derived from this software without specific prior written
  43  *    permission.
  44  *
  45  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  46  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
  47  * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
  48  * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  49  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  50  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  51  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  52  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  53  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  54  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  55  *
  56  *
  57  * The purpose of this copyright, based upon the Berkeley copyright, is to
  58  * ensure that the covered software remains freely available to everyone.
  59  *
  60  * The software (with necessary differences) is also available, but under
  61  * the terms of the GNU copyleft, for the Linux operating system and for
  62  * the djgpp ms-dos extender.
  63  *
  64  * W. Metzenthen   June 1994.
  65  *
  66  *
  67  * $FreeBSD: src/sys/gnu/i386/fpemul/polynomial.s,v 1.8 1999/08/28 00:42:55 peter Exp $
  68  *
  69  */
  70
  71         .file   "fpolynom.s"
  72
  73 #include <gnu/i386/fpemul/fpu_asm.h>
  74
  75
  76 /*      #define EXTRA_PRECISE*/
  77
  78 #define TERM_SIZE       $8
  79
  80
  81 .text
  82 ENTRY(polynomial)
  83         pushl   %ebp
  84         movl    %esp,%ebp
  85         subl    $32,%esp
  86         pushl   %esi
  87         pushl   %edi
  88         pushl   %ebx
  89
  90         movl    PARAM1,%esi             /* accum */
  91         movl    PARAM2,%edi             /* x */
  92         movl    PARAM3,%ebx             /* terms */
  93         movl    PARAM4,%ecx             /* n */
  94
  95         movl    TERM_SIZE,%eax
  96         mull    %ecx
  97         movl    %eax,%ecx
  98
  99         movl    4(%ebx,%ecx,1),%edx     /* terms[n] */
 100         movl    %edx,-20(%ebp)
 101         movl    (%ebx,%ecx,1),%edx      /* terms[n] */
 102         movl    %edx,-24(%ebp)
 103         xor     %eax,%eax
 104         movl    %eax,-28(%ebp)
 105
 106         subl    TERM_SIZE,%ecx
 107         js      L_accum_done
 108
 109 L_accum_loop:
 110         xor     %eax,%eax
 111         movl    %eax,-4(%ebp)
 112         movl    %eax,-8(%ebp)
 113
 114 #ifdef EXTRA_PRECISE
 115         movl    -28(%ebp),%eax
 116         mull    4(%edi)                 /* x ms long */
 117         movl    %edx,-12(%ebp)
 118 #endif EXTRA_PRECISE
 119
 120         movl    -24(%ebp),%eax
 121         mull    (%edi)                  /* x ls long */
 122 /*      movl    %eax,-16(%ebp)  */      /* Not needed */
 123         addl    %edx,-12(%ebp)
 124         adcl    $0,-8(%ebp)
 125
 126         movl    -24(%ebp),%eax
 127         mull    4(%edi)                 /* x ms long */
 128         addl    %eax,-12(%ebp)
 129         adcl    %edx,-8(%ebp)
 130         adcl    $0,-4(%ebp)
 131
 132         movl    -20(%ebp),%eax
 133         mull    (%edi)
 134         addl    %eax,-12(%ebp)
 135         adcl    %edx,-8(%ebp)
 136         adcl    $0,-4(%ebp)
 137
 138         movl    -20(%ebp),%eax
 139         mull    4(%edi)
 140         addl    %eax,-8(%ebp)
 141         adcl    %edx,-4(%ebp)
 142
 143 /* Now add the next term */
 144         movl    (%ebx,%ecx,1),%eax
 145         addl    %eax,-8(%ebp)
 146         movl    4(%ebx,%ecx,1),%eax
 147         adcl    %eax,-4(%ebp)
 148
 149 /* And put into the second register */
 150         movl    -4(%ebp),%eax
 151         movl    %eax,-20(%ebp)
 152         movl    -8(%ebp),%eax
 153         movl    %eax,-24(%ebp)
 154
 155 #ifdef EXTRA_PRECISE
 156         movl    -12(%ebp),%eax
 157         movl    %eax,-28(%ebp)
 158 #else
 159         testb   $128,-25(%ebp)
 160         je      L_no_poly_round
 161
 162         addl    $1,-24(%ebp)
 163         adcl    $0,-20(%ebp)
 164 L_no_poly_round:
 165 #endif EXTRA_PRECISE
 166
 167         subl    TERM_SIZE,%ecx
 168         jns     L_accum_loop
 169
 170 L_accum_done:
 171 #ifdef EXTRA_PRECISE
 172 /* And round the result */
 173         testb   $128,-25(%ebp)
 174         je      L_poly_done
 175
 176         addl    $1,-24(%ebp)
 177         adcl    $0,-20(%ebp)
 178 #endif EXTRA_PRECISE
 179
 180 L_poly_done:
 181         movl    -24(%ebp),%eax
 182         movl    %eax,(%esi)
 183         movl    -20(%ebp),%eax
 184         movl    %eax,4(%esi)
 185
 186         popl    %ebx
 187         popl    %edi
 188         popl    %esi
 189         leave
 190         ret