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  * $DragonFly: src/sys/i386/gnu/fpemul/Attic/polynomial.s,v 1.3 2003/08/07 21:17:20 dillon Exp $
  69  *
  70  */
  71
  72         .file   "fpolynom.s"
  73
  74 #include "fpu_asm.h"
  75
  76
  77 /*      #define EXTRA_PRECISE*/
  78
  79 #define TERM_SIZE       $8
  80
  81
  82 .text
  83 ENTRY(polynomial)
  84         pushl   %ebp
  85         movl    %esp,%ebp
  86         subl    $32,%esp
  87         pushl   %esi
  88         pushl   %edi
  89         pushl   %ebx
  90
  91         movl    PARAM1,%esi             /* accum */
  92         movl    PARAM2,%edi             /* x */
  93         movl    PARAM3,%ebx             /* terms */
  94         movl    PARAM4,%ecx             /* n */
  95
  96         movl    TERM_SIZE,%eax
  97         mull    %ecx
  98         movl    %eax,%ecx
  99
 100         movl    4(%ebx,%ecx,1),%edx     /* terms[n] */
 101         movl    %edx,-20(%ebp)
 102         movl    (%ebx,%ecx,1),%edx      /* terms[n] */
 103         movl    %edx,-24(%ebp)
 104         xor     %eax,%eax
 105         movl    %eax,-28(%ebp)
 106
 107         subl    TERM_SIZE,%ecx
 108         js      L_accum_done
 109
 110 L_accum_loop:
 111         xor     %eax,%eax
 112         movl    %eax,-4(%ebp)
 113         movl    %eax,-8(%ebp)
 114
 115 #ifdef EXTRA_PRECISE
 116         movl    -28(%ebp),%eax
 117         mull    4(%edi)                 /* x ms long */
 118         movl    %edx,-12(%ebp)
 119 #endif EXTRA_PRECISE
 120
 121         movl    -24(%ebp),%eax
 122         mull    (%edi)                  /* x ls long */
 123 /*      movl    %eax,-16(%ebp)  */      /* Not needed */
 124         addl    %edx,-12(%ebp)
 125         adcl    $0,-8(%ebp)
 126
 127         movl    -24(%ebp),%eax
 128         mull    4(%edi)                 /* x ms long */
 129         addl    %eax,-12(%ebp)
 130         adcl    %edx,-8(%ebp)
 131         adcl    $0,-4(%ebp)
 132
 133         movl    -20(%ebp),%eax
 134         mull    (%edi)
 135         addl    %eax,-12(%ebp)
 136         adcl    %edx,-8(%ebp)
 137         adcl    $0,-4(%ebp)
 138
 139         movl    -20(%ebp),%eax
 140         mull    4(%edi)
 141         addl    %eax,-8(%ebp)
 142         adcl    %edx,-4(%ebp)
 143
 144 /* Now add the next term */
 145         movl    (%ebx,%ecx,1),%eax
 146         addl    %eax,-8(%ebp)
 147         movl    4(%ebx,%ecx,1),%eax
 148         adcl    %eax,-4(%ebp)
 149
 150 /* And put into the second register */
 151         movl    -4(%ebp),%eax
 152         movl    %eax,-20(%ebp)
 153         movl    -8(%ebp),%eax
 154         movl    %eax,-24(%ebp)
 155
 156 #ifdef EXTRA_PRECISE
 157         movl    -12(%ebp),%eax
 158         movl    %eax,-28(%ebp)
 159 #else
 160         testb   $128,-25(%ebp)
 161         je      L_no_poly_round
 162
 163         addl    $1,-24(%ebp)
 164         adcl    $0,-20(%ebp)
 165 L_no_poly_round:
 166 #endif EXTRA_PRECISE
 167
 168         subl    TERM_SIZE,%ecx
 169         jns     L_accum_loop
 170
 171 L_accum_done:
 172 #ifdef EXTRA_PRECISE
 173 /* And round the result */
 174         testb   $128,-25(%ebp)
 175         je      L_poly_done
 176
 177         addl    $1,-24(%ebp)
 178         adcl    $0,-20(%ebp)
 179 #endif EXTRA_PRECISE
 180
 181 L_poly_done:
 182         movl    -24(%ebp),%eax
 183         movl    %eax,(%esi)
 184         movl    -20(%ebp),%eax
 185         movl    %eax,4(%esi)
 186
 187         popl    %ebx
 188         popl    %edi
 189         popl    %esi
 190         leave
 191         ret