1
// Copyright (c) 1999-2013 Regents of the University of California
3
// FFTW: Copyright (c) 2003,2006 Matteo Frigo
4
// Copyright (c) 2003,2006 Massachusets Institute of Technology
6
// fft8g.[cpp,h]: Copyright (c) 1995-2001 Takya Ooura
8
// ASMLIB: Copyright (c) 2004 Agner Fog
10
// This program is free software; you can redistribute it and/or modify it
11
// under the terms of the GNU General Public License as published by the
12
// Free Software Foundation; either version 2, or (at your option) any later
15
// This program is distributed in the hope that it will be useful, but WITHOUT
16
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
20
// You should have received a copy of the GNU General Public License along
21
// with this program; see the file COPYING. If not, write to the Free Software
22
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24
// In addition, as a special exception, the Regents of the University of
25
// California give permission to link the code of this program with libraries
26
// that provide specific optimized fast Fourier transform (FFT) functions
27
// as an alternative to FFTW and distribute a linked executable and
28
// source code. You must obey the GNU General Public License in all
29
// respects for all of the code used other than the FFT library itself.
30
// Any modification required to support these libraries must be distributed
31
// under the terms of this license. If you modify this program, you may extend
32
// this exception to your version of the program, but you are not obligated to
33
// do so. If you do not wish to do so, delete this exception statement from
34
// your version. Please be aware that FFTW and ASMLIB are not covered by
35
// this exception, therefore you may not use FFTW and ASMLIB in any derivative
36
// work so modified without permission of the authors of those packages.
39
#if defined(__arm__) && !defined(_RC_CHOP)
42
#include "sighandler.h"
45
#define _EM_INVALID 0x00000100 // Invalid Operation
46
#define _EM_ZERODIVIDE 0x00000200 // Divide by zero
47
#define _EM_OVERFLOW 0x00000400 // Overflow
48
#define _EM_UNDERFLOW 0x00000800 // Underflow
49
#define _EM_INEXACT 0x00001000 // Inexact result
50
#define _EM_DENORMAL 0x00008000 // Denormal result
51
#define _MCW_EM (_EM_INVALID|_EM_ZERODIVIDE|_EM_OVERFLOW|_EM_UNDERFLOW|_EM_INEXACT|_EM_DENORMAL)
54
#define _RC_CHOP 0x000c0000 // Truncate
55
#define _RC_DOWN 0x00080000 // Round down
56
#define _RC_UP 0x00040000 // Round up
57
#define _RC_NEAREST 0x00000000 // Round to nearest
58
#define _MCW_RC (_RC_CHOP|_RC_UP|_RC_DOWN|_RC_NEAREST)
60
#define _NAN_DEFAULT 0x00200000 // Use default NaN rather than propogating NaN
61
#define _MCW_NAN _NAN_DEFAULT
63
#define _FLUSH_TO_ZERO 0x00100000 // Replaces denormalized numbers with zero
64
#define _MCW_FLUSH_TO_ZERO _FLUSH_TO_ZERO
67
These VFP vector mode has been deprecated
68
#define _VECTOR_LEN(x) ((((x)-1)&7)<<16) // Set length of vectors to x
69
#define _MCW_VECTOR_LEN_MASK _VECTOR_LEN(8)
71
#define _VECTOR_STRIDE(x) (((x)==2)?(3<<20):0)
72
#define _MCW_VECTOR_STRIDE_MASK (3<<20)
75
static unsigned int save_cw;
76
static unsigned int cw;
78
inline static unsigned int controlfp(unsigned int flags, unsigned int mask) {
81
if (sigsetjmp(jb,1)) {
82
uninstall_sighandler();
85
fprintf(stderr,"User mode control of fp control register not allowed\n");
87
#if defined(__VFP_FP__) && !defined(__SOFTFP__)
88
__asm__ __volatile__ (
95
cw=(cw & ~mask) | (flags & mask);
97
#if defined(__VFP_FP__) && !defined(__SOFTFP__)
98
__asm__ __volatile__ (
103
uninstall_sighandler();
108
inline static unsigned int restorefp() {
110
#if defined(__VFP_FP__) && !defined(__SOFTFP__)
111
__asm__ __volatile__ (
122
static const uint64_t arm_TWO_TO_52(0x4330000000000000);
123
static const uint64_t arm_SIGN_BIT(0x8000000000000000);
124
static const uint32_t arm_TWO_TO_23(0x4b000000);
125
static const uint32_t arm_FSIGN_BIT(0x80000000);
127
inline static double arm_round(double x) {
128
register uint64_t s=*reinterpret_cast<uint64_t *>(&x) & arm_SIGN_BIT;
129
uint64_t a=s | arm_TWO_TO_52;
130
register double d=*reinterpret_cast<double *>(&a);
134
inline static float arm_round(float x) {
135
register uint32_t s=*reinterpret_cast<uint32_t *>(&x) & arm_FSIGN_BIT;
136
uint32_t a=s | arm_TWO_TO_23;
137
register float d=*reinterpret_cast<float *>(&a);
141
#define round(x) arm_round(x)
143
inline static double arm_floor(double x) {
147
inline static float arm_floor(float x) {
148
return round(x-0.5f);
151
#define floor(x) arm_floor(x)