    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.complex;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    022    import org.apache.commons.math3.exception.MathIllegalStateException;
    023    import org.apache.commons.math3.exception.OutOfRangeException;
    024    import org.apache.commons.math3.exception.ZeroException;
    025    import org.apache.commons.math3.exception.util.LocalizedFormats;
    026    import org.apache.commons.math3.util.FastMath;
    027
    028    /**
    029     * A helper class for the computation and caching of the {@code n}-th roots of
    030     * unity.
    031     *
    032     * @version $Id: RootsOfUnity.java 1416643 2012-12-03 19:37:14Z tn $
    033     * @since 3.0
    034     */
    035    public class RootsOfUnity implements Serializable {
    036
    037        /** Serializable version id. */
    038        private static final long serialVersionUID = 20120201L;
    039
    040        /** Number of roots of unity. */
    041        private int omegaCount;
    042
    043        /** Real part of the roots. */
    044        private double[] omegaReal;
    045
    046        /**
    047         * Imaginary part of the {@code n}-th roots of unity, for positive values
    048         * of {@code n}. In this array, the roots are stored in counter-clockwise
    049         * order.
    050         */
    051        private double[] omegaImaginaryCounterClockwise;
    052
    053        /**
    054         * Imaginary part of the {@code n}-th roots of unity, for negative values
    055         * of {@code n}. In this array, the roots are stored in clockwise order.
    056         */
    057        private double[] omegaImaginaryClockwise;
    058
    059        /**
    060         * {@code true} if {@link #computeRoots(int)} was called with a positive
    061         * value of its argument {@code n}. In this case, counter-clockwise ordering
    062         * of the roots of unity should be used.
    063         */
    064        private boolean isCounterClockWise;
    065
    066        /**
    067         * Build an engine for computing the {@code n}-th roots of unity.
    068         */
    069        public RootsOfUnity() {
    070
    071            omegaCount = 0;
    072            omegaReal = null;
    073            omegaImaginaryCounterClockwise = null;
    074            omegaImaginaryClockwise = null;
    075            isCounterClockWise = true;
    076        }
    077
    078        /**
    079         * Returns {@code true} if {@link #computeRoots(int)} was called with a
    080         * positive value of its argument {@code n}. If {@code true}, then
    081         * counter-clockwise ordering of the roots of unity should be used.
    082         *
    083         * @return {@code true} if the roots of unity are stored in
    084         * counter-clockwise order
    085         * @throws MathIllegalStateException if no roots of unity have been computed
    086         * yet
    087         */
    088        public synchronized boolean isCounterClockWise()
    089                throws MathIllegalStateException {
    090
    091            if (omegaCount == 0) {
    092                throw new MathIllegalStateException(
    093                        LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
    094            }
    095            return isCounterClockWise;
    096        }
    097
    098        /**
    099         * <p>
    100         * Computes the {@code n}-th roots of unity. The roots are stored in
    101         * {@code omega[]}, such that {@code omega[k] = w ^ k}, where
    102         * {@code k = 0, ..., n - 1}, {@code w = exp(2 * pi * i / n)} and
    103         * {@code i = sqrt(-1)}.
    104         * </p>
    105         * <p>
    106         * Note that {@code n} can be positive of negative
    107         * </p>
    108         * <ul>
    109         * <li>{@code abs(n)} is always the number of roots of unity.</li>
    110         * <li>If {@code n > 0}, then the roots are stored in counter-clockwise order.</li>
    111         * <li>If {@code n < 0}, then the roots are stored in clockwise order.</p>
    112         * </ul>
    113         *
    114         * @param n the (signed) number of roots of unity to be computed
    115         * @throws ZeroException if {@code n = 0}
    116         */
    117        public synchronized void computeRoots(int n) throws ZeroException {
    118
    119            if (n == 0) {
    120                throw new ZeroException(
    121                        LocalizedFormats.CANNOT_COMPUTE_0TH_ROOT_OF_UNITY);
    122            }
    123
    124            isCounterClockWise = n > 0;
    125
    126            // avoid repetitive calculations
    127            final int absN = FastMath.abs(n);
    128
    129            if (absN == omegaCount) {
    130                return;
    131            }
    132
    133            // calculate everything from scratch
    134            final double t = 2.0 * FastMath.PI / absN;
    135            final double cosT = FastMath.cos(t);
    136            final double sinT = FastMath.sin(t);
    137            omegaReal = new double[absN];
    138            omegaImaginaryCounterClockwise = new double[absN];
    139            omegaImaginaryClockwise = new double[absN];
    140            omegaReal[0] = 1.0;
    141            omegaImaginaryCounterClockwise[0] = 0.0;
    142            omegaImaginaryClockwise[0] = 0.0;
    143            for (int i = 1; i < absN; i++) {
    144                omegaReal[i] = omegaReal[i - 1] * cosT -
    145                        omegaImaginaryCounterClockwise[i - 1] * sinT;
    146                omegaImaginaryCounterClockwise[i] = omegaReal[i - 1] * sinT +
    147                        omegaImaginaryCounterClockwise[i - 1] * cosT;
    148                omegaImaginaryClockwise[i] = -omegaImaginaryCounterClockwise[i];
    149            }
    150            omegaCount = absN;
    151        }
    152
    153        /**
    154         * Get the real part of the {@code k}-th {@code n}-th root of unity.
    155         *
    156         * @param k index of the {@code n}-th root of unity
    157         * @return real part of the {@code k}-th {@code n}-th root of unity
    158         * @throws MathIllegalStateException if no roots of unity have been
    159         * computed yet
    160         * @throws MathIllegalArgumentException if {@code k} is out of range
    161         */
    162        public synchronized double getReal(int k)
    163                throws MathIllegalStateException, MathIllegalArgumentException {
    164
    165            if (omegaCount == 0) {
    166                throw new MathIllegalStateException(
    167                        LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
    168            }
    169            if ((k < 0) || (k >= omegaCount)) {
    170                throw new OutOfRangeException(
    171                        LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX,
    172                        Integer.valueOf(k),
    173                        Integer.valueOf(0),
    174                        Integer.valueOf(omegaCount - 1));
    175            }
    176
    177            return omegaReal[k];
    178        }
    179
    180        /**
    181         * Get the imaginary part of the {@code k}-th {@code n}-th root of unity.
    182         *
    183         * @param k index of the {@code n}-th root of unity
    184         * @return imaginary part of the {@code k}-th {@code n}-th root of unity
    185         * @throws MathIllegalStateException if no roots of unity have been
    186         * computed yet
    187         * @throws OutOfRangeException if {@code k} is out of range
    188         */
    189        public synchronized double getImaginary(int k)
    190                throws MathIllegalStateException, OutOfRangeException {
    191
    192            if (omegaCount == 0) {
    193                throw new MathIllegalStateException(
    194                        LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
    195            }
    196            if ((k < 0) || (k >= omegaCount)) {
    197                throw new OutOfRangeException(
    198                        LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX,
    199                        Integer.valueOf(k),
    200                        Integer.valueOf(0),
    201                        Integer.valueOf(omegaCount - 1));
    202            }
    203
    204            return isCounterClockWise ? omegaImaginaryCounterClockwise[k] :
    205                omegaImaginaryClockwise[k];
    206        }
    207
    208        /**
    209         * Returns the number of roots of unity currently stored. If
    210         * {@link #computeRoots(int)} was called with {@code n}, then this method
    211         * returns {@code abs(n)}. If no roots of unity have been computed yet, this
    212         * method returns 0.
    213         *
    214         * @return the number of roots of unity currently stored
    215         */
    216        public synchronized int getNumberOfRoots() {
    217            return omegaCount;
    218        }
    219    }