src/libm/k_rem_pio2.c
author Ryan C. Gordon <icculus@icculus.org>
Fri, 12 Aug 2016 19:59:00 -0400
changeset 10266 c09f06c4e8c8
parent 8670 0c15c8a2f8c3
permissions -rw-r--r--
emscripten: send fake mouse events for touches, like other targets do. (This really should be handled at the higher level and not in the individual targets, but this fixes the immediate bug.)
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     1
/* @(#)k_rem_pio2.c 5.1 93/09/24 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     2
/*
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     3
 * ====================================================
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     4
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     5
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     6
 * Developed at SunPro, a Sun Microsystems, Inc. business.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     7
 * Permission to use, copy, modify, and distribute this
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     8
 * software is freely granted, provided that this notice
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
     9
 * is preserved.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    10
 * ====================================================
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    11
 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    12
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    13
#if defined(LIBM_SCCS) && !defined(lint)
3162
dc1eb82ffdaa Von: Thomas Zimmermann
Sam Lantinga <slouken@libsdl.org>
parents: 2757
diff changeset
    14
static const char rcsid[] =
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    15
    "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    16
#endif
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    17
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    18
/*
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    19
 * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    20
 * double x[],y[]; int e0,nx,prec; int ipio2[];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    21
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    22
 * __kernel_rem_pio2 return the last three digits of N with
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    23
 *		y = x - N*pi/2
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    24
 * so that |y| < pi/2.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    25
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    26
 * The method is to compute the integer (mod 8) and fraction parts of
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    27
 * (2/pi)*x without doing the full multiplication. In general we
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    28
 * skip the part of the product that are known to be a huge integer (
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    29
 * more accurately, = 0 mod 8 ). Thus the number of operations are
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    30
 * independent of the exponent of the input.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    31
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    32
 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    33
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    34
 * Input parameters:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    35
 * 	x[]	The input value (must be positive) is broken into nx
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    36
 *		pieces of 24-bit integers in double precision format.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    37
 *		x[i] will be the i-th 24 bit of x. The scaled exponent
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    38
 *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    39
 *		match x's up to 24 bits.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    40
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    41
 *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    42
 *			e0 = ilogb(z)-23
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    43
 *			z  = scalbn(z,-e0)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    44
 *		for i = 0,1,2
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    45
 *			x[i] = floor(z)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    46
 *			z    = (z-x[i])*2**24
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    47
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    48
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    49
 *	y[]	ouput result in an array of double precision numbers.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    50
 *		The dimension of y[] is:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    51
 *			24-bit  precision	1
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    52
 *			53-bit  precision	2
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    53
 *			64-bit  precision	2
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    54
 *			113-bit precision	3
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    55
 *		The actual value is the sum of them. Thus for 113-bit
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    56
 *		precison, one may have to do something like:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    57
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    58
 *		long double t,w,r_head, r_tail;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    59
 *		t = (long double)y[2] + (long double)y[1];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    60
 *		w = (long double)y[0];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    61
 *		r_head = t+w;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    62
 *		r_tail = w - (r_head - t);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    63
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    64
 *	e0	The exponent of x[0]
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    65
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    66
 *	nx	dimension of x[]
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    67
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    68
 *  	prec	an integer indicating the precision:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    69
 *			0	24  bits (single)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    70
 *			1	53  bits (double)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    71
 *			2	64  bits (extended)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    72
 *			3	113 bits (quad)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    73
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    74
 *	ipio2[]
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    75
 *		integer array, contains the (24*i)-th to (24*i+23)-th
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    76
 *		bit of 2/pi after binary point. The corresponding
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    77
 *		floating value is
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    78
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    79
 *			ipio2[i] * 2^(-24(i+1)).
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    80
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    81
 * External function:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    82
 *	double scalbn(), floor();
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    83
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    84
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    85
 * Here is the description of some local variables:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    86
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    87
 * 	jk	jk+1 is the initial number of terms of ipio2[] needed
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    88
 *		in the computation. The recommended value is 2,3,4,
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    89
 *		6 for single, double, extended,and quad.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    90
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    91
 * 	jz	local integer variable indicating the number of
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    92
 *		terms of ipio2[] used.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    93
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    94
 *	jx	nx - 1
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    95
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    96
 *	jv	index for pointing to the suitable ipio2[] for the
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    97
 *		computation. In general, we want
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    98
 *			( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
    99
 *		is an integer. Thus
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   100
 *			e0-3-24*jv >= 0 or (e0-3)/24 >= jv
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   101
 *		Hence jv = max(0,(e0-3)/24).
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   102
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   103
 *	jp	jp+1 is the number of terms in PIo2[] needed, jp = jk.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   104
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   105
 * 	q[]	double array with integral value, representing the
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   106
 *		24-bits chunk of the product of x and 2/pi.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   107
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   108
 *	q0	the corresponding exponent of q[0]. Note that the
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   109
 *		exponent for q[i] would be q0-24*i.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   110
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   111
 *	PIo2[]	double precision array, obtained by cutting pi/2
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   112
 *		into 24 bits chunks.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   113
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   114
 *	f[]	ipio2[] in floating point
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   115
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   116
 *	iq[]	integer array by breaking up q[] in 24-bits chunk.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   117
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   118
 *	fq[]	final product of x*(2/pi) in fq[0],..,fq[jk]
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   119
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   120
 *	ih	integer. If >0 it indicates q[] is >= 0.5, hence
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   121
 *		it also indicates the *sign* of the result.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   122
 *
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   123
 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   124
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   125
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   126
/*
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   127
 * Constants:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   128
 * The hexadecimal values are the intended ones for the following
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   129
 * constants. The decimal values may be used, provided that the
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   130
 * compiler will convert from decimal to binary accurately enough
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   131
 * to produce the hexadecimal values shown.
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   132
 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   133
6044
35448a5ea044 Lots of fixes importing SDL source wholesale into a new iOS project
Sam Lantinga <slouken@libsdl.org>
parents: 3162
diff changeset
   134
#include "math_libm.h"
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   135
#include "math_private.h"
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   136
8670
0c15c8a2f8c3 Tossed in some SDL_asserts to make static analyzer happier.
Ryan C. Gordon <icculus@icculus.org>
parents: 6044
diff changeset
   137
#include "SDL_assert.h"
0c15c8a2f8c3 Tossed in some SDL_asserts to make static analyzer happier.
Ryan C. Gordon <icculus@icculus.org>
parents: 6044
diff changeset
   138
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   139
libm_hidden_proto(scalbn)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   140
    libm_hidden_proto(floor)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   141
#ifdef __STDC__
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   142
     static const int init_jk[] = { 2, 3, 4, 6 };       /* initial value for jk */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   143
