125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
3025c28e8Piotr Jasiukajtis#ifdef __RESTRICT
3125c28e8Piotr Jasiukajtis#define restrict _Restrict
3225c28e8Piotr Jasiukajtis#else
3325c28e8Piotr Jasiukajtis#define restrict
3425c28e8Piotr Jasiukajtis#endif
3525c28e8Piotr Jasiukajtis
3625c28e8Piotr Jasiukajtis/* float powf(float x, float y)
3725c28e8Piotr Jasiukajtis *
3825c28e8Piotr Jasiukajtis * Method :
3925c28e8Piotr Jasiukajtis *	1. Special cases:
4025c28e8Piotr Jasiukajtis *	for (anything) ** 0					=> 1
4125c28e8Piotr Jasiukajtis *		for (anything) ** NaN				=> QNaN + invalid
4225c28e8Piotr Jasiukajtis *		for NaN ** (anything)				=> QNaN + invalid
4325c28e8Piotr Jasiukajtis *		for +-1 ** +-Inf				=> QNaN + invalid
4425c28e8Piotr Jasiukajtis *		for +-(|x| < 1) ** +Inf				=> +0
4525c28e8Piotr Jasiukajtis *		for +-(|x| < 1) ** -Inf				=> +Inf
4625c28e8Piotr Jasiukajtis *		for +-(|x| > 1) ** +Inf				=> +Inf
4725c28e8Piotr Jasiukajtis *		for +-(|x| > 1) ** -Inf				=> +0
4825c28e8Piotr Jasiukajtis *		for +Inf ** (negative)				=> +0
4925c28e8Piotr Jasiukajtis *		for +Inf ** (positive)				=> +Inf
5025c28e8Piotr Jasiukajtis *		for -Inf ** (negative except odd integer)	=> +0
5125c28e8Piotr Jasiukajtis *		for -Inf ** (negative odd integer)		=> -0
5225c28e8Piotr Jasiukajtis *		for -Inf ** (positive except odd integer)	=> +Inf
5325c28e8Piotr Jasiukajtis *		for -Inf ** (positive odd integer)		=> -Inf
5425c28e8Piotr Jasiukajtis *		for (negative) ** (non-integer)			=> QNaN + invalid
5525c28e8Piotr Jasiukajtis *		for +0 ** (negative)				=> +Inf + overflow
5625c28e8Piotr Jasiukajtis *		for +0 ** (positive)				=> +0
5725c28e8Piotr Jasiukajtis *		for -0 ** (negative except odd integer)		=> +Inf + overflow
5825c28e8Piotr Jasiukajtis *		for -0 ** (negative odd integer)		=> -Inf + overflow
5925c28e8Piotr Jasiukajtis *		for -0 ** (positive except odd integer)		=> +0
6025c28e8Piotr Jasiukajtis *		for -0 ** (positive odd integer)		=> -0
6125c28e8Piotr Jasiukajtis *	2. Computes x**y from:
6225c28e8Piotr Jasiukajtis *		x**y = 2**(y*log2(x)) = 2**(w/256), where w = 256*log2(|x|)*y.
6325c28e8Piotr Jasiukajtis *	3. Computes w = 256 * log2(|x|) * y from
6425c28e8Piotr Jasiukajtis *		|x| = m * 2**n => log2(|x|) = n + log2(m).
6525c28e8Piotr Jasiukajtis *	Let m = m0 + dm, where m0 = 1 + k / 128,
6625c28e8Piotr Jasiukajtis *		k = [0, 128],
6725c28e8Piotr Jasiukajtis *		dm = [-1/256, 1/256].
6825c28e8Piotr Jasiukajtis *	Then 256*log2(m) = 256*log2(m0 + dm) = 256*log2(m0) + 256*log2(1+z),
6925c28e8Piotr Jasiukajtis *		where z = dm*(1/m0), z = [-1/258, 1/256].
7025c28e8Piotr Jasiukajtis *	Then
7125c28e8Piotr Jasiukajtis *		1/m0 is looked up in a table of 1, 1/(1+1/128), ..., 1/(1+128/128).
7225c28e8Piotr Jasiukajtis *		256*log2(m0) is looked up in a table of 256*log2(1), 256*log2(1+1/128),
7325c28e8Piotr Jasiukajtis *			..., 256*log2(1+128/128).
7425c28e8Piotr Jasiukajtis *		256*log2(1+z) is computed using approximation:
7525c28e8Piotr Jasiukajtis *			256*log2(1+z) = (((a3*z + a2)*z + a1)*z + a0)*z.
7625c28e8Piotr Jasiukajtis *	3. For w >= 32768
7725c28e8Piotr Jasiukajtis *		then for (negative) ** (odd integer)		=> -Inf + overflow
7825c28e8Piotr Jasiukajtis *		else						=> +Inf + overflow
7925c28e8Piotr Jasiukajtis *	For w <= -38400
8025c28e8Piotr Jasiukajtis *		then for (negative) ** (odd integer)		=> -0 + underflow
8125c28e8Piotr Jasiukajtis *		else						=> +0 + underflow
8225c28e8Piotr Jasiukajtis *	4. Computes 2 ** (w/256) from:
8325c28e8Piotr Jasiukajtis *		2 ** (w/256) = 2**a  *  2**(k/256)  *  2**(r/256)
8425c28e8Piotr Jasiukajtis *	Where:
8525c28e8Piotr Jasiukajtis *		a    =    int  ( w ) >> 8;
8625c28e8Piotr Jasiukajtis *		k    =    int  ( w ) & 0xFF;
8725c28e8Piotr Jasiukajtis *		r    =    frac ( w ).
8825c28e8Piotr Jasiukajtis *	Note that:
8925c28e8Piotr Jasiukajtis *		k = 0, 1, ..., 255;
9025c28e8Piotr Jasiukajtis *		r = (-1, 1).
9125c28e8Piotr Jasiukajtis *	Then:
9225c28e8Piotr Jasiukajtis *		2**(k/256) is looked up in a table of 2**0, 2**1/256, ...
9325c28e8Piotr Jasiukajtis *		2**(r/256) is computed using approximation:
9425c28e8Piotr Jasiukajtis *			2**(r/256) =  a0 + a1 * r + a2 * r**2
9525c28e8Piotr Jasiukajtis *		Multiplication by 2**a is done by adding "a" to
9625c28e8Piotr Jasiukajtis *		the biased exponent.
9725c28e8Piotr Jasiukajtis *	5. For (negative) ** (odd integer)	=> -(2**(w/256))
9825c28e8Piotr Jasiukajtis *	otherwise				=>   2**(w/256)
9925c28e8Piotr Jasiukajtis *
10025c28e8Piotr Jasiukajtis * Accuracy:
10125c28e8Piotr Jasiukajtis *	Max. relative aproximation error < 2**(-37.35) for 256*log2(1+z).
10225c28e8Piotr Jasiukajtis *	Max. relative aproximation error < 2**(-29.18) for 2**(r/256).
10325c28e8Piotr Jasiukajtis *	All calculations are done in double precision.
10425c28e8Piotr Jasiukajtis *	Maximum error observed: less than 0.528 ulp after 700.000.000
10525c28e8Piotr Jasiukajtis *	results.
10625c28e8Piotr Jasiukajtis */
10725c28e8Piotr Jasiukajtis
10825c28e8Piotr Jasiukajtisstatic void __vpowfx(int n, float * restrict px, float * restrict py,
10925c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez);
11025c28e8Piotr Jasiukajtis
11125c28e8Piotr Jasiukajtisstatic void __vpowf_n(int n, float * restrict px, int stridex, float * restrict py,
11225c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez);
11325c28e8Piotr Jasiukajtis
11425c28e8Piotr Jasiukajtisstatic void __vpowfx_n(int n, double yy, float * restrict py,
11525c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez);
11625c28e8Piotr Jasiukajtis
11725c28e8Piotr Jasiukajtis#pragma no_inline(__vpowfx)
11825c28e8Piotr Jasiukajtis#pragma no_inline(__vpowf_n)
11925c28e8Piotr Jasiukajtis#pragma no_inline(__vpowfx_n)
12025c28e8Piotr Jasiukajtis
12125c28e8Piotr Jasiukajtisstatic const double __TBL_exp2f[] = {
12225c28e8Piotr Jasiukajtis	/* 2^(i/256), i = [0, 255] */
12325c28e8Piotr Jasiukajtis1.000000000000000000e+00, 1.002711275050202522e+00, 1.005429901112802726e+00,
12425c28e8Piotr Jasiukajtis1.008155898118417548e+00, 1.010889286051700475e+00, 1.013630084951489430e+00,
12525c28e8Piotr Jasiukajtis1.016378314910953096e+00, 1.019133996077737914e+00, 1.021897148654116627e+00,
12625c28e8Piotr Jasiukajtis1.024667792897135721e+00, 1.027445949118763746e+00, 1.030231637686040980e+00,
12725c28e8Piotr Jasiukajtis1.033024879021228415e+00, 1.035825693601957198e+00, 1.038634101961378731e+00,
12825c28e8Piotr Jasiukajtis1.041450124688316103e+00, 1.044273782427413755e+00, 1.047105095879289793e+00,
12925c28e8Piotr Jasiukajtis1.049944085800687210e+00, 1.052790773004626423e+00, 1.055645178360557157e+00,
13025c28e8Piotr Jasiukajtis1.058507322794512762e+00, 1.061377227289262093e+00, 1.064254912884464499e+00,
13125c28e8Piotr Jasiukajtis1.067140400676823697e+00, 1.070033711820241873e+00, 1.072934867525975555e+00,
13225c28e8Piotr Jasiukajtis1.075843889062791048e+00, 1.078760797757119860e+00, 1.081685614993215250e+00,
13325c28e8Piotr Jasiukajtis1.084618362213309206e+00, 1.087559060917769660e+00, 1.090507732665257690e+00,
13425c28e8Piotr Jasiukajtis1.093464399072885840e+00, 1.096429081816376883e+00, 1.099401802630221914e+00,
13525c28e8Piotr Jasiukajtis1.102382583307840891e+00, 1.105371445701741173e+00, 1.108368411723678726e+00,
13625c28e8Piotr Jasiukajtis1.111373503344817548e+00, 1.114386742595892432e+00, 1.117408151567369279e+00,
13725c28e8Piotr Jasiukajtis1.120437752409606746e+00, 1.123475567333019898e+00, 1.126521618608241848e+00,
13825c28e8Piotr Jasiukajtis1.129575928566288079e+00, 1.132638519598719196e+00, 1.135709414157805464e+00,
13925c28e8Piotr Jasiukajtis1.138788634756691565e+00, 1.141876203969561576e+00, 1.144972144431804173e+00,
14025c28e8Piotr Jasiukajtis1.148076478840178938e+00, 1.151189229952982673e+00, 1.154310420590215935e+00,
14125c28e8Piotr Jasiukajtis1.157440073633751121e+00, 1.160578212027498779e+00, 1.163724858777577476e+00,
14225c28e8Piotr Jasiukajtis1.166880036952481658e+00, 1.170043769683250190e+00, 1.173216080163637320e+00,
14325c28e8Piotr Jasiukajtis1.176396991650281221e+00, 1.179586527462875845e+00, 1.182784710984341014e+00,
14425c28e8Piotr Jasiukajtis1.185991565660993841e+00, 1.189207115002721027e+00, 1.192431382583151178e+00,
14525c28e8Piotr Jasiukajtis1.195664392039827328e+00, 1.198906167074380580e+00, 1.202156731452703076e+00,
14625c28e8Piotr Jasiukajtis1.205416109005123859e+00, 1.208684323626581625e+00, 1.211961399276801243e+00,
14725c28e8Piotr Jasiukajtis1.215247359980468955e+00, 1.218542229827408452e+00, 1.221846032972757623e+00,
14825c28e8Piotr Jasiukajtis1.225158793637145527e+00, 1.228480536106870025e+00, 1.231811284734075862e+00,
14925c28e8Piotr Jasiukajtis1.235151063936933413e+00, 1.238499898199816540e+00, 1.241857812073484002e+00,
15025c28e8Piotr Jasiukajtis1.245224830175257980e+00, 1.248600977189204819e+00, 1.251986277866316222e+00,
15125c28e8Piotr Jasiukajtis1.255380757024691096e+00, 1.258784439549716527e+00, 1.262197350394250739e+00,
15225c28e8Piotr Jasiukajtis1.265619514578806282e+00, 1.269050957191733220e+00, 1.272491703389402762e+00,
15325c28e8Piotr Jasiukajtis1.275941778396392001e+00, 1.279401207505669325e+00, 1.282870016078778264e+00,
15425c28e8Piotr Jasiukajtis1.286348229546025568e+00, 1.289835873406665723e+00, 1.293332973229089466e+00,
15525c28e8Piotr Jasiukajtis1.296839554651009641e+00, 1.300355643379650594e+00, 1.303881265191935812e+00,
15625c28e8Piotr Jasiukajtis1.307416445934677318e+00, 1.310961211524764414e+00, 1.314515587949354636e+00,
15725c28e8Piotr Jasiukajtis1.318079601266064049e+00, 1.321653277603157539e+00, 1.325236643159741323e+00,
15825c28e8Piotr Jasiukajtis1.328829724205954355e+00, 1.332432547083161500e+00, 1.336045138204145832e+00,
15925c28e8Piotr Jasiukajtis1.339667524053302916e+00, 1.343299731186835322e+00, 1.346941786232945804e+00,
16025c28e8Piotr Jasiukajtis1.350593715892034474e+00, 1.354255546936892651e+00, 1.357927306212901142e+00,
16125c28e8Piotr Jasiukajtis1.361609020638224754e+00, 1.365300717204011915e+00, 1.369002422974590516e+00,
16225c28e8Piotr Jasiukajtis1.372714165087668414e+00, 1.376435970754530169e+00, 1.380167867260237990e+00,
16325c28e8Piotr Jasiukajtis1.383909881963832023e+00, 1.387662042298529075e+00, 1.391424375771926236e+00,
16425c28e8Piotr Jasiukajtis1.395196909966200272e+00, 1.398979672538311236e+00, 1.402772691220204759e+00,
16525c28e8Piotr Jasiukajtis1.406575993819015435e+00, 1.410389608217270663e+00, 1.414213562373095145e+00,
16625c28e8Piotr Jasiukajtis1.418047884320415175e+00, 1.421892602169165576e+00, 1.425747744105494208e+00,
16725c28e8Piotr Jasiukajtis1.429613338391970023e+00, 1.433489413367788901e+00, 1.437375997448982368e+00,
16825c28e8Piotr Jasiukajtis1.441273119128625657e+00, 1.445180806977046650e+00, 1.449099089642035043e+00,
16925c28e8Piotr Jasiukajtis1.453027995849052623e+00, 1.456967554401443765e+00, 1.460917794180647045e+00,
17025c28e8Piotr Jasiukajtis1.464878744146405731e+00, 1.468850433336981842e+00, 1.472832890869367528e+00,
17125c28e8Piotr Jasiukajtis1.476826145939499346e+00, 1.480830227822471867e+00, 1.484845165872752393e+00,
17225c28e8Piotr Jasiukajtis1.488870989524397004e+00, 1.492907728291264835e+00, 1.496955411767235455e+00,
17325c28e8Piotr Jasiukajtis1.501014069626425584e+00, 1.505083731623406473e+00, 1.509164427593422841e+00,
17425c28e8Piotr Jasiukajtis1.513256187452609813e+00, 1.517359041198214742e+00, 1.521473018908814590e+00,
17525c28e8Piotr Jasiukajtis1.525598150744538417e+00, 1.529734466947286986e+00, 1.533881997840955913e+00,
17625c28e8Piotr Jasiukajtis1.538040773831656827e+00, 1.542210825407940744e+00, 1.546392183141021448e+00,
17725c28e8Piotr Jasiukajtis1.550584877684999974e+00, 1.554788939777088652e+00, 1.559004400237836929e+00,
17825c28e8Piotr Jasiukajtis1.563231289971357629e+00, 1.567469639965552997e+00, 1.571719481292341403e+00,
17925c28e8Piotr Jasiukajtis1.575980845107886497e+00, 1.580253762652824578e+00, 1.584538265252493749e+00,
18025c28e8Piotr Jasiukajtis1.588834384317163950e+00, 1.593142151342266999e+00, 1.597461597908627073e+00,
18125c28e8Piotr Jasiukajtis1.601792755682693414e+00, 1.606135656416771029e+00, 1.610490331949254283e+00,
18225c28e8Piotr Jasiukajtis1.614856814204860713e+00, 1.619235135194863728e+00, 1.623625327017328868e+00,
18325c28e8Piotr Jasiukajtis1.628027421857347834e+00, 1.632441451987274972e+00, 1.636867449766964411e+00,
18425c28e8Piotr Jasiukajtis1.641305447644006321e+00, 1.645755478153964946e+00, 1.650217573920617742e+00,
18525c28e8Piotr Jasiukajtis1.654691767656194301e+00, 1.659178092161616158e+00, 1.663676580326736376e+00,
18625c28e8Piotr Jasiukajtis1.668187265130582464e+00, 1.672710179641596628e+00, 1.677245357017878469e+00,
18725c28e8Piotr Jasiukajtis1.681792830507429004e+00, 1.686352633448393368e+00, 1.690924799269305279e+00,
18825c28e8Piotr Jasiukajtis1.695509361489332623e+00, 1.700106353718523478e+00, 1.704715809658051251e+00,
18925c28e8Piotr Jasiukajtis1.709337763100462926e+00, 1.713972247929925974e+00, 1.718619298122477934e+00,
19025c28e8Piotr Jasiukajtis1.723278947746273992e+00, 1.727951230961837670e+00, 1.732636182022311067e+00,
19125c28e8Piotr Jasiukajtis1.737333835273706217e+00, 1.742044225155156445e+00, 1.746767386199169048e+00,
19225c28e8Piotr Jasiukajtis1.751503353031878207e+00, 1.756252160373299454e+00, 1.761013843037583904e+00,
19325c28e8Piotr Jasiukajtis1.765788435933272726e+00, 1.770575974063554714e+00, 1.775376492526521188e+00,
19425c28e8Piotr Jasiukajtis1.780190026515424462e+00, 1.785016611318934965e+00, 1.789856282321401038e+00,
19525c28e8Piotr Jasiukajtis1.794709075003107168e+00, 1.799575024940535117e+00, 1.804454167806623932e+00,
19625c28e8Piotr Jasiukajtis1.809346539371031959e+00, 1.814252175500398856e+00, 1.819171112158608494e+00,
19725c28e8Piotr Jasiukajtis1.824103385407053413e+00, 1.829049031404897274e+00, 1.834008086409342431e+00,
19825c28e8Piotr Jasiukajtis1.838980586775893711e+00, 1.843966568958625984e+00, 1.848966069510450838e+00,
19925c28e8Piotr Jasiukajtis1.853979125083385471e+00, 1.859005772428820480e+00, 1.864046048397788979e+00,
20025c28e8Piotr Jasiukajtis1.869099989941238604e+00, 1.874167634110299963e+00, 1.879249018056560194e+00,
20125c28e8Piotr Jasiukajtis1.884344179032334532e+00, 1.889453154390939194e+00, 1.894575981586965607e+00,
20225c28e8Piotr Jasiukajtis1.899712698176555303e+00, 1.904863341817674138e+00, 1.910027950270389852e+00,
20325c28e8Piotr Jasiukajtis1.915206561397147400e+00, 1.920399213163047403e+00, 1.925605943636125028e+00,
20425c28e8Piotr Jasiukajtis1.930826790987627106e+00, 1.936061793492294347e+00, 1.941310989528640452e+00,
20525c28e8Piotr Jasiukajtis1.946574417579233218e+00, 1.951852116230978318e+00, 1.957144124175400179e+00,
20625c28e8Piotr Jasiukajtis1.962450480208927317e+00, 1.967771223233175881e+00, 1.973106392255234320e+00,
20725c28e8Piotr Jasiukajtis1.978456026387950928e+00, 1.983820164850219392e+00, 1.989198846967266343e+00,
20825c28e8Piotr Jasiukajtis1.994592112170940235e+00
20925c28e8Piotr Jasiukajtis};
21025c28e8Piotr Jasiukajtis
21125c28e8Piotr Jasiukajtisstatic const double __TBL_log2f[] = {
21225c28e8Piotr Jasiukajtis	/* __TBL_log2f[2*i] = 256*log2(1+i/128), i = [0, 128] */
21325c28e8Piotr Jasiukajtis	/* __TBL_log2f[2*i+1] = 2**(-23)/(1+i/128), i = [0, 128] */
21425c28e8Piotr Jasiukajtis0.000000000000000000e+00, 1.192092895507812500e-07, 2.874177388353054585e+00,
21525c28e8Piotr Jasiukajtis1.182851865310077503e-07, 5.726160135284354524e+00, 1.173753004807692373e-07,
21625c28e8Piotr Jasiukajtis8.556288393587271557e+00, 1.164793058206106825e-07, 1.136489455576407970e+01,
21725c28e8Piotr Jasiukajtis1.155968868371212153e-07, 1.415230348830453799e+01, 1.147277373120300688e-07,
21825c28e8Piotr Jasiukajtis1.691883275718974389e+01, 1.138715601679104456e-07, 1.966479284501270897e+01,
21925c28e8Piotr Jasiukajtis1.130280671296296339e-07, 2.239048736008688678e+01, 1.121969784007352926e-07,
22025c28e8Piotr Jasiukajtis2.509621323789484038e+01, 1.113780223540145949e-07, 2.778226093521127638e+01,
22125c28e8Piotr Jasiukajtis1.105709352355072477e-07, 3.044891461721790193e+01, 1.097754608812949697e-07,
22225c28e8Piotr Jasiukajtis3.309645233791141550e+01, 1.089913504464285680e-07, 3.572514621409114710e+01,
22325c28e8Piotr Jasiukajtis1.082183621453900683e-07, 3.833526259319860685e+01, 1.074562610035211292e-07,
22425c28e8Piotr Jasiukajtis4.092706221526768928e+01, 1.067048186188811188e-07, 4.350080036923196758e+01,
22525c28e8Piotr Jasiukajtis1.059638129340277719e-07, 4.605672704382322280e+01, 1.052330280172413778e-07,
22625c28e8Piotr Jasiukajtis4.859508707328441091e+01, 1.045122538527397202e-07, 5.111612027810928538e+01,
22725c28e8Piotr Jasiukajtis1.038012861394557784e-07, 5.362006160101114460e+01, 1.030999260979729787e-07,
22825c28e8Piotr Jasiukajtis5.610714123831336053e+01, 1.024079802852348971e-07, 5.857758476694550609e+01,
22925c28e8Piotr Jasiukajtis1.017252604166666732e-07, 6.103161326722020164e+01, 1.010515831953642383e-07,
23025c28e8Piotr Jasiukajtis6.346944344155788542e+01, 1.003867701480263102e-07, 6.589128772931884725e+01,
23125c28e8Piotr Jasiukajtis9.973064746732026447e-08, 6.829735441789475203e+01, 9.908304586038961692e-08,
23225c28e8Piotr Jasiukajtis7.068784775020480993e+01, 9.844380040322580637e-08, 7.306296802873558249e+01,
23325c28e8Piotr Jasiukajtis9.781275040064102225e-08, 7.542291171625650748e+01, 9.718973925159236158e-08,
23425c28e8Piotr Jasiukajtis7.776787153333835079e+01, 9.657461431962025166e-08, 8.009803655279496581e+01,
23525c28e8Piotr Jasiukajtis9.596722680817610579e-08, 8.241359229116476115e+01, 9.536743164062500529e-08,
23625c28e8Piotr Jasiukajtis8.471472079734193983e+01, 9.477508734472049048e-08, 8.700160073846393516e+01,
23725c28e8Piotr Jasiukajtis9.419005594135801946e-08, 8.927440748315585495e+01, 9.361220283742331508e-08,
23825c28e8Piotr Jasiukajtis9.153331318222942059e+01, 9.304139672256097884e-08, 9.377848684692884262e+01,
23925c28e8Piotr Jasiukajtis9.247750946969696962e-08, 9.601009442481273481e+01, 9.192041603915663129e-08,
24025c28e8Piotr Jasiukajtis9.822829887335737453e+01, 9.136999438622755046e-08, 1.004332602313626381e+02,
24125c28e8Piotr Jasiukajtis9.082612537202380448e-08, 1.026251356882391832e+02, 9.028869267751479078e-08,
24225c28e8Piotr Jasiukajtis1.048040796512516550e+02, 8.975758272058823405e-08, 1.069702438107898530e+02,
24325c28e8Piotr Jasiukajtis8.923268457602338686e-08, 1.091237772037370775e+02, 8.871388989825581272e-08,
24425c28e8Piotr Jasiukajtis1.112648262750015107e+02, 8.820109284682080489e-08, 1.133935349372744383e+02,
24525c28e8Piotr Jasiukajtis8.769419001436781487e-08, 1.155100446290761766e+02, 8.719308035714285707e-08,
24625c28e8Piotr Jasiukajtis1.176144943711480977e+02, 8.669766512784091150e-08, 1.197070208212473403e+02,
24725c28e8Piotr Jasiukajtis8.620784781073446298e-08, 1.217877583273978246e+02, 8.572353405898876167e-08,
24825c28e8Piotr Jasiukajtis1.238568389796496376e+02, 8.524463163407821503e-08, 1.259143926603967287e+02,
24925c28e8Piotr Jasiukajtis8.477105034722222546e-08, 1.279605470933005762e+02, 8.430270200276242743e-08,
25025c28e8Piotr Jasiukajtis1.299954278908662388e+02, 8.383950034340659995e-08, 1.320191586007148601e+02,
25125c28e8Piotr Jasiukajtis8.338136099726775949e-08, 1.340318607505952855e+02, 8.292820142663043248e-08,
25225c28e8Piotr Jasiukajtis1.360336538921758915e+02, 8.247994087837838296e-08, 1.380246556436560468e+02,
25325c28e8Piotr Jasiukajtis8.203650033602151192e-08, 1.400049817312349774e+02, 8.159780247326202734e-08,
25425c28e8Piotr Jasiukajtis1.419747460294751704e+02, 8.116377160904255122e-08, 1.439340606005945915e+02,
25525c28e8Piotr Jasiukajtis8.073433366402115954e-08, 1.458830357327226466e+02, 8.030941611842105082e-08,
25625c28e8Piotr Jasiukajtis1.478217799771516638e+02, 7.988894797120419333e-08, 1.497504001846159838e+02,
25725c28e8Piotr Jasiukajtis7.947285970052082892e-08, 1.516690015406285852e+02, 7.906108322538860398e-08,
25825c28e8Piotr Jasiukajtis1.535776875999046922e+02, 7.865355186855669953e-08, 1.554765603199003294e+02,
25925c28e8Piotr Jasiukajtis7.825020032051282044e-08, 1.573657200934933087e+02, 7.785096460459183052e-08,
26025c28e8Piotr Jasiukajtis1.592452657808323124e+02, 7.745578204314720208e-08, 1.611152947403800511e+02,
26125c28e8Piotr Jasiukajtis7.706459122474748130e-08, 1.629759028591741128e+02, 7.667733197236181018e-08,
26225c28e8Piotr Jasiukajtis1.648271845823295223e+02, 7.629394531250000159e-08, 1.666692329418057170e+02,
26325c28e8Piotr Jasiukajtis7.591437344527363039e-08, 1.685021395844594565e+02, 7.553855971534653557e-08,
26425c28e8Piotr Jasiukajtis1.703259947994051231e+02, 7.516644858374384321e-08, 1.721408875447028777e+02,
26525c28e8Piotr Jasiukajtis7.479798560049019504e-08, 1.739469054733941960e+02, 7.443311737804878042e-08,
26625c28e8Piotr Jasiukajtis1.757441349589039135e+02, 7.407179156553397416e-08, 1.775326611198272531e+02,
26725c28e8Piotr Jasiukajtis7.371395682367149407e-08, 1.793125678441195987e+02, 7.335956280048077330e-08,
26825c28e8Piotr Jasiukajtis1.810839378127059831e+02, 7.300856010765549954e-08, 1.828468525225273993e+02,
26925c28e8Piotr Jasiukajtis7.266090029761905417e-08, 1.846013923090393973e+02, 7.231653584123223301e-08,
27025c28e8Piotr Jasiukajtis1.863476363681789962e+02, 7.197542010613207272e-08, 1.880856627778145764e+02,
27125c28e8Piotr Jasiukajtis7.163750733568075279e-08, 1.898155485186936176e+02, 7.130275262850466758e-08,
27225c28e8Piotr Jasiukajtis1.915373694949018386e+02, 7.097111191860465018e-08, 1.932512005538479514e+02,
27325c28e8Piotr Jasiukajtis7.064254195601851460e-08, 1.949571155057867031e+02, 7.031700028801843312e-08,
27425c28e8Piotr Jasiukajtis1.966551871428931406e+02, 6.999444524082569196e-08, 1.983454872579004018e+02,
27525c28e8Piotr Jasiukajtis6.967483590182648015e-08, 2.000280866623128588e+02, 6.935813210227272390e-08,
27625c28e8Piotr Jasiukajtis2.017030552042064926e+02, 6.904429440045249486e-08, 2.033704617856271284e+02,
27725c28e8Piotr Jasiukajtis6.873328406531531472e-08, 2.050303743795980154e+02, 6.842506306053811558e-08,
27825c28e8Piotr Jasiukajtis2.066828600467466401e+02, 6.811959402901785336e-08, 2.083279849515614899e+02,
27925c28e8Piotr Jasiukajtis6.781684027777777772e-08, 2.099658143782880586e+02, 6.751676576327433535e-08,
28025c28e8Piotr Jasiukajtis2.115964127464742432e+02, 6.721933507709251725e-08, 2.132198436261738550e+02,
28125c28e8Piotr Jasiukajtis6.692451343201754014e-08, 2.148361697528176535e+02, 6.663226664847161225e-08,
28225c28e8Piotr Jasiukajtis2.164454530417600608e+02, 6.634256114130434863e-08, 2.180477546025107358e+02,
28325c28e8Piotr Jasiukajtis6.605536390692640687e-08, 2.196431347526584545e+02, 6.577064251077586116e-08,
28425c28e8Piotr Jasiukajtis2.212316530314957390e+02, 6.548836507510729591e-08, 2.228133682133515663e+02,
28525c28e8Piotr Jasiukajtis6.520850026709402365e-08, 2.243883383206399174e+02, 6.493101728723404362e-08,
28625c28e8Piotr Jasiukajtis2.259566206366313565e+02, 6.465588585805084723e-08, 2.275182717179543204e+02,
28725c28e8Piotr Jasiukajtis6.438307621308016336e-08, 2.290733474068335340e+02, 6.411255908613445100e-08,
28825c28e8Piotr Jasiukajtis2.306219028430716378e+02, 6.384430570083681460e-08, 2.321639924757807307e+02,
28925c28e8Piotr Jasiukajtis6.357828776041666578e-08, 2.336996700748701699e+02, 6.331447743775933615e-08,
29025c28e8Piotr Jasiukajtis2.352289887422961954e+02, 6.305284736570248109e-08, 2.367520009230799189e+02,
29125c28e8Piotr Jasiukajtis6.279337062757202180e-08, 2.382687584160988763e+02, 6.253602074795082293e-08,
29225c28e8Piotr Jasiukajtis2.397793123846580556e+02, 6.228077168367347501e-08, 2.412837133668454044e+02,
29325c28e8Piotr Jasiukajtis6.202759781504065697e-08, 2.427820112856774699e+02, 6.177647393724696421e-08,
29425c28e8Piotr Jasiukajtis2.442742554590400630e+02, 6.152737525201612732e-08, 2.457604946094287186e+02,
29525c28e8Piotr Jasiukajtis6.128027735943774537e-08, 2.472407768734942692e+02, 6.103515625000000127e-08,
29625c28e8Piotr Jasiukajtis2.487151498113976231e+02, 6.079198829681274795e-08, 2.501836604159786077e+02,
29725c28e8Piotr Jasiukajtis6.055075024801586965e-08, 2.516463551217433974e+02, 6.031141921936758485e-08,
29825c28e8Piotr Jasiukajtis2.531032798136744475e+02, 6.007397268700787318e-08, 2.545544798358676246e+02,
29925c28e8Piotr Jasiukajtis5.983838848039215603e-08, 2.560000000000000000e+02, 5.960464477539062500e-08
30025c28e8Piotr Jasiukajtis};
30125c28e8Piotr Jasiukajtis
30225c28e8Piotr Jasiukajtisstatic const double __TBL_expfb[] = {
30325c28e8Piotr Jasiukajtis7.006492321624085355e-46, 1.401298464324817071e-45, 2.802596928649634142e-45,
30425c28e8Piotr Jasiukajtis5.605193857299268284e-45, 1.121038771459853657e-44, 2.242077542919707313e-44,
30525c28e8Piotr Jasiukajtis4.484155085839414627e-44, 8.968310171678829254e-44, 1.793662034335765851e-43,
30625c28e8Piotr Jasiukajtis3.587324068671531702e-43, 7.174648137343063403e-43, 1.434929627468612681e-42,
30725c28e8Piotr Jasiukajtis2.869859254937225361e-42, 5.739718509874450723e-42, 1.147943701974890145e-41,
30825c28e8Piotr Jasiukajtis2.295887403949780289e-41, 4.591774807899560578e-41, 9.183549615799121156e-41,
30925c28e8Piotr Jasiukajtis1.836709923159824231e-40, 3.673419846319648462e-40, 7.346839692639296925e-40,
31025c28e8Piotr Jasiukajtis1.469367938527859385e-39, 2.938735877055718770e-39, 5.877471754111437540e-39,
31125c28e8Piotr Jasiukajtis1.175494350822287508e-38, 2.350988701644575016e-38, 4.701977403289150032e-38,
31225c28e8Piotr Jasiukajtis9.403954806578300064e-38, 1.880790961315660013e-37, 3.761581922631320025e-37,
31325c28e8Piotr Jasiukajtis7.523163845262640051e-37, 1.504632769052528010e-36, 3.009265538105056020e-36,
31425c28e8Piotr Jasiukajtis6.018531076210112041e-36, 1.203706215242022408e-35, 2.407412430484044816e-35,
31525c28e8Piotr Jasiukajtis4.814824860968089633e-35, 9.629649721936179265e-35, 1.925929944387235853e-34,
31625c28e8Piotr Jasiukajtis3.851859888774471706e-34, 7.703719777548943412e-34, 1.540743955509788682e-33,
31725c28e8Piotr Jasiukajtis3.081487911019577365e-33, 6.162975822039154730e-33, 1.232595164407830946e-32,
31825c28e8Piotr Jasiukajtis2.465190328815661892e-32, 4.930380657631323784e-32, 9.860761315262647568e-32,
31925c28e8Piotr Jasiukajtis1.972152263052529514e-31, 3.944304526105059027e-31, 7.888609052210118054e-31,
32025c28e8Piotr Jasiukajtis1.577721810442023611e-30, 3.155443620884047222e-30, 6.310887241768094443e-30,
32125c28e8Piotr Jasiukajtis1.262177448353618889e-29, 2.524354896707237777e-29, 5.048709793414475555e-29,
32225c28e8Piotr Jasiukajtis1.009741958682895111e-28, 2.019483917365790222e-28, 4.038967834731580444e-28,
32325c28e8Piotr Jasiukajtis8.077935669463160887e-28, 1.615587133892632177e-27, 3.231174267785264355e-27,
32425c28e8Piotr Jasiukajtis6.462348535570528710e-27, 1.292469707114105742e-26, 2.584939414228211484e-26,
32525c28e8Piotr Jasiukajtis5.169878828456422968e-26, 1.033975765691284594e-25, 2.067951531382569187e-25,
32625c28e8Piotr Jasiukajtis4.135903062765138374e-25, 8.271806125530276749e-25, 1.654361225106055350e-24,
32725c28e8Piotr Jasiukajtis3.308722450212110699e-24, 6.617444900424221399e-24, 1.323488980084844280e-23,
32825c28e8Piotr Jasiukajtis2.646977960169688560e-23, 5.293955920339377119e-23, 1.058791184067875424e-22,
32925c28e8Piotr Jasiukajtis2.117582368135750848e-22, 4.235164736271501695e-22, 8.470329472543003391e-22,
33025c28e8Piotr Jasiukajtis1.694065894508600678e-21, 3.388131789017201356e-21, 6.776263578034402713e-21,
33125c28e8Piotr Jasiukajtis1.355252715606880543e-20, 2.710505431213761085e-20, 5.421010862427522170e-20,
33225c28e8Piotr Jasiukajtis1.084202172485504434e-19, 2.168404344971008868e-19, 4.336808689942017736e-19,
33325c28e8Piotr Jasiukajtis8.673617379884035472e-19, 1.734723475976807094e-18, 3.469446951953614189e-18,
33425c28e8Piotr Jasiukajtis6.938893903907228378e-18, 1.387778780781445676e-17, 2.775557561562891351e-17,
33525c28e8Piotr Jasiukajtis5.551115123125782702e-17, 1.110223024625156540e-16, 2.220446049250313081e-16,
33625c28e8Piotr Jasiukajtis4.440892098500626162e-16, 8.881784197001252323e-16, 1.776356839400250465e-15,
33725c28e8Piotr Jasiukajtis3.552713678800500929e-15, 7.105427357601001859e-15, 1.421085471520200372e-14,
33825c28e8Piotr Jasiukajtis2.842170943040400743e-14, 5.684341886080801487e-14, 1.136868377216160297e-13,
33925c28e8Piotr Jasiukajtis2.273736754432320595e-13, 4.547473508864641190e-13, 9.094947017729282379e-13,
34025c28e8Piotr Jasiukajtis1.818989403545856476e-12, 3.637978807091712952e-12, 7.275957614183425903e-12,
34125c28e8Piotr Jasiukajtis1.455191522836685181e-11, 2.910383045673370361e-11, 5.820766091346740723e-11,
34225c28e8Piotr Jasiukajtis1.164153218269348145e-10, 2.328306436538696289e-10, 4.656612873077392578e-10,
34325c28e8Piotr Jasiukajtis9.313225746154785156e-10, 1.862645149230957031e-09, 3.725290298461914062e-09,
34425c28e8Piotr Jasiukajtis7.450580596923828125e-09, 1.490116119384765625e-08, 2.980232238769531250e-08,
34525c28e8Piotr Jasiukajtis5.960464477539062500e-08, 1.192092895507812500e-07, 2.384185791015625000e-07,
34625c28e8Piotr Jasiukajtis4.768371582031250000e-07, 9.536743164062500000e-07, 1.907348632812500000e-06,
34725c28e8Piotr Jasiukajtis3.814697265625000000e-06, 7.629394531250000000e-06, 1.525878906250000000e-05,
34825c28e8Piotr Jasiukajtis3.051757812500000000e-05, 6.103515625000000000e-05, 1.220703125000000000e-04,
34925c28e8Piotr Jasiukajtis2.441406250000000000e-04, 4.882812500000000000e-04, 9.765625000000000000e-04,
35025c28e8Piotr Jasiukajtis1.953125000000000000e-03, 3.906250000000000000e-03, 7.812500000000000000e-03,
35125c28e8Piotr Jasiukajtis1.562500000000000000e-02, 3.125000000000000000e-02, 6.250000000000000000e-02,
35225c28e8Piotr Jasiukajtis1.250000000000000000e-01, 2.500000000000000000e-01, 5.000000000000000000e-01,
35325c28e8Piotr Jasiukajtis1.000000000000000000e+00, 2.000000000000000000e+00, 4.000000000000000000e+00,
35425c28e8Piotr Jasiukajtis8.000000000000000000e+00, 1.600000000000000000e+01, 3.200000000000000000e+01,
35525c28e8Piotr Jasiukajtis6.400000000000000000e+01, 1.280000000000000000e+02, 2.560000000000000000e+02,
35625c28e8Piotr Jasiukajtis5.120000000000000000e+02, 1.024000000000000000e+03, 2.048000000000000000e+03,
35725c28e8Piotr Jasiukajtis4.096000000000000000e+03, 8.192000000000000000e+03, 1.638400000000000000e+04,
35825c28e8Piotr Jasiukajtis3.276800000000000000e+04, 6.553600000000000000e+04, 1.310720000000000000e+05,
35925c28e8Piotr Jasiukajtis2.621440000000000000e+05, 5.242880000000000000e+05, 1.048576000000000000e+06,
36025c28e8Piotr Jasiukajtis2.097152000000000000e+06, 4.194304000000000000e+06, 8.388608000000000000e+06,
36125c28e8Piotr Jasiukajtis1.677721600000000000e+07, 3.355443200000000000e+07, 6.710886400000000000e+07,
36225c28e8Piotr Jasiukajtis1.342177280000000000e+08, 2.684354560000000000e+08, 5.368709120000000000e+08,
36325c28e8Piotr Jasiukajtis1.073741824000000000e+09, 2.147483648000000000e+09, 4.294967296000000000e+09,
36425c28e8Piotr Jasiukajtis8.589934592000000000e+09, 1.717986918400000000e+10, 3.435973836800000000e+10,
36525c28e8Piotr Jasiukajtis6.871947673600000000e+10, 1.374389534720000000e+11, 2.748779069440000000e+11,
36625c28e8Piotr Jasiukajtis5.497558138880000000e+11, 1.099511627776000000e+12, 2.199023255552000000e+12,
36725c28e8Piotr Jasiukajtis4.398046511104000000e+12, 8.796093022208000000e+12, 1.759218604441600000e+13,
36825c28e8Piotr Jasiukajtis3.518437208883200000e+13, 7.036874417766400000e+13, 1.407374883553280000e+14,
36925c28e8Piotr Jasiukajtis2.814749767106560000e+14, 5.629499534213120000e+14, 1.125899906842624000e+15,
37025c28e8Piotr Jasiukajtis2.251799813685248000e+15, 4.503599627370496000e+15, 9.007199254740992000e+15,
37125c28e8Piotr Jasiukajtis1.801439850948198400e+16, 3.602879701896396800e+16, 7.205759403792793600e+16,
37225c28e8Piotr Jasiukajtis1.441151880758558720e+17, 2.882303761517117440e+17, 5.764607523034234880e+17,
37325c28e8Piotr Jasiukajtis1.152921504606846976e+18, 2.305843009213693952e+18, 4.611686018427387904e+18,
37425c28e8Piotr Jasiukajtis9.223372036854775808e+18, 1.844674407370955162e+19, 3.689348814741910323e+19,
37525c28e8Piotr Jasiukajtis7.378697629483820646e+19, 1.475739525896764129e+20, 2.951479051793528259e+20,
37625c28e8Piotr Jasiukajtis5.902958103587056517e+20, 1.180591620717411303e+21, 2.361183241434822607e+21,
37725c28e8Piotr Jasiukajtis4.722366482869645214e+21, 9.444732965739290427e+21, 1.888946593147858085e+22,
37825c28e8Piotr Jasiukajtis3.777893186295716171e+22, 7.555786372591432342e+22, 1.511157274518286468e+23,
37925c28e8Piotr Jasiukajtis3.022314549036572937e+23, 6.044629098073145874e+23, 1.208925819614629175e+24,
38025c28e8Piotr Jasiukajtis2.417851639229258349e+24, 4.835703278458516699e+24, 9.671406556917033398e+24,
38125c28e8Piotr Jasiukajtis1.934281311383406680e+25, 3.868562622766813359e+25, 7.737125245533626718e+25,
38225c28e8Piotr Jasiukajtis1.547425049106725344e+26, 3.094850098213450687e+26, 6.189700196426901374e+26,
38325c28e8Piotr Jasiukajtis1.237940039285380275e+27, 2.475880078570760550e+27, 4.951760157141521100e+27,
38425c28e8Piotr Jasiukajtis9.903520314283042199e+27, 1.980704062856608440e+28, 3.961408125713216880e+28,
38525c28e8Piotr Jasiukajtis7.922816251426433759e+28, 1.584563250285286752e+29, 3.169126500570573504e+29,
38625c28e8Piotr Jasiukajtis6.338253001141147007e+29, 1.267650600228229401e+30, 2.535301200456458803e+30,
38725c28e8Piotr Jasiukajtis5.070602400912917606e+30, 1.014120480182583521e+31, 2.028240960365167042e+31,
38825c28e8Piotr Jasiukajtis4.056481920730334085e+31, 8.112963841460668170e+31, 1.622592768292133634e+32,
38925c28e8Piotr Jasiukajtis3.245185536584267268e+32, 6.490371073168534536e+32, 1.298074214633706907e+33,
39025c28e8Piotr Jasiukajtis2.596148429267413814e+33, 5.192296858534827629e+33, 1.038459371706965526e+34,
39125c28e8Piotr Jasiukajtis2.076918743413931051e+34, 4.153837486827862103e+34, 8.307674973655724206e+34,
39225c28e8Piotr Jasiukajtis1.661534994731144841e+35, 3.323069989462289682e+35, 6.646139978924579365e+35,
39325c28e8Piotr Jasiukajtis1.329227995784915873e+36, 2.658455991569831746e+36, 5.316911983139663492e+36,
39425c28e8Piotr Jasiukajtis1.063382396627932698e+37, 2.126764793255865397e+37, 4.253529586511730793e+37,
39525c28e8Piotr Jasiukajtis8.507059173023461587e+37, 1.701411834604692317e+38, 3.402823669209384635e+38
39625c28e8Piotr Jasiukajtis};
39725c28e8Piotr Jasiukajtis
39825c28e8Piotr Jasiukajtisstatic const double
39925c28e8Piotr Jasiukajtis	KA3 = -3.60659926599003171364e-01*256.0,
40025c28e8Piotr Jasiukajtis	KA2 =  4.80902715189356683026e-01*256.0,
40125c28e8Piotr Jasiukajtis	KA1 = -7.21347520569871841065e-01*256.0,
40225c28e8Piotr Jasiukajtis	KA0 =  1.44269504088069658645e+00*256.0,
40325c28e8Piotr Jasiukajtis	KB2 =  3.66556671660783833261e-06,
40425c28e8Piotr Jasiukajtis	KB1 =  2.70760782821392980564e-03,
40525c28e8Piotr Jasiukajtis	DONE = 1.0,
40625c28e8Piotr Jasiukajtis	HTHRESH = 32768.0,
40725c28e8Piotr Jasiukajtis	LTHRESH = -38400.0;
40825c28e8Piotr Jasiukajtis
40925c28e8Piotr Jasiukajtis#define RETURN(ret)						\
41025c28e8Piotr Jasiukajtis{								\
41125c28e8Piotr Jasiukajtis	*pz = (ret);						\
41225c28e8Piotr Jasiukajtis	px += stridex;						\
41325c28e8Piotr Jasiukajtis	py += stridey;						\
41425c28e8Piotr Jasiukajtis	pz += stridez;						\
41525c28e8Piotr Jasiukajtis	if (n_n == 0)						\
41625c28e8Piotr Jasiukajtis	{							\
41725c28e8Piotr Jasiukajtis		spx = px; spy = py; spz = pz;			\
41825c28e8Piotr Jasiukajtis		continue;					\
41925c28e8Piotr Jasiukajtis	}							\
42025c28e8Piotr Jasiukajtis	n--;							\
42125c28e8Piotr Jasiukajtis	break;							\
42225c28e8Piotr Jasiukajtis}
42325c28e8Piotr Jasiukajtis
42425c28e8Piotr Jasiukajtisvoid
42525c28e8Piotr Jasiukajtis__vpowf(int n, float * restrict px, int stridex, float * restrict py,
42625c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez)
42725c28e8Piotr Jasiukajtis{
42825c28e8Piotr Jasiukajtis	float		*spx, *spy, *spz;
42925c28e8Piotr Jasiukajtis	double		y0, yy0;
43025c28e8Piotr Jasiukajtis	long long	di0;
43125c28e8Piotr Jasiukajtis	unsigned	ux, sx, uy, ay, ax0;
43225c28e8Piotr Jasiukajtis	int		exp, i0, ind0, exp0, yisint0, n_n;
43325c28e8Piotr Jasiukajtis
43425c28e8Piotr Jasiukajtis#ifndef NOPOWFIX
43525c28e8Piotr Jasiukajtis	if (stridex == 0)
43625c28e8Piotr Jasiukajtis	{
43725c28e8Piotr Jasiukajtis		unsigned	hx = *(unsigned*)px;
43825c28e8Piotr Jasiukajtis
43925c28e8Piotr Jasiukajtis		if ((hx >= 0x00800000) &&	/* x not zero or subnormal		*/
44025c28e8Piotr Jasiukajtis		(hx < 0x7f800000) &&		/* x not inf, nan or negative sign bit	*/
44125c28e8Piotr Jasiukajtis		(hx != 0x3f800000))		/* x not 1				*/
44225c28e8Piotr Jasiukajtis		{
44325c28e8Piotr Jasiukajtis			__vpowfx(n, px, py, stridey, pz, stridez);
44425c28e8Piotr Jasiukajtis			return;
44525c28e8Piotr Jasiukajtis		}
44625c28e8Piotr Jasiukajtis	}
44725c28e8Piotr Jasiukajtis#endif
44825c28e8Piotr Jasiukajtis
44925c28e8Piotr Jasiukajtis	while (n > 0)
45025c28e8Piotr Jasiukajtis	{
45125c28e8Piotr Jasiukajtis		n_n = 0;
45225c28e8Piotr Jasiukajtis		spx = px;
45325c28e8Piotr Jasiukajtis		spy = py;
45425c28e8Piotr Jasiukajtis		spz = pz;
45525c28e8Piotr Jasiukajtis		for (; n > 0 ; n--)
45625c28e8Piotr Jasiukajtis		{
45725c28e8Piotr Jasiukajtis			uy = *(unsigned int*)py;
45825c28e8Piotr Jasiukajtis			ux = *(unsigned int*)px;
45925c28e8Piotr Jasiukajtis			ay = uy & 0x7fffffff;
46025c28e8Piotr Jasiukajtis			ax0 = ux & 0x7fffffff;
46125c28e8Piotr Jasiukajtis			sx = ux >> 31;
46225c28e8Piotr Jasiukajtis			yisint0 = 0;	/* Y - non-integer */
46325c28e8Piotr Jasiukajtis
46425c28e8Piotr Jasiukajtis			/* |X| or |Y| = Inf,Nan */
46525c28e8Piotr Jasiukajtis			if (ax0 >= 0x7f800000 || ay >= 0x7f800000)
46625c28e8Piotr Jasiukajtis			{
46725c28e8Piotr Jasiukajtis				if (ay == 0)
46825c28e8Piotr Jasiukajtis					RETURN(1.0f)	/* pow(X,0) */
46925c28e8Piotr Jasiukajtis				/* |X| or |Y| = Nan */
47025c28e8Piotr Jasiukajtis				if (ax0 > 0x7f800000 || ay > 0x7f800000)
47125c28e8Piotr Jasiukajtis					RETURN (*px + *py)
47225c28e8Piotr Jasiukajtis				if (ay == 0x7f800000)		/* |Y| = Inf */
47325c28e8Piotr Jasiukajtis				{
47425c28e8Piotr Jasiukajtis					float fy;
47525c28e8Piotr Jasiukajtis					if (ax0 == 0x3f800000)
47625c28e8Piotr Jasiukajtis						fy = *py - *py;		/* +-1 ** +-Inf = NaN */
47725c28e8Piotr Jasiukajtis					else
47825c28e8Piotr Jasiukajtis						fy = ((ax0 < 0x3f800000) != (uy >> 31)) ? 0.0f : *(float*) &ay;
47925c28e8Piotr Jasiukajtis					RETURN(fy)
48025c28e8Piotr Jasiukajtis				}
48125c28e8Piotr Jasiukajtis				if (sx)	/* X = -Inf */
48225c28e8Piotr Jasiukajtis				{
48325c28e8Piotr Jasiukajtis					exp = ay >> 23;
48425c28e8Piotr Jasiukajtis					if (exp >= 0x97)	/* |Y| >= 2^24 */
48525c28e8Piotr Jasiukajtis						yisint0 = 2;	/* Y - even */
48625c28e8Piotr Jasiukajtis					else if (exp >= 0x7f)	/* |Y| >= 1 */
48725c28e8Piotr Jasiukajtis					{
48825c28e8Piotr Jasiukajtis						i0 = ay >> ((0x7f + 23) - exp);
48925c28e8Piotr Jasiukajtis						if ((i0 << ((0x7f + 23) - exp)) == ay)
49025c28e8Piotr Jasiukajtis							yisint0 = 2 - (i0 & 1);
49125c28e8Piotr Jasiukajtis					}
49225c28e8Piotr Jasiukajtis				}
49325c28e8Piotr Jasiukajtis				if (uy >> 31)
49425c28e8Piotr Jasiukajtis					ax0 = 0;
49525c28e8Piotr Jasiukajtis				ax0 += yisint0 << 31;
49625c28e8Piotr Jasiukajtis				RETURN(*(float*)&ax0)
49725c28e8Piotr Jasiukajtis			}
49825c28e8Piotr Jasiukajtis
49925c28e8Piotr Jasiukajtis			if ((int)ux < 0x00800000)	/* X = denormal or negative */
50025c28e8Piotr Jasiukajtis			{
50125c28e8Piotr Jasiukajtis				if (ay == 0)
50225c28e8Piotr Jasiukajtis					RETURN(1.0f)	/* pow(X,0) */
50325c28e8Piotr Jasiukajtis				exp0 = (ax0 >> 23) - 127;
50425c28e8Piotr Jasiukajtis
50525c28e8Piotr Jasiukajtis				if ((int)ax0 < 0x00800000)	/* X = denormal */
50625c28e8Piotr Jasiukajtis				{
50725c28e8Piotr Jasiukajtis					*((float*) &ax0) = (float) (int)ax0;
50825c28e8Piotr Jasiukajtis					exp0 = (ax0 >> 23) - (127 + 149);
50925c28e8Piotr Jasiukajtis				}
51025c28e8Piotr Jasiukajtis
51125c28e8Piotr Jasiukajtis				if ((int)ux <= 0)	/* X <= 0 */
51225c28e8Piotr Jasiukajtis				{
51325c28e8Piotr Jasiukajtis					exp = ay >> 23;
51425c28e8Piotr Jasiukajtis					if (exp >= 0x97)	/* |Y| >= 2^24 */
51525c28e8Piotr Jasiukajtis						yisint0 = 2;	/* Y - even */
51625c28e8Piotr Jasiukajtis					else if (exp >= 0x7f)	/* |Y| >= 1 */
51725c28e8Piotr Jasiukajtis					{
51825c28e8Piotr Jasiukajtis						i0 = ay >> ((0x7f + 23) - exp);
51925c28e8Piotr Jasiukajtis						if ((i0 << ((0x7f + 23) - exp)) == ay)
52025c28e8Piotr Jasiukajtis							yisint0 = 2 - (i0 & 1);
52125c28e8Piotr Jasiukajtis					}
52225c28e8Piotr Jasiukajtis
52325c28e8Piotr Jasiukajtis					if (ax0 == 0)		/* pow(0,Y) */
52425c28e8Piotr Jasiukajtis					{
52525c28e8Piotr Jasiukajtis						float fy;
52625c28e8Piotr Jasiukajtis						fy = (uy >> 31) ? 1.0f / 0.0f : 0.0f;
52725c28e8Piotr Jasiukajtis						if (sx & yisint0)
52825c28e8Piotr Jasiukajtis							fy = -fy;
52925c28e8Piotr Jasiukajtis						RETURN(fy)
53025c28e8Piotr Jasiukajtis					}
53125c28e8Piotr Jasiukajtis
53225c28e8Piotr Jasiukajtis					if (yisint0 == 0)	/* pow(neg,non-integer) */
53325c28e8Piotr Jasiukajtis						RETURN(0.0f / 0.0f)	/* NaN */
53425c28e8Piotr Jasiukajtis				}
53525c28e8Piotr Jasiukajtis
53625c28e8Piotr Jasiukajtis				/* perform yy0 = 256*log2(xi)*yi */
53725c28e8Piotr Jasiukajtis				ax0 &= 0x007fffff;
53825c28e8Piotr Jasiukajtis				i0 = (ax0 + 0x8000) & 0xffff0000;
53925c28e8Piotr Jasiukajtis				ind0 = i0 >> 15;
54025c28e8Piotr Jasiukajtis				i0 = ax0 - i0;
54125c28e8Piotr Jasiukajtis				y0 = (double) i0 * __TBL_log2f[ind0 + 1];
54225c28e8Piotr Jasiukajtis				yy0 = __TBL_log2f[ind0] + (double) (exp0 << 8);
54325c28e8Piotr Jasiukajtis				yy0 += (((KA3 * y0 + KA2) * y0 + KA1) * y0 + KA0) * y0;
54425c28e8Piotr Jasiukajtis				yy0 = (double)py[0] * yy0;
54525c28e8Piotr Jasiukajtis
54625c28e8Piotr Jasiukajtis				/* perform 2 ** (yy0/256) */
54725c28e8Piotr Jasiukajtis				if (yy0 >= HTHRESH)
54825c28e8Piotr Jasiukajtis					yy0 = HTHRESH;
54925c28e8Piotr Jasiukajtis				if (yy0 <= LTHRESH)
55025c28e8Piotr Jasiukajtis					yy0 = LTHRESH;
55125c28e8Piotr Jasiukajtis				ind0 = (int) yy0;
55225c28e8Piotr Jasiukajtis				y0 = yy0 - (double)ind0;
55325c28e8Piotr Jasiukajtis				yy0 = (KB2 * y0 + KB1) * y0 + DONE;
55425c28e8Piotr Jasiukajtis				di0 = ((long long)((ind0 >> 8) + (yisint0 << 11))) << 52;
55525c28e8Piotr Jasiukajtis				di0 += ((long long*)__TBL_exp2f)[ind0 & 255];
55625c28e8Piotr Jasiukajtis				RETURN((float) (yy0 * *(double*)&di0))
55725c28e8Piotr Jasiukajtis			}
55825c28e8Piotr Jasiukajtis			px += stridex;
55925c28e8Piotr Jasiukajtis			py += stridey;
56025c28e8Piotr Jasiukajtis			pz += stridez;
56125c28e8Piotr Jasiukajtis			n_n++;
56225c28e8Piotr Jasiukajtis		}
56325c28e8Piotr Jasiukajtis		if (n_n > 0)
56425c28e8Piotr Jasiukajtis			__vpowf_n(n_n, spx, stridex, spy, stridey, spz, stridez);
56525c28e8Piotr Jasiukajtis	}
56625c28e8Piotr Jasiukajtis}
56725c28e8Piotr Jasiukajtis
56825c28e8Piotr Jasiukajtis
56925c28e8Piotr Jasiukajtisstatic void
57025c28e8Piotr Jasiukajtis__vpowf_n(int n, float * restrict px, int stridex, float * restrict py,
57125c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez)
57225c28e8Piotr Jasiukajtis{
57325c28e8Piotr Jasiukajtis	double		y0, yy0;
57425c28e8Piotr Jasiukajtis	double		di0;
57525c28e8Piotr Jasiukajtis	int		ind0, i0, exp0;
57625c28e8Piotr Jasiukajtis	unsigned	ax0;
57725c28e8Piotr Jasiukajtis	double		y1, yy1;
57825c28e8Piotr Jasiukajtis	double		di1;
57925c28e8Piotr Jasiukajtis	int		ind1, i1, exp1;
58025c28e8Piotr Jasiukajtis	unsigned	ax1;
58125c28e8Piotr Jasiukajtis	double		y2, yy2;
58225c28e8Piotr Jasiukajtis	double		di2;
58325c28e8Piotr Jasiukajtis	int		ind2, i2, exp2;
58425c28e8Piotr Jasiukajtis	unsigned	ax2;
58525c28e8Piotr Jasiukajtis
58625c28e8Piotr Jasiukajtis	for (; n > 2 ; n -= 3)
58725c28e8Piotr Jasiukajtis	{
58825c28e8Piotr Jasiukajtis		/* perform yy0 = 256*log2(xi)*yi */
58925c28e8Piotr Jasiukajtis		ax0 = ((int*)px)[0];
59025c28e8Piotr Jasiukajtis		px += stridex;
59125c28e8Piotr Jasiukajtis		ax1 = ((int*)px)[0];
59225c28e8Piotr Jasiukajtis		px += stridex;
59325c28e8Piotr Jasiukajtis		ax2 = ((int*)px)[0];
59425c28e8Piotr Jasiukajtis		px += stridex;
59525c28e8Piotr Jasiukajtis		exp0 = ((ax0 & 0x7fffffff) >> 23) - 127;
59625c28e8Piotr Jasiukajtis		exp1 = ((ax1 & 0x7fffffff) >> 23) - 127;
59725c28e8Piotr Jasiukajtis		exp2 = ((ax2 & 0x7fffffff) >> 23) - 127;
59825c28e8Piotr Jasiukajtis		ax0 &= 0x007fffff;
59925c28e8Piotr Jasiukajtis		ax1 &= 0x007fffff;
60025c28e8Piotr Jasiukajtis		ax2 &= 0x007fffff;
60125c28e8Piotr Jasiukajtis		i0 = (ax0 + 0x8000) & 0xffff0000;
60225c28e8Piotr Jasiukajtis		i1 = (ax1 + 0x8000) & 0xffff0000;
60325c28e8Piotr Jasiukajtis		i2 = (ax2 + 0x8000) & 0xffff0000;
60425c28e8Piotr Jasiukajtis		ind0 = i0 >> 15;
60525c28e8Piotr Jasiukajtis		ind1 = i1 >> 15;
60625c28e8Piotr Jasiukajtis		ind2 = i2 >> 15;
60725c28e8Piotr Jasiukajtis		i0 = ax0 - i0;
60825c28e8Piotr Jasiukajtis		i1 = ax1 - i1;
60925c28e8Piotr Jasiukajtis		i2 = ax2 - i2;
61025c28e8Piotr Jasiukajtis		y0 = (double) i0 * __TBL_log2f[ind0 + 1];
61125c28e8Piotr Jasiukajtis		y1 = (double) i1 * __TBL_log2f[ind1 + 1];
61225c28e8Piotr Jasiukajtis		y2 = (double) i2 * __TBL_log2f[ind2 + 1];
61325c28e8Piotr Jasiukajtis		yy0 = __TBL_log2f[ind0] + (double) (exp0 << 8);
61425c28e8Piotr Jasiukajtis		yy1 = __TBL_log2f[ind1] + (double) (exp1 << 8);
61525c28e8Piotr Jasiukajtis		yy2 = __TBL_log2f[ind2] + (double) (exp2 << 8);
61625c28e8Piotr Jasiukajtis		yy0 += (((KA3 * y0 + KA2) * y0 + KA1) * y0 + KA0) * y0;
61725c28e8Piotr Jasiukajtis		yy1 += (((KA3 * y1 + KA2) * y1 + KA1) * y1 + KA0) * y1;
61825c28e8Piotr Jasiukajtis		yy2 += (((KA3 * y2 + KA2) * y2 + KA1) * y2 + KA0) * y2;
61925c28e8Piotr Jasiukajtis		yy0 = (double)py[0] * yy0;
62025c28e8Piotr Jasiukajtis		py += stridey;
62125c28e8Piotr Jasiukajtis		yy1 = (double)py[0] * yy1;
62225c28e8Piotr Jasiukajtis		py += stridey;
62325c28e8Piotr Jasiukajtis		yy2 = (double)py[0] * yy2;
62425c28e8Piotr Jasiukajtis		py += stridey;
62525c28e8Piotr Jasiukajtis
62625c28e8Piotr Jasiukajtis		/* perform 2 ** (yy0/256) */
62725c28e8Piotr Jasiukajtis		if (yy0 >= HTHRESH)
62825c28e8Piotr Jasiukajtis			yy0 = HTHRESH;
62925c28e8Piotr Jasiukajtis		if (yy0 <= LTHRESH)
63025c28e8Piotr Jasiukajtis			yy0 = LTHRESH;
63125c28e8Piotr Jasiukajtis		if (yy1 >= HTHRESH)
63225c28e8Piotr Jasiukajtis			yy1 = HTHRESH;
63325c28e8Piotr Jasiukajtis		if (yy1 <= LTHRESH)
63425c28e8Piotr Jasiukajtis			yy1 = LTHRESH;
63525c28e8Piotr Jasiukajtis		if (yy2 >= HTHRESH)
63625c28e8Piotr Jasiukajtis			yy2 = HTHRESH;
63725c28e8Piotr Jasiukajtis		if (yy2 <= LTHRESH)
63825c28e8Piotr Jasiukajtis			yy2 = LTHRESH;
63925c28e8Piotr Jasiukajtis
64025c28e8Piotr Jasiukajtis		ind0 = (int) yy0;
64125c28e8Piotr Jasiukajtis		ind1 = (int) yy1;
64225c28e8Piotr Jasiukajtis		ind2 = (int) yy2;
64325c28e8Piotr Jasiukajtis		y0 = yy0 - (double)ind0;
64425c28e8Piotr Jasiukajtis		y1 = yy1 - (double)ind1;
64525c28e8Piotr Jasiukajtis		y2 = yy2 - (double)ind2;
64625c28e8Piotr Jasiukajtis		yy0 = (KB2 * y0 + KB1) * y0 + DONE;
64725c28e8Piotr Jasiukajtis		yy1 = (KB2 * y1 + KB1) * y1 + DONE;
64825c28e8Piotr Jasiukajtis		yy2 = (KB2 * y2 + KB1) * y2 + DONE;
64925c28e8Piotr Jasiukajtis		di0 = (__TBL_expfb + 150)[ind0 >> 8];
65025c28e8Piotr Jasiukajtis		di1 = (__TBL_expfb + 150)[ind1 >> 8];
65125c28e8Piotr Jasiukajtis		di2 = (__TBL_expfb + 150)[ind2 >> 8];
65225c28e8Piotr Jasiukajtis		di0 *= __TBL_exp2f[ind0 & 255];
65325c28e8Piotr Jasiukajtis		di1 *= __TBL_exp2f[ind1 & 255];
65425c28e8Piotr Jasiukajtis		di2 *= __TBL_exp2f[ind2 & 255];
65525c28e8Piotr Jasiukajtis		pz[0] = (float) (yy0 * di0);
65625c28e8Piotr Jasiukajtis		pz += stridez;
65725c28e8Piotr Jasiukajtis		pz[0] = (float) (yy1 * di1);
65825c28e8Piotr Jasiukajtis		pz += stridez;
65925c28e8Piotr Jasiukajtis		pz[0] = (float) (yy2 * di2);
66025c28e8Piotr Jasiukajtis		pz += stridez;
66125c28e8Piotr Jasiukajtis	}
66225c28e8Piotr Jasiukajtis
66325c28e8Piotr Jasiukajtis	for (; n > 0 ; n--)
66425c28e8Piotr Jasiukajtis	{
66525c28e8Piotr Jasiukajtis		/* perform yy0 = 256*log2(xi)*yi */
66625c28e8Piotr Jasiukajtis		ax0 = ((int*)px)[0];
66725c28e8Piotr Jasiukajtis		exp0 = ((ax0 & 0x7fffffff) >> 23) - 127;
66825c28e8Piotr Jasiukajtis		ax0 &= 0x007fffff;
66925c28e8Piotr Jasiukajtis		i0 = (ax0 + 0x8000) & 0xffff0000;
67025c28e8Piotr Jasiukajtis		ind0 = i0 >> 15;
67125c28e8Piotr Jasiukajtis		i0 = ax0 - i0;
67225c28e8Piotr Jasiukajtis		y0 = (double) i0 * __TBL_log2f[ind0 + 1];
67325c28e8Piotr Jasiukajtis		yy0 = __TBL_log2f[ind0] + (double) (exp0 << 8);
67425c28e8Piotr Jasiukajtis		yy0 += (((KA3 * y0 + KA2) * y0 + KA1) * y0 + KA0) * y0;
67525c28e8Piotr Jasiukajtis		yy0 = (double)py[0] * yy0;
67625c28e8Piotr Jasiukajtis
67725c28e8Piotr Jasiukajtis		/* perform 2 ** (yy0/256) */
67825c28e8Piotr Jasiukajtis		if (yy0 >= HTHRESH)
67925c28e8Piotr Jasiukajtis			yy0 = HTHRESH;
68025c28e8Piotr Jasiukajtis		if (yy0 <= LTHRESH)
68125c28e8Piotr Jasiukajtis			yy0 = LTHRESH;
68225c28e8Piotr Jasiukajtis		ind0 = (int) yy0;
68325c28e8Piotr Jasiukajtis		y0 = yy0 - (double)ind0;
68425c28e8Piotr Jasiukajtis		yy0 = (KB2 * y0 + KB1) * y0 + DONE;
68525c28e8Piotr Jasiukajtis		di0 = (__TBL_expfb + 150)[ind0 >> 8];
68625c28e8Piotr Jasiukajtis		di0 *= __TBL_exp2f[ind0 & 255];
68725c28e8Piotr Jasiukajtis		pz[0] = (float) (yy0 * di0);
68825c28e8Piotr Jasiukajtis		px += stridex;
68925c28e8Piotr Jasiukajtis		py += stridey;
69025c28e8Piotr Jasiukajtis		pz += stridez;
69125c28e8Piotr Jasiukajtis	}
69225c28e8Piotr Jasiukajtis}
69325c28e8Piotr Jasiukajtis
69425c28e8Piotr Jasiukajtis
69525c28e8Piotr Jasiukajtisstatic void
69625c28e8Piotr Jasiukajtis__vpowfx(int n, float * restrict px, float * restrict py,
69725c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez)
69825c28e8Piotr Jasiukajtis{
69925c28e8Piotr Jasiukajtis	float		*spy, *spz;
70025c28e8Piotr Jasiukajtis	double		yy, y0;
70125c28e8Piotr Jasiukajtis	int		ind0, exp0, i0, n_n;
70225c28e8Piotr Jasiukajtis	unsigned	ux, ax, ax0, uy, ay;
70325c28e8Piotr Jasiukajtis
70425c28e8Piotr Jasiukajtis	/* perform yy = 256*log2(xi)*yi */
70525c28e8Piotr Jasiukajtis	ux = *(unsigned int*)px;
70625c28e8Piotr Jasiukajtis	ax = ux & 0x7fffffff;
70725c28e8Piotr Jasiukajtis	exp0 = (ax >> 23) - 127;
70825c28e8Piotr Jasiukajtis	ax0 = ux & 0x007fffff;
70925c28e8Piotr Jasiukajtis	i0 = (ax0 + 0x8000) & 0xffff0000;
71025c28e8Piotr Jasiukajtis	ind0 = i0 >> 15;
71125c28e8Piotr Jasiukajtis	i0 = ax0 - i0;
71225c28e8Piotr Jasiukajtis	y0 = (double) i0 * __TBL_log2f[ind0 + 1];
71325c28e8Piotr Jasiukajtis	yy = __TBL_log2f[ind0] + (double) (exp0 << 8);
71425c28e8Piotr Jasiukajtis	yy += (((KA3 * y0 + KA2) * y0 + KA1) * y0 + KA0) * y0;
71525c28e8Piotr Jasiukajtis
71625c28e8Piotr Jasiukajtis	while (n > 0)
71725c28e8Piotr Jasiukajtis	{
71825c28e8Piotr Jasiukajtis		n_n = 0;
71925c28e8Piotr Jasiukajtis		spy = py;
72025c28e8Piotr Jasiukajtis		spz = pz;
72125c28e8Piotr Jasiukajtis		for (; n > 0 ; n--)
72225c28e8Piotr Jasiukajtis		{
72325c28e8Piotr Jasiukajtis			uy = *(unsigned int*)py;
72425c28e8Piotr Jasiukajtis			ay = uy & 0x7fffffff;
72525c28e8Piotr Jasiukajtis
72625c28e8Piotr Jasiukajtis			if (ay >= 0x7f800000)		/* |Y| = Inf or Nan */
72725c28e8Piotr Jasiukajtis			{
72825c28e8Piotr Jasiukajtis				float fy;
72925c28e8Piotr Jasiukajtis				if (ay > 0x7f800000)
73025c28e8Piotr Jasiukajtis					fy = *py + *py;	/* |Y| = Nan */
73125c28e8Piotr Jasiukajtis				else
73225c28e8Piotr Jasiukajtis					fy = ((ax < 0x3f800000) != (uy >> 31)) ? 0.0f : *(float*)&ay;
73325c28e8Piotr Jasiukajtis				*pz = fy;
73425c28e8Piotr Jasiukajtis				py += stridey;
73525c28e8Piotr Jasiukajtis				pz += stridez;
73625c28e8Piotr Jasiukajtis				if (n_n == 0)
73725c28e8Piotr Jasiukajtis				{
73825c28e8Piotr Jasiukajtis					spy = py;
73925c28e8Piotr Jasiukajtis					spz = pz;
74025c28e8Piotr Jasiukajtis					continue;
74125c28e8Piotr Jasiukajtis				}
74225c28e8Piotr Jasiukajtis				n--;
74325c28e8Piotr Jasiukajtis				break;
74425c28e8Piotr Jasiukajtis			}
74525c28e8Piotr Jasiukajtis			py += stridey;
74625c28e8Piotr Jasiukajtis			pz += stridez;
74725c28e8Piotr Jasiukajtis			n_n++;
74825c28e8Piotr Jasiukajtis		}
74925c28e8Piotr Jasiukajtis		if (n_n > 0)
75025c28e8Piotr Jasiukajtis			__vpowfx_n(n_n, yy, spy, stridey, spz, stridez);
75125c28e8Piotr Jasiukajtis	}
75225c28e8Piotr Jasiukajtis}
75325c28e8Piotr Jasiukajtis
75425c28e8Piotr Jasiukajtis
75525c28e8Piotr Jasiukajtisstatic void
75625c28e8Piotr Jasiukajtis__vpowfx_n(int n, double yy, float * restrict py,
75725c28e8Piotr Jasiukajtis	int stridey, float * restrict pz, int stridez)
75825c28e8Piotr Jasiukajtis{
75925c28e8Piotr Jasiukajtis	double		y0, yy0, di0;
76025c28e8Piotr Jasiukajtis	double		y1, yy1, di1;
76125c28e8Piotr Jasiukajtis	double		y2, yy2, di2;
76225c28e8Piotr Jasiukajtis	int		ind0, ind1, ind2;
76325c28e8Piotr Jasiukajtis
76425c28e8Piotr Jasiukajtis	for (; n > 2 ; n-= 3)
76525c28e8Piotr Jasiukajtis	{
76625c28e8Piotr Jasiukajtis		/* perform 2 ** (yy/256) */
76725c28e8Piotr Jasiukajtis		yy0 = (double)py[0] * yy;
76825c28e8Piotr Jasiukajtis		py += stridey;
76925c28e8Piotr Jasiukajtis		yy1 = (double)py[0] * yy;
77025c28e8Piotr Jasiukajtis		py += stridey;
77125c28e8Piotr Jasiukajtis		yy2 = (double)py[0] * yy;
77225c28e8Piotr Jasiukajtis		py += stridey;
77325c28e8Piotr Jasiukajtis		if (yy0 >= HTHRESH)
77425c28e8Piotr Jasiukajtis			yy0 = HTHRESH;
77525c28e8Piotr Jasiukajtis		if (yy0 <= LTHRESH)
77625c28e8Piotr Jasiukajtis			yy0 = LTHRESH;
77725c28e8Piotr Jasiukajtis		if (yy1 >= HTHRESH)
77825c28e8Piotr Jasiukajtis			yy1 = HTHRESH;
77925c28e8Piotr Jasiukajtis		if (yy1 <= LTHRESH)
78025c28e8Piotr Jasiukajtis			yy1 = LTHRESH;
78125c28e8Piotr Jasiukajtis		if (yy2 >= HTHRESH)
78225c28e8Piotr Jasiukajtis			yy2 = HTHRESH;
78325c28e8Piotr Jasiukajtis		if (yy2 <= LTHRESH)
78425c28e8Piotr Jasiukajtis			yy2 = LTHRESH;
78525c28e8Piotr Jasiukajtis		ind0 = (int) yy0;
78625c28e8Piotr Jasiukajtis		ind1 = (int) yy1;
78725c28e8Piotr Jasiukajtis		ind2 = (int) yy2;
78825c28e8Piotr Jasiukajtis		y0 = yy0 - (double)ind0;
78925c28e8Piotr Jasiukajtis		y1 = yy1 - (double)ind1;
79025c28e8Piotr Jasiukajtis		y2 = yy2 - (double)ind2;
79125c28e8Piotr Jasiukajtis		yy0 = (KB2 * y0 + KB1) * y0 + DONE;
79225c28e8Piotr Jasiukajtis		yy1 = (KB2 * y1 + KB1) * y1 + DONE;
79325c28e8Piotr Jasiukajtis		yy2 = (KB2 * y2 + KB1) * y2 + DONE;
79425c28e8Piotr Jasiukajtis		di0 = (__TBL_expfb + 150)[ind0 >> 8];
79525c28e8Piotr Jasiukajtis		di1 = (__TBL_expfb + 150)[ind1 >> 8];
79625c28e8Piotr Jasiukajtis		di2 = (__TBL_expfb + 150)[ind2 >> 8];
79725c28e8Piotr Jasiukajtis		di0 *= __TBL_exp2f[ind0 & 255];
79825c28e8Piotr Jasiukajtis		di1 *= __TBL_exp2f[ind1 & 255];
79925c28e8Piotr Jasiukajtis		di2 *= __TBL_exp2f[ind2 & 255];
80025c28e8Piotr Jasiukajtis		pz[0] = (float) (yy0 * di0);
80125c28e8Piotr Jasiukajtis		pz += stridez;
80225c28e8Piotr Jasiukajtis		pz[0] = (float) (yy1 * di1);
80325c28e8Piotr Jasiukajtis		pz += stridez;
80425c28e8Piotr Jasiukajtis		pz[0] = (float) (yy2 * di2);
80525c28e8Piotr Jasiukajtis		pz += stridez;
80625c28e8Piotr Jasiukajtis	}
80725c28e8Piotr Jasiukajtis	for (; n > 0 ; n--)
80825c28e8Piotr Jasiukajtis	{
80925c28e8Piotr Jasiukajtis		/* perform 2 ** (yy/256) */
81025c28e8Piotr Jasiukajtis		yy0 = (double)py[0] * yy;
81125c28e8Piotr Jasiukajtis		if (yy0 >= HTHRESH)
81225c28e8Piotr Jasiukajtis			yy0 = HTHRESH;
81325c28e8Piotr Jasiukajtis		if (yy0 <= LTHRESH)
81425c28e8Piotr Jasiukajtis			yy0 = LTHRESH;
81525c28e8Piotr Jasiukajtis		ind0 = (int) yy0;
81625c28e8Piotr Jasiukajtis		y0 = yy0 - (double)ind0;
81725c28e8Piotr Jasiukajtis		yy0 = (KB2 * y0 + KB1) * y0 + DONE;
81825c28e8Piotr Jasiukajtis		di0 = (__TBL_expfb + 150)[ind0 >> 8];
81925c28e8Piotr Jasiukajtis		di0 *= __TBL_exp2f[ind0 & 255];
82025c28e8Piotr Jasiukajtis		pz[0] = (float) (yy0 * di0);
82125c28e8Piotr Jasiukajtis		py += stridey;
82225c28e8Piotr Jasiukajtis		pz += stridez;
82325c28e8Piotr Jasiukajtis	}
82425c28e8Piotr Jasiukajtis}
825