#include <iostream>
#include <cstdlib>
#include "lilib.h"

using namespace std;

// Jn ^ (m) (x)
LongInterval besselj(int n, LongInterval x, int m){
  int sign, absN, s, k;
  LongFloat epsilon;
  LongInterval x2, x22, j, deltaJ, deltaJOld;

  if(m < 0){
    cerr << "[ERROR] besselj(int n, LongInterval x, int m) : m must be not negative." << endl;
    exit(1);
  }

  epsilon = 10;
  epsilon = pow(epsilon, -lilib::getPrecision());

  if(n % 2 == -1){
    sign = -1;
  }
  else{
    sign = 1;
  }

  if(n > 0){
    absN = n;
  }
  else{
    absN = -n;
  }

  x2 = x / 2;
  x22 = x2 * x2;

  if(m < absN){
    s = 0;
  }
  else{
    s = (m - absN + 1) / 2;
  }

  deltaJ = pow(x2, absN + 2 * s - m);
  for(k = 0; k < absN; k++){
    deltaJ /= k + 1;
  }
  for(k = 0; k < s; k++){
    deltaJ /= -(k + 1) * (absN + k + 1);
  }
  for(k = 0; k < m; k++){
    deltaJ *= (absN + 2 * s - k);
    deltaJ /= 2;
  }

  j = deltaJ;
  k = s + 1;

  deltaJOld = 0;
  while(deltaJ.mag() > deltaJOld.mig()){
    deltaJOld = deltaJ;

    deltaJ *= -(absN + 2 * k) * (absN + 2 * k - 1) * x22 / (k * (absN + k) * (absN + 2 * k - m) * (absN + 2 * k - m - 1));

    j += deltaJ;
    k++;
  }

  while(deltaJ.mag() > (deltaJOld * epsilon).mig()){
    deltaJ *= -(absN + 2 * k) * (absN + 2 * k - 1) * x22 / (k * (absN + k) * (absN + 2 * k - m) * (absN + 2 * k - m - 1));

    j += deltaJ;
    k++;
  }

  deltaJ *= -(absN + 2 * k) * (absN + 2 * k - 1) * x22 / (k * (absN + k) * (absN + 2 * k - m) * (absN + 2 * k - m - 1));

  j = sign * (j + deltaJ.mag());
//  j = longInfSup(max(inf(j), -1), min(sup(j), 1));

  return j;
}

int main(){
  lilib::setPrecision(100);
  cout << "Long precision is " << lilib::getPrecision() << "." << endl;
  cout << endl;

  cout << besselj(-3, 30, 5) << endl;

  return 0;
}