/*
    * Licensed to the Hipparchus project under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      https://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.hipparchus.analysis.differentiation;
  
  import java.util.Arrays;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.Field;
  import org.hipparchus.exception.LocalizedCoreFormats;
  import org.hipparchus.exception.MathIllegalArgumentException;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.FieldSinCos;
  import org.hipparchus.util.FieldSinhCosh;
  import org.hipparchus.util.MathArrays;
  import org.hipparchus.util.MathUtils;
  
  /** Class representing both the value and the differentials of a function.
   * <p>This class is a stripped-down version of {@link FieldDerivativeStructure}
   * with {@link FieldDerivativeStructure#getOrder() derivation order} limited to one.
   * It should have less overhead than {@link FieldDerivativeStructure} in its domain.</p>
   * <p>This class is an implementation of Rall's numbers. Rall's numbers are an
   * extension to the real numbers used throughout mathematical expressions; they hold
   * the derivative together with the value of a function.</p>
   * <p>{@link FieldGradient} instances can be used directly thanks to
   * the arithmetic operators to the mathematical functions provided as
   * methods by this class (+, -, *, /, %, sin, cos ...).</p>
   * <p>Implementing complex expressions by hand using {@link Derivative}-based
   * classes (or in fact any {@link org.hipparchus.CalculusFieldElement} class) is
   * a tedious and error-prone task but has the advantage of not requiring users
   * to compute the derivatives by themselves and allowing to switch for one
   * derivative implementation to another as they all share the same filed API.</p>
   * <p>Instances of this class are guaranteed to be immutable.</p>
   * @param <T> the type of the function parameters and value
   * @see DerivativeStructure
   * @see UnivariateDerivative1
   * @see UnivariateDerivative2
   * @see Gradient
   * @see FieldDerivativeStructure
   * @see FieldUnivariateDerivative1
   * @see FieldUnivariateDerivative2
   * @since 1.7
   */
  public class FieldGradient<T extends CalculusFieldElement<T>> implements FieldDerivative1<T, FieldGradient<T>> {
  
      /** Value of the function. */
      private final T value;
  
      /** Gradient of the function. */
      private final T[] grad;
  
      /** Build an instance with values and unitialized derivatives array.
       * @param value value of the function
       * @param freeParameters number of free parameters
       */
      private FieldGradient(final T value, int freeParameters) {
          this.value = value;
          this.grad  = MathArrays.buildArray(value.getField(), freeParameters);
      }
  
      /** Build an instance with values and derivative.
       * @param value value of the function
       * @param gradient gradient of the function
       */
      @SafeVarargs
      public FieldGradient(final T value, final T... gradient) {
          this(value, gradient.length);
          System.arraycopy(gradient, 0, grad, 0, grad.length);
      }
  
      /** Build an instance from a {@link FieldDerivativeStructure}.
       * @param ds derivative structure
       * @exception MathIllegalArgumentException if {@code ds} order
       * is not 1
       */
      public FieldGradient(final FieldDerivativeStructure<T> ds) throws MathIllegalArgumentException {
          this(ds.getValue(), ds.getFreeParameters());
          MathUtils.checkDimension(ds.getOrder(), 1);
          System.arraycopy(ds.getAllDerivatives(), 1, grad, 0, grad.length);
      }
  
      /** Build an instance corresponding to a constant value.
       * @param freeParameters number of free parameters (i.e. dimension of the gradient)
       * @param value constant value of the function
       * @param <T> the type of the function parameters and value
       * @return a {@code FieldGradient} with a constant value and all derivatives set to 0.0
      */
     public static <T extends CalculusFieldElement<T>> FieldGradient<T> constant(final int freeParameters, final T value) {
         final FieldGradient<T> g = new FieldGradient<>(value, freeParameters);
         Arrays.fill(g.grad, value.getField().getZero());
         return g;
     }
 
     /** Build a {@code Gradient} representing a variable.
      * <p>Instances built using this method are considered
      * to be the free variables with respect to which differentials
      * are computed. As such, their differential with respect to
      * themselves is +1.</p>
      * @param freeParameters number of free parameters (i.e. dimension of the gradient)
      * @param index index of the variable (from 0 to {@link #getFreeParameters() getFreeParameters()} - 1)
      * @param value value of the variable
      * @param <T> the type of the function parameters and value
      * @return a {@code FieldGradient} with a constant value and all derivatives set to 0.0 except the
      * one at {@code index} which will be set to 1.0
      */
     public static <T extends CalculusFieldElement<T>> FieldGradient<T> variable(final int freeParameters,
                                                                                 final int index, final T value) {
         final FieldGradient<T> g = new FieldGradient<>(value, freeParameters);
         final Field<T> field = value.getField();
         Arrays.fill(g.grad, field.getZero());
         g.grad[index] = field.getOne();
         return g;
     }
 
     /** Get the {@link Field} the value and parameters of the function belongs to.
      * @return {@link Field} the value and parameters of the function belongs to
      */
     public Field<T> getValueField() {
         return value.getField();
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> newInstance(final double c) {
         return newInstance(getValueField().getZero().newInstance(c));
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> newInstance(final T c) {
         return new FieldGradient<>(c, MathArrays.buildArray(value.getField(), grad.length));
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> withValue(final T v) {
         return new FieldGradient<>(v, grad);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> getAddendum() {
         return new FieldGradient<>(value.getField().getZero(), grad);
     }
 
     /** Get the value part of the function.
      * @return value part of the value of the function
      */
     @Override
     public T getValue() {
         return value;
     }
 
     /** Get the gradient part of the function.
      * @return gradient part of the value of the function
      */
     public T[] getGradient() {
         return grad.clone();
     }
 
     /** Get the number of free parameters.
      * @return number of free parameters
      */
     @Override
     public int getFreeParameters() {
         return grad.length;
     }
 
     /** {@inheritDoc} */
     @Override
     public T getPartialDerivative(final int ... orders)
         throws MathIllegalArgumentException {
 
         // check the number of components
         if (orders.length != grad.length) {
             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
                                                    orders.length, grad.length);
         }
 
         // check that either all derivation orders are set to 0,
         // or that only one is set to 1 and all other ones are set to 0
         int selected = -1;
         for (int i = 0; i < orders.length; ++i) {
             if (orders[i] != 0) {
                 if (selected >= 0 || orders[i] != 1) {
                      throw new MathIllegalArgumentException(LocalizedCoreFormats.DERIVATION_ORDER_NOT_ALLOWED,
                                                            orders[i]);
                 }
                 // found the component set to derivation order 1
                 selected = i;
             }
         }
 
         return (selected < 0) ? value : grad[selected];
 
     }
 
     /** Get the partial derivative with respect to one parameter.
      * @param n index of the parameter (counting from 0)
      * @return partial derivative with respect to the n<sup>th</sup> parameter
      * @exception MathIllegalArgumentException if n is either negative or larger
      * or equal to {@link #getFreeParameters()}
      */
     public T getPartialDerivative(final int n) throws MathIllegalArgumentException {
         if (n < 0 || n >= grad.length) {
             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, n, 0, grad.length - 1);
         }
         return grad[n];
     }
 
     /** Convert the instance to a {@link FieldDerivativeStructure}.
      * @return derivative structure with same value and derivative as the instance
      */
     public FieldDerivativeStructure<T> toDerivativeStructure() {
         final T[] derivatives = MathArrays.buildArray(getValueField(), 1 + grad.length);
         derivatives[0] = value;
         System.arraycopy(grad, 0, derivatives, 1, grad.length);
         return getField().getConversionFactory().build(derivatives);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> add(final double a) {
         return new FieldGradient<>(value.add(a), grad);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> add(final FieldGradient<T> a) {
         final FieldGradient<T> result = newInstance(value.add(a.value));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].add(a.grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> subtract(final double a) {
         return new FieldGradient<>(value.subtract(a), grad);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> subtract(final FieldGradient<T> a) {
         final FieldGradient<T> result = newInstance(value.subtract(a.value));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].subtract(a.grad[i]);
         }
         return result;
     }
 
     /** '&times;' operator.
      * @param n right hand side parameter of the operator
      * @return this&times;n
      */
     public FieldGradient<T> multiply(final T n) {
         final FieldGradient<T> result = newInstance(value.multiply(n));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].multiply(n);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> multiply(final int n) {
         final FieldGradient<T> result = newInstance(value.multiply(n));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].multiply(n);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> multiply(final double a) {
         final FieldGradient<T> result = newInstance(value.multiply(a));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].multiply(a);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> multiply(final FieldGradient<T> a) {
         final FieldGradient<T> result = newInstance(value.multiply(a.value));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].multiply(a.value).add(value.multiply(a.grad[i]));
         }
         return result;
     }
 
     /** '&divide;' operator.
      * @param a right hand side parameter of the operator
      * @return this&divide;a
      */
     public FieldGradient<T> divide(final T a) {
         final FieldGradient<T> result = newInstance(value.divide(a));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].divide(a);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> divide(final double a) {
         final FieldGradient<T> result = newInstance(value.divide(a));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].divide(a);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> divide(final FieldGradient<T> a) {
         final T inv1 = a.value.reciprocal();
         final T inv2 = inv1.multiply(inv1);
         final FieldGradient<T> result = newInstance(value.multiply(inv1));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].multiply(a.value).subtract(value.multiply(a.grad[i])).multiply(inv2);
         }
         return result;
     }
 
     /** IEEE remainder operator.
      * @param a right hand side parameter of the operator
      * @return this - n &times; a where n is the closest integer to this/a
      * (the even integer is chosen for n if this/a is halfway between two integers)
      */
     public FieldGradient<T> remainder(final T a) {
         return new FieldGradient<>(FastMath.IEEEremainder(value, a), grad);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> remainder(final double a) {
         return new FieldGradient<>(FastMath.IEEEremainder(value, a), grad);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> remainder(final FieldGradient<T> a) {
 
         // compute k such that lhs % rhs = lhs - k rhs
         final T rem = FastMath.IEEEremainder(value, a.value);
         final T k   = FastMath.rint(value.subtract(rem).divide(a.value));
 
         final FieldGradient<T> result = newInstance(rem);
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].subtract(k.multiply(a.grad[i]));
         }
         return result;
 
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> negate() {
         final FieldGradient<T> result = newInstance(value.negate());
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = grad[i].negate();
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> abs() {
         if (Double.doubleToLongBits(value.getReal()) < 0) {
             // we use the bits representation to also handle -0.0
             return negate();
         } else {
             return this;
         }
     }
 
     /**
      * Returns the instance with the sign of the argument.
      * A NaN {@code sign} argument is treated as positive.
      *
      * @param sign the sign for the returned value
      * @return the instance with the same sign as the {@code sign} argument
      */
     public FieldGradient<T> copySign(final T sign) {
         long m = Double.doubleToLongBits(value.getReal());
         long s = Double.doubleToLongBits(sign.getReal());
         if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
             return this;
         }
         return negate(); // flip sign
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> copySign(final FieldGradient<T> sign) {
         long m = Double.doubleToLongBits(value.getReal());
         long s = Double.doubleToLongBits(sign.value.getReal());
         if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
             return this;
         }
         return negate(); // flip sign
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> copySign(final double sign) {
         long m = Double.doubleToLongBits(value.getReal());
         long s = Double.doubleToLongBits(sign);
         if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
             return this;
         }
         return negate(); // flip sign
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> scalb(final int n) {
         final FieldGradient<T> result = newInstance(FastMath.scalb(value, n));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = FastMath.scalb(grad[i], n);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> hypot(final FieldGradient<T> y) {
 
         if (Double.isInfinite(value.getReal()) || Double.isInfinite(y.value.getReal())) {
             return newInstance(Double.POSITIVE_INFINITY);
         } else if (Double.isNaN(value.getReal()) || Double.isNaN(y.value.getReal())) {
             return newInstance(Double.NaN);
         } else {
 
             final int expX = getExponent();
             final int expY = y.getExponent();
             if (expX > expY + 27) {
                 // y is negligible with respect to x
                 return abs();
             } else if (expY > expX + 27) {
                 // x is negligible with respect to y
                 return y.abs();
             } else {
 
                 // find an intermediate scale to avoid both overflow and underflow
                 final int middleExp = (expX + expY) / 2;
 
                 // scale parameters without losing precision
                 final FieldGradient<T> scaledX = scalb(-middleExp);
                 final FieldGradient<T> scaledY = y.scalb(-middleExp);
 
                 // compute scaled hypotenuse
                 final FieldGradient<T> scaledH =
                         scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
 
                 // remove scaling
                 return scaledH.scalb(middleExp);
 
             }
 
         }
     }
 
     /** Compute composition of the instance by a function.
      * @param g0 value of the function at the current point (i.e. at {@code g(getValue())})
      * @param g1 first derivative of the function at the current point (i.e. at {@code g'(getValue())})
      * @return g(this)
      */
     @Override
     public FieldGradient<T> compose(final T g0, final T g1) {
         final FieldGradient<T> result = newInstance(g0);
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = g1.multiply(grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> rootN(final int n) {
         if (n == 2) {
             return sqrt();
         } else if (n == 3) {
             return cbrt();
         } else {
             final T r = FastMath.pow(value, 1.0 / n);
             return compose(r, FastMath.pow(r, n - 1).multiply(n).reciprocal());
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradientField<T> getField() {
         return FieldGradientField.getField(getValueField(), getFreeParameters());
     }
 
     /** Compute a<sup>x</sup> where a is a double and x a {@link FieldGradient}
      * @param a number to exponentiate
      * @param x power to apply
      * @param <T> the type of the function parameters and value
      * @return a<sup>x</sup>
      */
     public static <T extends CalculusFieldElement<T>> FieldGradient<T> pow(final double a, final FieldGradient<T> x) {
         if (a == 0) {
             return x.getField().getZero();
         } else {
             final T aX    = FastMath.pow(x.value.newInstance(a), x.value);
             final T aXlnA = aX.multiply(FastMath.log(a));
             final FieldGradient<T> result = x.newInstance(aX);
             for (int i = 0; i < x.grad.length; ++i) {
                 result.grad[i] =  aXlnA.multiply(x.grad[i]);
             }
             return result;
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> pow(final double p) {
         if (p == 0) {
             return getField().getOne();
         } else {
             final T f0Pm1 = FastMath.pow(value, p - 1);
             return compose(f0Pm1.multiply(value), f0Pm1.multiply(p));
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> pow(final int n) {
         if (n == 0) {
             return getField().getOne();
         } else {
             final T f0Nm1 = FastMath.pow(value, n - 1);
             return compose(f0Nm1.multiply(value), f0Nm1.multiply(n));
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldSinCos<FieldGradient<T>> sinCos() {
         final FieldSinCos<T> sinCos = FastMath.sinCos(value);
         final FieldGradient<T> sin = newInstance(sinCos.sin());
         final FieldGradient<T> cos = newInstance(sinCos.cos());
         final T mSin = sinCos.sin().negate();
         for (int i = 0; i < grad.length; ++i) {
             sin.grad[i] =  grad[i].multiply(sinCos.cos());
             cos.grad[i] =  grad[i].multiply(mSin);
         }
         return new FieldSinCos<>(sin, cos);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> atan2(final FieldGradient<T> x) {
         final T inv = value.square().add(x.value.multiply(x.value)).reciprocal();
         final FieldGradient<T> result = newInstance(FastMath.atan2(value, x.value));
         final T xValueInv = x.value.multiply(inv);
         final T mValueInv = value.negate().multiply(inv);
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = xValueInv.multiply(grad[i]).add(x.grad[i].multiply(mValueInv));
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldSinhCosh<FieldGradient<T>> sinhCosh() {
         final FieldSinhCosh<T> sinhCosh = FastMath.sinhCosh(value);
         final FieldGradient<T> sinh = newInstance(sinhCosh.sinh());
         final FieldGradient<T> cosh = newInstance(sinhCosh.cosh());
         for (int i = 0; i < grad.length; ++i) {
             sinh.grad[i] = grad[i].multiply(sinhCosh.cosh());
             cosh.grad[i] = grad[i].multiply(sinhCosh.sinh());
         }
         return new FieldSinhCosh<>(sinh, cosh);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> toDegrees() {
         final FieldGradient<T> result = newInstance(FastMath.toDegrees(value));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = FastMath.toDegrees(grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> toRadians() {
         final FieldGradient<T> result = newInstance(FastMath.toRadians(value));
         for (int i = 0; i < grad.length; ++i) {
             result.grad[i] = FastMath.toRadians(grad[i]);
         }
         return result;
     }
 
     /** Evaluate Taylor expansion of a gradient.
      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
      */
     public T taylor(final double... delta) {
         T result = value;
         for (int i = 0; i < grad.length; ++i) {
             result = result.add(grad[i].multiply(delta[i]));
         }
         return result;
     }
 
     /** Evaluate Taylor expansion of a gradient.
      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
      */
     public T taylor(@SuppressWarnings("unchecked") final T... delta) {
         T result = value;
         for (int i = 0; i < grad.length; ++i) {
             result = result.add(grad[i].multiply(delta[i]));
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final FieldGradient<T>[] a, final FieldGradient<T>[] b) {
 
         // extract values and first derivatives
         final Field<T> field = a[0].value.getField();
         final int n  = a.length;
         final T[] a0 = MathArrays.buildArray(field, n);
         final T[] b0 = MathArrays.buildArray(field, n);
         final T[] a1 = MathArrays.buildArray(field, 2 * n);
         final T[] b1 = MathArrays.buildArray(field, 2 * n);
         for (int i = 0; i < n; ++i) {
             final FieldGradient<T> ai = a[i];
             final FieldGradient<T> bi = b[i];
             a0[i]         = ai.value;
             b0[i]         = bi.value;
             a1[2 * i]     = ai.value;
             b1[2 * i + 1] = bi.value;
         }
 
         final FieldGradient<T> result = newInstance(a[0].value.linearCombination(a0, b0));
         for (int k = 0; k < grad.length; ++k) {
             for (int i = 0; i < n; ++i) {
                 a1[2 * i + 1] = a[i].grad[k];
                 b1[2 * i]     = b[i].grad[k];
             }
             result.grad[k] = a[0].value.linearCombination(a1, b1);
         }
         return result;
 
     }
 
     /**
      * Compute a linear combination.
      * @param a Factors.
      * @param b Factors.
      * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
      * @throws MathIllegalArgumentException if arrays dimensions don't match
      */
     public FieldGradient<T> linearCombination(final T[] a, final FieldGradient<T>[] b) {
 
         // extract values and first derivatives
         final Field<T> field = b[0].value.getField();
         final int      n  = b.length;
         final T[] b0 = MathArrays.buildArray(field, n);
         final T[] b1 = MathArrays.buildArray(field, n);
         for (int i = 0; i < n; ++i) {
             b0[i] = b[i].value;
         }
 
         final FieldGradient<T> result = newInstance(b[0].value.linearCombination(a, b0));
         for (int k = 0; k < grad.length; ++k) {
             for (int i = 0; i < n; ++i) {
                 b1[i] = b[i].grad[k];
             }
             result.grad[k] = b[0].value.linearCombination(a, b1);
         }
         return result;
 
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final double[] a, final FieldGradient<T>[] b) {
 
         // extract values and first derivatives
         final Field<T> field = b[0].value.getField();
         final int      n  = b.length;
         final T[] b0 = MathArrays.buildArray(field, n);
         final T[] b1 = MathArrays.buildArray(field, n);
         for (int i = 0; i < n; ++i) {
             b0[i] = b[i].value;
         }
 
         final FieldGradient<T> result = newInstance(b[0].value.linearCombination(a, b0));
         for (int k = 0; k < grad.length; ++k) {
             for (int i = 0; i < n; ++i) {
                 b1[i] = b[i].grad[k];
             }
             result.grad[k] = b[0].value.linearCombination(a, b1);
         }
         return result;
 
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final FieldGradient<T> a1, final FieldGradient<T> b1,
                                               final FieldGradient<T> a2, final FieldGradient<T> b2) {
         final FieldGradient<T> result = newInstance(a1.value.linearCombination(a1.value, b1.value,
                                                                                a2.value, b2.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             result.grad[i] = a1.value.linearCombination(a1.value,       b1.grad[i],
                                                             a1.grad[i], b1.value,
                                                             a2.value,       b2.grad[i],
                                                             a2.grad[i], b2.value);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final double a1, final FieldGradient<T> b1,
                                               final double a2, final FieldGradient<T> b2) {
         final FieldGradient<T> result = newInstance(b1.value.linearCombination(a1, b1.value,
                                                                                a2, b2.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             result.grad[i] = b1.value.linearCombination(a1, b1.grad[i],
                                                             a2, b2.grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final FieldGradient<T> a1, final FieldGradient<T> b1,
                                               final FieldGradient<T> a2, final FieldGradient<T> b2,
                                               final FieldGradient<T> a3, final FieldGradient<T> b3) {
         final Field<T> field = a1.value.getField();
         final T[] a = MathArrays.buildArray(field, 6);
         final T[] b = MathArrays.buildArray(field, 6);
         a[0] = a1.value;
         a[2] = a2.value;
         a[4] = a3.value;
         b[1] = b1.value;
         b[3] = b2.value;
         b[5] = b3.value;
         final FieldGradient<T> result = newInstance(a1.value.linearCombination(a1.value, b1.value,
                                                                                a2.value, b2.value,
                                                                                a3.value, b3.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             a[1] = a1.grad[i];
             a[3] = a2.grad[i];
             a[5] = a3.grad[i];
             b[0] = b1.grad[i];
             b[2] = b2.grad[i];
             b[4] = b3.grad[i];
             result.grad[i] = a1.value.linearCombination(a, b);
         }
         return result;
     }
 
     /**
      * Compute a linear combination.
      * @param a1 first factor of the first term
      * @param b1 second factor of the first term
      * @param a2 first factor of the second term
      * @param b2 second factor of the second term
      * @param a3 first factor of the third term
      * @param b3 second factor of the third term
      * @return a<sub>1</sub>&times;b<sub>1</sub> +
      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
      * @see #linearCombination(double, FieldGradient, double, FieldGradient)
      * @see #linearCombination(double, FieldGradient, double, FieldGradient, double, FieldGradient, double, FieldGradient)
      * @exception MathIllegalArgumentException if number of free parameters or orders are inconsistent
      */
     public FieldGradient<T> linearCombination(final T a1, final FieldGradient<T> b1,
                                               final T a2, final FieldGradient<T> b2,
                                               final T a3, final FieldGradient<T> b3) {
         final FieldGradient<T> result = newInstance(b1.value.linearCombination(a1, b1.value,
                                                                                a2, b2.value,
                                                                                a3, b3.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             result.grad[i] = b1.value.linearCombination(a1, b1.grad[i],
                                                         a2, b2.grad[i],
                                                         a3, b3.grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final double a1, final FieldGradient<T> b1,
                                               final double a2, final FieldGradient<T> b2,
                                               final double a3, final FieldGradient<T> b3) {
         final FieldGradient<T> result = newInstance(b1.value.linearCombination(a1, b1.value,
                                                                                a2, b2.value,
                                                                                a3, b3.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             result.grad[i] = b1.value.linearCombination(a1, b1.grad[i],
                                                             a2, b2.grad[i],
                                                             a3, b3.grad[i]);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final FieldGradient<T> a1, final FieldGradient<T> b1,
                                               final FieldGradient<T> a2, final FieldGradient<T> b2,
                                               final FieldGradient<T> a3, final FieldGradient<T> b3,
                                               final FieldGradient<T> a4, final FieldGradient<T> b4) {
         final Field<T> field = a1.value.getField();
         final T[] a = MathArrays.buildArray(field, 8);
         final T[] b = MathArrays.buildArray(field, 8);
         a[0] = a1.value;
         a[2] = a2.value;
         a[4] = a3.value;
         a[6] = a4.value;
         b[1] = b1.value;
         b[3] = b2.value;
         b[5] = b3.value;
         b[7] = b4.value;
         final FieldGradient<T> result = newInstance(a1.value.linearCombination(a1.value, b1.value,
                                                                                a2.value, b2.value,
                                                                                a3.value, b3.value,
                                                                                a4.value, b4.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             a[1] = a1.grad[i];
             a[3] = a2.grad[i];
             a[5] = a3.grad[i];
             a[7] = a4.grad[i];
             b[0] = b1.grad[i];
             b[2] = b2.grad[i];
             b[4] = b3.grad[i];
             b[6] = b4.grad[i];
             result.grad[i] = a1.value.linearCombination(a, b);
         }
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> linearCombination(final double a1, final FieldGradient<T> b1,
                                               final double a2, final FieldGradient<T> b2,
                                               final double a3, final FieldGradient<T> b3,
                                               final double a4, final FieldGradient<T> b4) {
         final FieldGradient<T> result = newInstance(b1.value.linearCombination(a1, b1.value,
                                                                                a2, b2.value,
                                                                                a3, b3.value,
                                                                                a4, b4.value));
         for (int i = 0; i < b1.grad.length; ++i) {
             result.grad[i] = b1.value.linearCombination(a1, b1.grad[i],
                                                             a2, b2.grad[i],
                                                             a3, b3.grad[i],
                                                             a4, b4.grad[i]);
         }
         return result;
     }
 
     /**
      * Add an independent variable to the Taylor expansion.
      * @return object with one more variable
      * @since 4.0
      */
     public FieldGradient<T> stackVariable() {
         final T[] gradient = MathArrays.buildArray(getValueField(), this.getFreeParameters() + 1);
         System.arraycopy(this.grad, 0, gradient, 0, this.getFreeParameters());
         return new FieldGradient<>(this.value, gradient);
     }
 
     /** {@inheritDoc} */
     @Override
     public FieldGradient<T> getPi() {
         return new FieldGradient<>(getValueField().getZero().getPi(), getFreeParameters());
     }
 
     /** Test for the equality of two univariate derivatives.
      * <p>
      * univariate derivatives are considered equal if they have the same derivatives.
      * </p>
      * @param other Object to test for equality to this
      * @return true if two univariate derivatives are equal
      */
     @Override
     public boolean equals(Object other) {
 
         if (this == other) {
             return true;
         }
 
         if (other instanceof FieldGradient) {
             @SuppressWarnings("unchecked")
             final FieldGradient<T> rhs = (FieldGradient<T>) other;
             if (!value.equals(rhs.value) || grad.length != rhs.grad.length) {
                 return false;
             }
             for (int i = 0; i < grad.length; ++i) {
                 if (!grad[i].equals(rhs.grad[i])) {
                     return false;
                 }
             }
             return true;
         }
 
         return false;
 
     }
 
     /** Get a hashCode for the univariate derivative.
      * @return a hash code value for this object
      */
     @Override
     public int hashCode() {
         return 129 + 7 *value.hashCode() - 15 * Arrays.hashCode(grad);
     }
 
 }