#else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   144
     static int init_jk[] = { 2, 3, 4, 6 };
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   145
#endif
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   146
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   147
#ifdef __STDC__
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   148
static const double PIo2[] = {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   149
#else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   150
static double PIo2[] = {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   151
#endif
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   152
    1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   153
    7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   154
    5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   155
    3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   156
    1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   157
    1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   158
    2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   159
    2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   160
};
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   161
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   162
#ifdef __STDC__
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   163
static const double
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   164
#else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   165
static double
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   166
#endif
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   167
  zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07,    /* 0x41700000, 0x00000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   168
    twon24 = 5.96046447753906250000e-08;        /* 0x3E700000, 0x00000000 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   169
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   170
#ifdef __STDC__
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   171
int attribute_hidden
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   172
__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   173
                  const int32_t * ipio2)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   174
#else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   175
int attribute_hidden
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   176
__kernel_rem_pio2(x, y, e0, nx, prec, ipio2)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   177
     double x[], y[];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   178
     int e0, nx, prec;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   179
     int32_t ipio2[];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   180
#endif
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   181
{
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   182
    int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   183
    double z, fw, f[20], fq[20], q[20];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   184
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   185
    /* initialize jk */
8670
0c15c8a2f8c3 Tossed in some SDL_asserts to make static analyzer happier.
Ryan C. Gordon <icculus@icculus.org>
parents: 6044
diff changeset
   186
    SDL_assert((prec >= 0) && (prec < SDL_arraysize(init_jk)));
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   187
    jk = init_jk[prec];
8670
0c15c8a2f8c3 Tossed in some SDL_asserts to make static analyzer happier.
Ryan C. Gordon <icculus@icculus.org>
parents: 6044
diff changeset
   188
    SDL_assert((jk >= 2) && (jk <= 6));
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   189
    jp = jk;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   190
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   191
    /* determine jx,jv,q0, note that 3>q0 */
8670
0c15c8a2f8c3 Tossed in some SDL_asserts to make static analyzer happier.
Ryan C. Gordon <icculus@icculus.org>
parents: 6044
diff changeset
   192
    SDL_assert(nx > 0);
2757
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   193
    jx = nx - 1;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   194
    jv = (e0 - 3) / 24;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   195
    if (jv < 0)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   196
        jv = 0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   197
    q0 = e0 - 24 * (jv + 1);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   198
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   199
    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   200
    j = jv - jx;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   201
    m = jx + jk;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   202
    for (i = 0; i <= m; i++, j++)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   203
        f[i] = (j < 0) ? zero : (double) ipio2[j];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   204
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   205
    /* compute q[0],q[1],...q[jk] */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   206
    for (i = 0; i <= jk; i++) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   207
        for (j = 0, fw = 0.0; j <= jx; j++)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   208
            fw += x[j] * f[jx + i - j];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   209
        q[i] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   210
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   211
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   212
    jz = jk;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   213
  recompute:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   214
    /* distill q[] into iq[] reversingly */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   215
    for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   216
        fw = (double) ((int32_t) (twon24 * z));
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   217
        iq[i] = (int32_t) (z - two24 * fw);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   218
        z = q[j - 1] + fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   219
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   220
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   221
    /* compute n */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   222
    z = scalbn(z, q0);          /* actual value of z */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   223
    z -= 8.0 * floor(z * 0.125);        /* trim off integer >= 8 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   224
    n = (int32_t) z;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   225
    z -= (double) n;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   226
    ih = 0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   227
    if (q0 > 0) {               /* need iq[jz-1] to determine n */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   228
        i = (iq[jz - 1] >> (24 - q0));
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   229
        n += i;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   230
        iq[jz - 1] -= i << (24 - q0);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   231
        ih = iq[jz - 1] >> (23 - q0);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   232
    } else if (q0 == 0)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   233
        ih = iq[jz - 1] >> 23;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   234
    else if (z >= 0.5)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   235
        ih = 2;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   236
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   237
    if (ih > 0) {               /* q > 0.5 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   238
        n += 1;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   239
        carry = 0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   240
        for (i = 0; i < jz; i++) {      /* compute 1-q */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   241
            j = iq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   242
            if (carry == 0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   243
                if (j != 0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   244
                    carry = 1;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   245
                    iq[i] = 0x1000000 - j;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   246
                }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   247
            } else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   248
                iq[i] = 0xffffff - j;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   249
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   250
        if (q0 > 0) {           /* rare case: chance is 1 in 12 */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   251
            switch (q0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   252
            case 1:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   253
                iq[jz - 1] &= 0x7fffff;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   254
                break;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   255
            case 2:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   256
                iq[jz - 1] &= 0x3fffff;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   257
                break;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   258
            }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   259
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   260
        if (ih == 2) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   261
            z = one - z;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   262
            if (carry != 0)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   263
                z -= scalbn(one, q0);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   264
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   265
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   266
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   267
    /* check if recomputation is needed */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   268
    if (z == zero) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   269
        j = 0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   270
        for (i = jz - 1; i >= jk; i--)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   271
            j |= iq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   272
        if (j == 0) {           /* need recomputation */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   273
            for (k = 1; iq[jk - k] == 0; k++);  /* k = no. of terms needed */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   274
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   275
            for (i = jz + 1; i <= jz + k; i++) {        /* add q[jz+1] to q[jz+k] */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   276
                f[jx + i] = (double) ipio2[jv + i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   277
                for (j = 0, fw = 0.0; j <= jx; j++)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   278
                    fw += x[j] * f[jx + i - j];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   279
                q[i] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   280
            }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   281
            jz += k;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   282
            goto recompute;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   283
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   284
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   285
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   286
    /* chop off zero terms */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   287
    if (z == 0.0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   288
        jz -= 1;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   289
        q0 -= 24;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   290
        while (iq[jz] == 0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   291
            jz--;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   292
            q0 -= 24;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   293
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   294
    } else {                    /* break z into 24-bit if necessary */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   295
        z = scalbn(z, -q0);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   296
        if (z >= two24) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   297
            fw = (double) ((int32_t) (twon24 * z));
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   298
            iq[jz] = (int32_t) (z - two24 * fw);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   299
            jz += 1;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   300
            q0 += 24;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   301
            iq[jz] = (int32_t) fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   302
        } else
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   303
            iq[jz] = (int32_t) z;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   304
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   305
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   306
    /* convert integer "bit" chunk to floating-point value */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   307
    fw = scalbn(one, q0);
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   308
    for (i = jz; i >= 0; i--) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   309
        q[i] = fw * (double) iq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   310
        fw *= twon24;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   311
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   312
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   313
    /* compute PIo2[0,...,jp]*q[jz,...,0] */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   314
    for (i = jz; i >= 0; i--) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   315
        for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   316
            fw += PIo2[k] * q[i + k];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   317
        fq[jz - i] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   318
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   319
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   320
    /* compress fq[] into y[] */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   321
    switch (prec) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   322
    case 0:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   323
        fw = 0.0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   324
        for (i = jz; i >= 0; i--)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   325
            fw += fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   326
        y[0] = (ih == 0) ? fw : -fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   327
        break;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   328
    case 1:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   329
    case 2:
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   330
        fw = 0.0;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   331
        for (i = jz; i >= 0; i--)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   332
            fw += fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   333
        y[0] = (ih == 0) ? fw : -fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   334
        fw = fq[0] - fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   335
        for (i = 1; i <= jz; i++)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   336
            fw += fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   337
        y[1] = (ih == 0) ? fw : -fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   338
        break;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   339
    case 3:                    /* painful */
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   340
        for (i = jz; i > 0; i--) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   341
            fw = fq[i - 1] + fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   342
            fq[i] += fq[i - 1] - fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   343
            fq[i - 1] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   344
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   345
        for (i = jz; i > 1; i--) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   346
            fw = fq[i - 1] + fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   347
            fq[i] += fq[i - 1] - fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   348
            fq[i - 1] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   349
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   350
        for (fw = 0.0, i = jz; i >= 2; i--)
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   351
            fw += fq[i];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   352
        if (ih == 0) {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   353
            y[0] = fq[0];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   354
            y[1] = fq[1];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   355
            y[2] = fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   356
        } else {
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   357
            y[0] = -fq[0];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   358
            y[1] = -fq[1];
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   359
            y[2] = -fw;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   360
        }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   361
    }
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   362
    return n & 7;
0581f49c9294 Whoops, missed a file...
Sam Lantinga <slouken@libsdl.org>
parents:
diff changeset
   363
}