1*b1cdbd2cSJim Jagielski /**************************************************************
2*b1cdbd2cSJim Jagielski  *
3*b1cdbd2cSJim Jagielski  * Licensed to the Apache Software Foundation (ASF) under one
4*b1cdbd2cSJim Jagielski  * or more contributor license agreements.  See the NOTICE file
5*b1cdbd2cSJim Jagielski  * distributed with this work for additional information
6*b1cdbd2cSJim Jagielski  * regarding copyright ownership.  The ASF licenses this file
7*b1cdbd2cSJim Jagielski  * to you under the Apache License, Version 2.0 (the
8*b1cdbd2cSJim Jagielski  * "License"); you may not use this file except in compliance
9*b1cdbd2cSJim Jagielski  * with the License.  You may obtain a copy of the License at
10*b1cdbd2cSJim Jagielski  *
11*b1cdbd2cSJim Jagielski  *   http://www.apache.org/licenses/LICENSE-2.0
12*b1cdbd2cSJim Jagielski  *
13*b1cdbd2cSJim Jagielski  * Unless required by applicable law or agreed to in writing,
14*b1cdbd2cSJim Jagielski  * software distributed under the License is distributed on an
15*b1cdbd2cSJim Jagielski  * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16*b1cdbd2cSJim Jagielski  * KIND, either express or implied.  See the License for the
17*b1cdbd2cSJim Jagielski  * specific language governing permissions and limitations
18*b1cdbd2cSJim Jagielski  * under the License.
19*b1cdbd2cSJim Jagielski  *
20*b1cdbd2cSJim Jagielski  *************************************************************/
21*b1cdbd2cSJim Jagielski 
22*b1cdbd2cSJim Jagielski 
23*b1cdbd2cSJim Jagielski 
24*b1cdbd2cSJim Jagielski // MARKER(update_precomp.py): autogen include statement, do not remove
25*b1cdbd2cSJim Jagielski #include "precompiled_svtools.hxx"
26*b1cdbd2cSJim Jagielski 
27*b1cdbd2cSJim Jagielski #include <mcvmath.hxx>
28*b1cdbd2cSJim Jagielski 
29*b1cdbd2cSJim Jagielski // ---------------------------------------------------------------------
30*b1cdbd2cSJim Jagielski // die folgenden Tabellen enthalten     sin(phi) * 2**14
31*b1cdbd2cSJim Jagielski // fuer phi= 360Grad*2**-32 bis 360 Grad
32*b1cdbd2cSJim Jagielski // def. fuer x: phi=360Grad * 2**(x-16)
33*b1cdbd2cSJim Jagielski //           d.h. x =  16 -> 360 Grad
34*b1cdbd2cSJim Jagielski //                x = -16 -> (2**-16) * 360 Grad
35*b1cdbd2cSJim Jagielski //         x:         -16 ... 0 ... 15
36*b1cdbd2cSJim Jagielski //x=    0,     1,     2,     3,     4,     5,     6,      7,
37*b1cdbd2cSJim Jagielski //      8,     9,    10,    11,    12,    13,    14,     15
38*b1cdbd2cSJim Jagielski 
39*b1cdbd2cSJim Jagielski static const short CosTab[16] =
40*b1cdbd2cSJim Jagielski {
41*b1cdbd2cSJim Jagielski 	16384, 16384, 16384, 16384, 16384, 16384, 16384,  16383,
42*b1cdbd2cSJim Jagielski 	16379, 16364, 16305, 16069, 15137, 11585,     0, -16383
43*b1cdbd2cSJim Jagielski };
44*b1cdbd2cSJim Jagielski static const short SinTab[16]=
45*b1cdbd2cSJim Jagielski {
46*b1cdbd2cSJim Jagielski 		2,     3,      6,    13,    25,     50,   101,  201,
47*b1cdbd2cSJim Jagielski 	  402,   804,   1606,  3196,  6270,  11585, 16384,    0
48*b1cdbd2cSJim Jagielski };
49*b1cdbd2cSJim Jagielski 
50*b1cdbd2cSJim Jagielski /**************************************************************************
51*b1cdbd2cSJim Jagielski |*
52*b1cdbd2cSJim Jagielski |*    ImpMultBig2()
53*b1cdbd2cSJim Jagielski |*
54*b1cdbd2cSJim Jagielski |*    Beschreibung       Multiplikation fuer FixPoint-Berechnungen
55*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
56*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
57*b1cdbd2cSJim Jagielski |*
58*b1cdbd2cSJim Jagielski **************************************************************************/
59*b1cdbd2cSJim Jagielski 
60*b1cdbd2cSJim Jagielski //  first parameter should be the bigger one
61*b1cdbd2cSJim Jagielski 
ImpMultBig2(const Fix & a,const Fix & b)62*b1cdbd2cSJim Jagielski Fix ImpMultBig2( const Fix& a, const Fix& b )
63*b1cdbd2cSJim Jagielski {
64*b1cdbd2cSJim Jagielski 	Fix f;
65*b1cdbd2cSJim Jagielski 	f.x = (((b.x+FIX_A2)>>FIX_P2)*a.x+FIX_A3)>>FIX_P3;
66*b1cdbd2cSJim Jagielski 	return f;
67*b1cdbd2cSJim Jagielski }
68*b1cdbd2cSJim Jagielski 
69*b1cdbd2cSJim Jagielski /**************************************************************************
70*b1cdbd2cSJim Jagielski |*
71*b1cdbd2cSJim Jagielski |*    ImpMultBig2()
72*b1cdbd2cSJim Jagielski |*
73*b1cdbd2cSJim Jagielski |*    Beschreibung       Multiplikation fuer FixPoint-Berechnungen
74*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
75*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
76*b1cdbd2cSJim Jagielski |*
77*b1cdbd2cSJim Jagielski **************************************************************************/
78*b1cdbd2cSJim Jagielski 
79*b1cdbd2cSJim Jagielski //  first parameter should be the bigger one
80*b1cdbd2cSJim Jagielski 
ImpMultBig2(const FixCpx & ra,const FixCpx & rb)81*b1cdbd2cSJim Jagielski FixCpx ImpMultBig2( const FixCpx& ra, const FixCpx& rb )
82*b1cdbd2cSJim Jagielski {
83*b1cdbd2cSJim Jagielski 	Fix rr = ImpMultBig2(ra.r,rb.r)-ImpMultBig2(ra.i,rb.i);
84*b1cdbd2cSJim Jagielski 	Fix ii = ImpMultBig2(ra.r,rb.i)+ImpMultBig2(ra.i,rb.r);
85*b1cdbd2cSJim Jagielski 	return FixCpx( rr,ii );
86*b1cdbd2cSJim Jagielski }
87*b1cdbd2cSJim Jagielski 
88*b1cdbd2cSJim Jagielski /**************************************************************************
89*b1cdbd2cSJim Jagielski |*
90*b1cdbd2cSJim Jagielski |*    ImpSqrt()
91*b1cdbd2cSJim Jagielski |*
92*b1cdbd2cSJim Jagielski |*    Beschreibung       Wurzelfunktion fuer FixPoint-Berechnungen
93*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
94*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
95*b1cdbd2cSJim Jagielski |*
96*b1cdbd2cSJim Jagielski **************************************************************************/
97*b1cdbd2cSJim Jagielski 
ImpSqrt(sal_uLong nRadi)98*b1cdbd2cSJim Jagielski sal_uInt16 ImpSqrt( sal_uLong nRadi )
99*b1cdbd2cSJim Jagielski {
100*b1cdbd2cSJim Jagielski 	register sal_uLong  inf = 1;
101*b1cdbd2cSJim Jagielski 	register sal_uLong  sup = nRadi;
102*b1cdbd2cSJim Jagielski 	register sal_uLong sqr;
103*b1cdbd2cSJim Jagielski 
104*b1cdbd2cSJim Jagielski 	if ( !nRadi )
105*b1cdbd2cSJim Jagielski 		return 0;
106*b1cdbd2cSJim Jagielski 
107*b1cdbd2cSJim Jagielski 	while ( (inf<<1) <= sup )
108*b1cdbd2cSJim Jagielski 	{
109*b1cdbd2cSJim Jagielski 		sup >>= 1;
110*b1cdbd2cSJim Jagielski 		inf <<= 1;
111*b1cdbd2cSJim Jagielski 	}
112*b1cdbd2cSJim Jagielski 	sqr = (sup+inf) >> 1;               // Anfangswert der Iteration
113*b1cdbd2cSJim Jagielski 
114*b1cdbd2cSJim Jagielski 	sqr = (nRadi/sqr + sqr) >> 1;       // 2 Newton-Iterationen reichen fuer
115*b1cdbd2cSJim Jagielski 	sqr = (nRadi/sqr + sqr) >> 1;       // +- 1 Digit
116*b1cdbd2cSJim Jagielski 
117*b1cdbd2cSJim Jagielski 	return sal::static_int_cast< sal_uInt16 >(sqr);
118*b1cdbd2cSJim Jagielski }
119*b1cdbd2cSJim Jagielski 
120*b1cdbd2cSJim Jagielski /**************************************************************************
121*b1cdbd2cSJim Jagielski |*
122*b1cdbd2cSJim Jagielski |*    ImpExPI()
123*b1cdbd2cSJim Jagielski |*
124*b1cdbd2cSJim Jagielski |*    Beschreibung       EXPI-Funktion fuer FixPoint-Berechnungen
125*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
126*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
127*b1cdbd2cSJim Jagielski |*
128*b1cdbd2cSJim Jagielski **************************************************************************/
129*b1cdbd2cSJim Jagielski 
130*b1cdbd2cSJim Jagielski // e**(i*nPhi), Einheit nPhi: 2**16 == 360 Grad
131*b1cdbd2cSJim Jagielski 
ImpExPI(sal_uInt16 nPhi)132*b1cdbd2cSJim Jagielski FixCpx ImpExPI( sal_uInt16 nPhi )
133*b1cdbd2cSJim Jagielski {
134*b1cdbd2cSJim Jagielski 	short i;
135*b1cdbd2cSJim Jagielski 	FixCpx aIter(1L);                   // e**(0*i)
136*b1cdbd2cSJim Jagielski 	FixCpx Mul;
137*b1cdbd2cSJim Jagielski 	const char Sft=14-FIX_POST;
138*b1cdbd2cSJim Jagielski 
139*b1cdbd2cSJim Jagielski 	for ( i = 15; i >= 0; i-- )
140*b1cdbd2cSJim Jagielski 	{
141*b1cdbd2cSJim Jagielski 		if ( (1L<<i) & nPhi )
142*b1cdbd2cSJim Jagielski 		{
143*b1cdbd2cSJim Jagielski 			Mul.r.x = CosTab[i]>>Sft;   // e**(i(phi1+phi2)) =
144*b1cdbd2cSJim Jagielski 			Mul.i.x = SinTab[i]>>Sft;   // e**(i*phi1)) * e**(i*phi2))
145*b1cdbd2cSJim Jagielski 			aIter  *= Mul;
146*b1cdbd2cSJim Jagielski 		}
147*b1cdbd2cSJim Jagielski 	}
148*b1cdbd2cSJim Jagielski 
149*b1cdbd2cSJim Jagielski 	return aIter;
150*b1cdbd2cSJim Jagielski }
151*b1cdbd2cSJim Jagielski 
152*b1cdbd2cSJim Jagielski /**************************************************************************
153*b1cdbd2cSJim Jagielski |*
154*b1cdbd2cSJim Jagielski |*    ImpATanx2()
155*b1cdbd2cSJim Jagielski |*
156*b1cdbd2cSJim Jagielski |*    Beschreibung       ATANX2-Funktion fuer FixPoint-Berechnungen
157*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
158*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
159*b1cdbd2cSJim Jagielski |*
160*b1cdbd2cSJim Jagielski **************************************************************************/
161*b1cdbd2cSJim Jagielski 
162*b1cdbd2cSJim Jagielski // use for x*x+y*y==1 only
163*b1cdbd2cSJim Jagielski 
ImpATanx2(const Fix & rX,const Fix & rY)164*b1cdbd2cSJim Jagielski static sal_uInt16 ImpATanx2( const Fix& rX, const Fix& rY )
165*b1cdbd2cSJim Jagielski {
166*b1cdbd2cSJim Jagielski 	sal_uInt16      phi0 = 0;           // result angel higher part
167*b1cdbd2cSJim Jagielski 	sal_uInt16      phi = 0;            // dito lower part
168*b1cdbd2cSJim Jagielski 	long        x = rX.x;
169*b1cdbd2cSJim Jagielski 	long        y = rY.x;
170*b1cdbd2cSJim Jagielski 	long        z;
171*b1cdbd2cSJim Jagielski 	const char  Sft=14-FIX_POST;
172*b1cdbd2cSJim Jagielski 	short       i;
173*b1cdbd2cSJim Jagielski 	FixCpx      aTry;
174*b1cdbd2cSJim Jagielski 	FixCpx      aInc;
175*b1cdbd2cSJim Jagielski 	FixCpx      aIter(1L);
176*b1cdbd2cSJim Jagielski 	sal_Bool        Small = sal_False;
177*b1cdbd2cSJim Jagielski 
178*b1cdbd2cSJim Jagielski 	if ( (x==0) && (y==0) )
179*b1cdbd2cSJim Jagielski 		return 0;
180*b1cdbd2cSJim Jagielski 
181*b1cdbd2cSJim Jagielski 	if ( y < 0)
182*b1cdbd2cSJim Jagielski 	{
183*b1cdbd2cSJim Jagielski 		// reduce 3. to 1. quadrant (0..90 Degree)
184*b1cdbd2cSJim Jagielski 		phi0 += 180L * 65536L / 360L;
185*b1cdbd2cSJim Jagielski 		// turn 180 degree
186*b1cdbd2cSJim Jagielski 		y    *= -1;
187*b1cdbd2cSJim Jagielski 		x    *= -1;
188*b1cdbd2cSJim Jagielski 	}
189*b1cdbd2cSJim Jagielski 
190*b1cdbd2cSJim Jagielski 	if ( x < 0)
191*b1cdbd2cSJim Jagielski 	{
192*b1cdbd2cSJim Jagielski 		// 2. to 1. q.
193*b1cdbd2cSJim Jagielski 		phi0 += 90L * 65536L / 360L;
194*b1cdbd2cSJim Jagielski 		// turn 90 degree clockwise
195*b1cdbd2cSJim Jagielski 		z = y;
196*b1cdbd2cSJim Jagielski 		y = -x;
197*b1cdbd2cSJim Jagielski 		x = z;
198*b1cdbd2cSJim Jagielski 	}
199*b1cdbd2cSJim Jagielski 
200*b1cdbd2cSJim Jagielski 	for ( i = 13; i >= 0; i-- )
201*b1cdbd2cSJim Jagielski 	{
202*b1cdbd2cSJim Jagielski 		aInc.r.x = CosTab[i]>>Sft; // e**(i(phi1+phi2)) =
203*b1cdbd2cSJim Jagielski 		aInc.i.x = SinTab[i]>>Sft; // e**(i*phi1)) * e**(i*phi2))
204*b1cdbd2cSJim Jagielski 		aTry     = aIter*aInc;
205*b1cdbd2cSJim Jagielski 
206*b1cdbd2cSJim Jagielski 		if ( Small )
207*b1cdbd2cSJim Jagielski 		{
208*b1cdbd2cSJim Jagielski 			// is try ok
209*b1cdbd2cSJim Jagielski 		   if ( aTry.r.x >= x )
210*b1cdbd2cSJim Jagielski 		   {
211*b1cdbd2cSJim Jagielski 				aIter =  aTry;
212*b1cdbd2cSJim Jagielski 				phi   += (1<<i);
213*b1cdbd2cSJim Jagielski 			}
214*b1cdbd2cSJim Jagielski 		}
215*b1cdbd2cSJim Jagielski 		else
216*b1cdbd2cSJim Jagielski 		{
217*b1cdbd2cSJim Jagielski 			// is try ok
218*b1cdbd2cSJim Jagielski 			if ( aTry.i.x <= y )
219*b1cdbd2cSJim Jagielski 			{
220*b1cdbd2cSJim Jagielski 				aIter = aTry;
221*b1cdbd2cSJim Jagielski 				phi  += (1<<i);
222*b1cdbd2cSJim Jagielski 
223*b1cdbd2cSJim Jagielski 				if ( i > 11 )
224*b1cdbd2cSJim Jagielski 					Small=sal_True;
225*b1cdbd2cSJim Jagielski 			}
226*b1cdbd2cSJim Jagielski 		}
227*b1cdbd2cSJim Jagielski 	}
228*b1cdbd2cSJim Jagielski 
229*b1cdbd2cSJim Jagielski 	return phi0+phi;
230*b1cdbd2cSJim Jagielski }
231*b1cdbd2cSJim Jagielski 
232*b1cdbd2cSJim Jagielski /**************************************************************************
233*b1cdbd2cSJim Jagielski |*
234*b1cdbd2cSJim Jagielski |*    ImpATan2()
235*b1cdbd2cSJim Jagielski |*
236*b1cdbd2cSJim Jagielski |*    Beschreibung       ATAN-Funktion fuer FixPoint-Berechnungen
237*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
238*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
239*b1cdbd2cSJim Jagielski |*
240*b1cdbd2cSJim Jagielski **************************************************************************/
241*b1cdbd2cSJim Jagielski 
ImpATan2(const short x,const short y)242*b1cdbd2cSJim Jagielski sal_uInt16 ImpATan2( const short x, const short y )
243*b1cdbd2cSJim Jagielski {
244*b1cdbd2cSJim Jagielski 	Fix rRad = ImpSqrt(sal_uLong(long(x)*x+long(y)*y));
245*b1cdbd2cSJim Jagielski 
246*b1cdbd2cSJim Jagielski 	if ( !rRad.x )
247*b1cdbd2cSJim Jagielski 		return 0;
248*b1cdbd2cSJim Jagielski 	Fix fx = x;
249*b1cdbd2cSJim Jagielski 	fx.DivBig( rRad );            // Normiere auf Einheitskreis
250*b1cdbd2cSJim Jagielski 	Fix fy = y;
251*b1cdbd2cSJim Jagielski 	fy.DivBig( rRad );
252*b1cdbd2cSJim Jagielski 
253*b1cdbd2cSJim Jagielski 	return ImpATanx2( fx, fy );
254*b1cdbd2cSJim Jagielski }
255*b1cdbd2cSJim Jagielski 
256*b1cdbd2cSJim Jagielski /**************************************************************************
257*b1cdbd2cSJim Jagielski |*
258*b1cdbd2cSJim Jagielski |*    ImpCartToPolar()
259*b1cdbd2cSJim Jagielski |*
260*b1cdbd2cSJim Jagielski |*    Beschreibung       Koordinaaten-Wandlung
261*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
262*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
263*b1cdbd2cSJim Jagielski |*
264*b1cdbd2cSJim Jagielski **************************************************************************/
265*b1cdbd2cSJim Jagielski 
ImpCartToPolar(const short x,const short y,Fix & rRad,sal_uInt16 & rPhi)266*b1cdbd2cSJim Jagielski void ImpCartToPolar( const short x, const short y, Fix& rRad, sal_uInt16& rPhi )
267*b1cdbd2cSJim Jagielski {
268*b1cdbd2cSJim Jagielski 	rRad = Fix( ImpSqrt( sal_uLong( long(x)*x+long(y)*y ) ) );
269*b1cdbd2cSJim Jagielski 
270*b1cdbd2cSJim Jagielski 	if ( !rRad.x )
271*b1cdbd2cSJim Jagielski 		rPhi=0;
272*b1cdbd2cSJim Jagielski 	else
273*b1cdbd2cSJim Jagielski 	{
274*b1cdbd2cSJim Jagielski 		// Normiere auf Einheitskreis
275*b1cdbd2cSJim Jagielski 		Fix fx = x;
276*b1cdbd2cSJim Jagielski 		fx.DivBig(rRad);
277*b1cdbd2cSJim Jagielski 		Fix fy = y;
278*b1cdbd2cSJim Jagielski 		fy.DivBig(rRad);
279*b1cdbd2cSJim Jagielski 		rPhi = ImpATanx2(fx, fy);
280*b1cdbd2cSJim Jagielski 	}
281*b1cdbd2cSJim Jagielski }
282*b1cdbd2cSJim Jagielski 
283*b1cdbd2cSJim Jagielski /**************************************************************************
284*b1cdbd2cSJim Jagielski |*
285*b1cdbd2cSJim Jagielski |*    ImpPolarToCart()
286*b1cdbd2cSJim Jagielski |*
287*b1cdbd2cSJim Jagielski |*    Beschreibung       Koordinaaten-Wandlung
288*b1cdbd2cSJim Jagielski |*    Ersterstellung     SH 01.07.93
289*b1cdbd2cSJim Jagielski |*    Letzte Aenderung   SH 01.07.93
290*b1cdbd2cSJim Jagielski |*
291*b1cdbd2cSJim Jagielski **************************************************************************/
292*b1cdbd2cSJim Jagielski 
ImpPolarToCart(const Fix & rR,const sal_uInt16 Phi,short & rX,short & rY)293*b1cdbd2cSJim Jagielski void ImpPolarToCart( const Fix& rR, const sal_uInt16 Phi, short& rX, short& rY )
294*b1cdbd2cSJim Jagielski {
295*b1cdbd2cSJim Jagielski 	FixCpx fc = ImpExPI( Phi );  // calculate sin() & cos()
296*b1cdbd2cSJim Jagielski 	fc.GetReal().MultBig( rR );
297*b1cdbd2cSJim Jagielski 	rX = sal::static_int_cast< short >(long( fc.GetReal() ));
298*b1cdbd2cSJim Jagielski 	fc.GetImag().MultBig( rR );
299*b1cdbd2cSJim Jagielski 	rY = sal::static_int_cast< short >(long( fc.GetImag() ));
300*b1cdbd2cSJim Jagielski }
301*b1cdbd2cSJim Jagielski 
302