development/compilers/freepascal/rtl/i386/math.inc

   1 {
   2     $Id$
   3     This file is part of the Free Pascal run time library.
   4     Copyright (c) 1999-2000 by the Free Pascal development team
   5
   6     Implementation of mathamatical Routines (only for real)
   7
   8     See the file COPYING.FPC, included in this distribution,
   9     for details about the copyright.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  14
  15  **********************************************************************}
  16
  17
  18 {****************************************************************************
  19                        EXTENDED data type routines
  20  ****************************************************************************}
  21
  22 {$ifdef hasinternmath}
  23     function pi : extended;[internproc:in_pi];
  24     function abs(d : extended) : extended;[internproc:in_abs_extended];
  25     function sqr(d : extended) : extended;[internproc:in_sqr_extended];
  26     function sqrt(d : extended) : extended;[internproc:in_sqrt_extended];
  27     function arctan(d : extended) : extended;[internproc:in_arctan_extended];
  28     function ln(d : extended) : extended;[internproc:in_ln_extended];
  29     function sin(d : extended) : extended;[internproc:in_sin_extended];
  30     function cos(d : extended) : extended;[internproc:in_cos_extended];
  31 {$else hasinternmath}
  32     function pi : extended;assembler;[internconst:in_const_pi];
  33       asm
  34             fldpi
  35       end [];
  36
  37
  38     function abs(d : extended) : extended;assembler;[internconst:in_const_abs];
  39       asm
  40             fldt d
  41             fabs
  42       end [];
  43
  44
  45     function sqr(d : extended) : extended;assembler;[internconst:in_const_sqr];
  46       asm
  47             fldt d
  48             fldt d
  49             fmulp %st(1)
  50       end [];
  51
  52
  53     function sqrt(d : extended) : extended;assembler;[internconst:in_const_sqrt];
  54       asm
  55             fldt d
  56             fsqrt
  57       end [];
  58
  59
  60     function arctan(d : extended) : extended;assembler;[internconst:in_const_arctan];
  61       asm
  62             fldt d
  63             fld1
  64             fpatan
  65       end [];
  66
  67     function cos(d : extended) : extended;assembler;[internconst:in_const_cos];
  68       asm
  69             fldt d
  70             fcos
  71             fstsw
  72             sahf
  73             jnp .LCOS1
  74             fstp %st(0)
  75             fldt .LCOS0
  76             jmp .LCOS1
  77          .data
  78          .LCOS0:
  79             .long       0xffffffff
  80             .long       0xffffffff
  81             .long       0xffffffff
  82          .text
  83          .LCOS1:
  84       end ['EAX'];
  85
  86     function ln(d : extended) : extended;assembler;[internconst:in_const_ln];
  87       asm
  88             fldln2
  89             fldt d
  90             fyl2x
  91       end [];
  92
  93
  94     function sin(d : extended) : extended;assembler;[internconst:in_const_sin];
  95       asm
  96             fldt d
  97             fsin
  98             fstsw
  99             sahf
 100             jnp .LSIN1
 101             fstp %st(0)
 102             fldt .LSIN0
 103             jmp .LSIN1
 104          .data
 105          .LSIN0:
 106             .long       0xffffffff
 107             .long       0xffffffff
 108             .long       0xffffffff
 109          .text
 110          .LSIN1:
 111       end ['EAX'];
 112
 113
 114 {$endif hasinternmath}
 115
 116     function exp(d : extended) : extended;assembler;[internconst:in_const_exp];
 117        asm
 118             // comes from DJ GPP
 119             fldt        d
 120             fldl2e
 121             fmulp       %st(1)
 122             fstcw      .LCW1
 123             fstcw      .LCW2
 124             andw        $0xf3ff,.LCW2
 125             orw         $0x0400,.LCW2
 126             fldcw      .LCW2
 127             fld         %st(0)
 128             frndint
 129             fldcw      .LCW1
 130             fxch        %st(1)
 131             fsub        %st(1),%st
 132             f2xm1
 133             fld1
 134             faddp       %st(1)
 135             fscale
 136             fstp        %st(1)
 137             jmp         .LCW3
 138             // store some help data in the data segment
 139         .data
 140         .LCW1:
 141             .word       0
 142         .LCW2:
 143             .word       0
 144         .text
 145         .LCW3:
 146       end;
 147
 148
 149     function frac(d : extended) : extended;assembler;[internconst:in_const_frac];
 150       asm
 151             subl $16,%esp
 152             fnstcw -4(%ebp)
 153             fwait
 154             movw -4(%ebp),%cx
 155             orw $0x0c3f,%cx
 156             movw %cx,-8(%ebp)
 157             fldcw -8(%ebp)
 158             fwait
 159             fldt d
 160             frndint
 161             fldt d
 162             fsub %st(1)
 163             fstp %st(1)
 164             fclex
 165             fldcw -4(%ebp)
 166       end ['ECX'];
 167
 168
 169     function int(d : extended) : extended;assembler;[internconst:in_const_int];
 170       asm
 171             subl $16,%esp
 172             fnstcw -4(%ebp)
 173             fwait
 174             movw -4(%ebp),%cx
 175             orw $0x0c3f,%cx
 176             movw %cx,-8(%ebp)
 177             fldcw -8(%ebp)
 178             fwait
 179             fldt d
 180             frndint
 181             fclex
 182             fldcw -4(%ebp)
 183       end ['ECX'];
 184
 185
 186     function trunc(d : extended) : longint;assembler;[internconst:in_const_trunc];
 187       asm
 188             subl $16,%esp
 189             fnstcw -4(%ebp)
 190             fwait
 191             movw -4(%ebp),%cx
 192             orw $0x0c3f,%cx
 193             movw %cx,-8(%ebp)
 194             fldcw -8(%ebp)
 195             fwait
 196             fldt d
 197             fistpl -8(%ebp)
 198             movl -8(%ebp),%eax
 199             fldcw -4(%ebp)
 200       end ['EAX','ECX'];
 201
 202
 203     function round(d : extended) : longint;assembler;[internconst:in_const_round];
 204       asm
 205             subl $8,%esp
 206             fnstcw -4(%ebp)
 207             fwait
 208             movw $0x1372,-8(%ebp)
 209             fldcw -8(%ebp)
 210             fwait
 211             fldt d
 212             fistpl -8(%ebp)
 213             movl -8(%ebp),%eax
 214             fldcw -4(%ebp)
 215       end ['EAX','ECX'];
 216
 217
 218    function power(bas,expo : extended) : extended;
 219      begin
 220         if bas=0 then
 221           begin
 222             if expo<>0 then
 223               power:=0.0
 224             else
 225               HandleError(207);
 226           end
 227         else if expo=0 then
 228          power:=1
 229         else
 230         { bas < 0 is not allowed }
 231          if bas<0 then
 232           handleerror(207)
 233          else
 234           power:=exp(ln(bas)*expo);
 235      end;
 236
 237
 238 {****************************************************************************
 239                        Longint data type routines
 240  ****************************************************************************}
 241
 242    function power(bas,expo : longint) : longint;
 243      begin
 244         if bas=0 then
 245           begin
 246             if expo<>0 then
 247               power:=0
 248             else
 249               HandleError(207);
 250           end
 251         else if expo=0 then
 252          power:=1
 253         else
 254          begin
 255            if bas<0 then
 256             begin
 257               if odd(expo) then
 258                power:=-round(exp(ln(-bas)*expo))
 259               else
 260                power:=round(exp(ln(-bas)*expo));
 261             end
 262            else
 263             power:=round(exp(ln(bas)*expo));
 264          end;
 265      end;
 266
 267
 268 {****************************************************************************
 269                          Fixed data type routines
 270  ****************************************************************************}
 271
 272 {$ifdef HASFIXED} { Not yet allowed }
 273
 274     function sqrt(d : fixed) : fixed;
 275
 276       begin
 277          asm
 278             movl d,%eax
 279             movl %eax,%ebx
 280             movl %eax,%ecx
 281             jecxz .L_kl
 282             xorl %esi,%esi
 283          .L_it:
 284             xorl %edx,%edx
 285             idivl %ebx
 286             addl %ebx,%eax
 287             shrl $1,%eax
 288             subl %eax,%esi
 289             cmpl $1,%esi
 290             jbe .L_kl
 291             movl %eax,%esi
 292             movl %eax,%ebx
 293             movl %ecx,%eax
 294             jmp .L_it
 295          .L_kl:
 296             shl $8,%eax
 297             leave
 298             ret $4
 299          end;
 300       end;
 301
 302
 303     function int(d : fixed) : fixed;
 304     {*****************************************************************}
 305     { Returns the integral part of d                                  }
 306     {*****************************************************************}
 307     begin
 308       int:=d and $ffff0000;       { keep only upper bits      }
 309     end;
 310
 311
 312     function trunc(d : fixed) : longint;
 313     {*****************************************************************}
 314     { Returns the Truncated integral part of d                        }
 315     {*****************************************************************}
 316     begin
 317       trunc:=longint(integer(d shr 16));   { keep only upper 16 bits  }
 318     end;
 319
 320     function frac(d : fixed) : fixed;
 321     {*****************************************************************}
 322     { Returns the Fractional part of d                                }
 323     {*****************************************************************}
 324     begin
 325       frac:=d AND $ffff;         { keep only decimal parts - lower 16 bits }
 326     end;
 327
 328     function abs(d : fixed) : fixed;
 329     {*****************************************************************}
 330     { Returns the Absolute value of d                                 }
 331     {*****************************************************************}
 332     begin
 333        asm
 334            movl d,%eax
 335            rol $16,%eax             { Swap high & low word.}
 336            {Absolute value: Invert all bits and increment when <0 .}
 337            cwd                      { When ax<0, dx contains $ffff}
 338            xorw %dx,%ax             { Inverts all bits when dx=$ffff.}
 339            subw %dx,%ax             { Increments when dx=$ffff.}
 340            rol $16,%eax             { Swap high & low word.}
 341            leave
 342            ret $4
 343        end;
 344     end;
 345
 346
 347     function sqr(d : fixed) : fixed;
 348     {*****************************************************************}
 349     { Returns the Absolute squared value of d                         }
 350     {*****************************************************************}
 351     begin
 352             {16-bit precision needed, not 32 =)}
 353        sqr := d*d;
 354 {       sqr := (d SHR 8 * d) SHR 8; }
 355     end;
 356
 357
 358     function Round(x: fixed): longint;
 359     {*****************************************************************}
 360     { Returns the Rounded value of d as a longint                     }
 361     {*****************************************************************}
 362     var
 363      lowf:integer;
 364      highf:integer;
 365     begin
 366       lowf:=x and $ffff;       { keep decimal part ... }
 367       highf :=integer(x shr 16);
 368       if lowf > 5 then
 369         highf:=highf+1
 370       else
 371       if lowf = 5 then
 372       begin
 373         { here we must check the sign ...       }
 374         { if greater or equal to zero, then     }
 375         { greater value will be found by adding }
 376         { one...                                }
 377          if highf >= 0 then
 378            Highf:=Highf+1;
 379       end;
 380       Round:= longint(highf);
 381     end;
 382
 383 {$endif HASFIXED}
 384
 385 {
 386   $Log$
 387   Revision 1.1  2002/02/19 08:25:18  sasu
 388   Initial revision
 389
 390   Revision 1.1  2000/07/13 06:30:42  michael
 391   + Initial import
 392
 393   Revision 1.23  2000/05/02 10:37:50  pierre
 394    * 0**n where n<>0 is 0; 0**0 generates RTE 207
 395
 396   Revision 1.22  2000/04/07 21:29:00  pierre
 397    changed to get nasm to compile system
 398
 399   Revision 1.21  2000/02/15 14:37:36  florian
 400     * disabled FIXED data type per default
 401
 402   Revision 1.20  2000/02/09 16:59:29  peter
 403     * truncated log
 404
 405   Revision 1.19  2000/01/07 16:41:33  daniel
 406     * copyright 2000
 407
 408   Revision 1.18  2000/01/07 16:32:24  daniel
 409     * copyright 2000 added
 410
 411   Revision 1.17  1999/10/06 17:44:43  peter
 412     * fixed power(int,int) with negative base
 413     * power(ext,ext) with negative base gives rte 207
 414
 415   Revision 1.16  1999/09/15 20:24:11  florian
 416     * some math functions are now coded inline by the compiler
 417 }