#+LATEX_CLASS_OPTIONS: [a4paper,12pt]
#+LATEX_HEADER: \usepackage[margin=0.5in]{geometry}
#+LATEX_HEADER: \usepackage{fontspec}
#+LATEX_HAEDER: \usepackage{graphicx}
#+LATEX_HEADER: \setmainfont{LiberationSerif}
#+LATEX_HEADER: \date{}
#+OPTIONS: toc:nil

#+INCLUDE: front.tex
#+LATEX: \clearpage

#+INCLUDE: toc.tex
#+LATEX: \clearpage

|------------------------------------------------------------------------------------------------------------|
| Programs are followed by their respective inputs and outputs i.e, both stdin and stdout are shown together |
|------------------------------------------------------------------------------------------------------------|

* Program for finding roots of f(x) = 0 using Newton Raphson Method
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C :tangle 1.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>
#define EPSILON 0.0000001
#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)
#define f1(x) ((1056 * x * x) - (128 * x) + 198)
double newtonRaphson(double x) {
  double h = f(x) / f1(x);
  if (f(x) == 0 || fabs(h) < EPSILON)
    return x;
  return newtonRaphson(x - h);
}
int main() {
  printf("Root for f(x) = 352x^3 - 64x^2 + 198x - 36 is %lf",
         newtonRaphson(-4));
}
#+end_src

* Program for nding roots of f(x) = 0 using bisection method
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C :tangle 2.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define EPSILON 0.0000001
#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)
double bisection(double a, double b) {
  double x;
  if (f(a) * f(b) > 0) {
    printf("The values of function at the respective initial guesses must have "
           "opposite signs");
    exit(1);
  }
  x = (a + b) / 2;
  if (f(x) == 0 || fabs(b - a) < EPSILON)
    return x;
  if (f(x) > 0)
    return bisection(x, b);
  return bisection(a, x);
}
int main() {
  printf("Root for f(x) = 352x^3 - 64x^2 + 198x - 36 is %lf",
         bisection(1.6, -4));
}
#+end_src

#+LATEX: \clearpage

* Program for finding roots of f(x) = 0 using secant method
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C :tangle 3.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>
#define EPSILON 0.0000001
#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)
double secant(double a, double b) {
  double x1;
  x1 = (a * f(b) - b * f(a)) / (f(b) - f(a));
  if (f(x1) == 0 || fabs(a - b) < EPSILON)
    return x1;
  return secant(b, x1);
}
int main() {
  printf("Root for f(x) = 352x^3 - 64x^2 + 198x - 36 is %lf", secant(1.6, -4));
}
#+end_src

#+LATEX: \clearpage

* Program to implement Langrange Interpolation.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C :tangle 4.c :results output :wrap src text
#include <stdio.h>

int main() {
  double xp, yp = 0, p;
  int i, j, n;

  printf("Number of inputs: ");
  scanf("%d", &n);

  double x[n], y[n];
  printf("Input sample space:\n");

  for (i = 0; i < n; i++) {
    printf("x%d: ", i);
    scanf("%lf", x + i);
    printf("y%d: ", i);
    scanf("%lf", y + i);
  }

  printf("Enter interpolation point x: ");
  scanf("%lf", &xp);

  for (i = 0; i < n; i++) {
    p = 1;
    for (j = 0; j < n; j++) {
      if (i != j) {
        p *= (xp - x[j]) / (x[i] - x[j]);
      }
    }
    yp += p * y[i];
  }
  printf("Interpolated value for %lf is %lf.", xp, yp);
}
#+end_src

#+begin_src text
Number of inputs: 4
Input sample space:
x0: 0
y0: 2
x1: 1
y1: 3
x2: 2
y2: 12
x3: 5
y3: 147
Enter interpolation point x: 3
Interpolated value for 3.000000 is 35.000000.
#+end_src

#+LATEX: \clearpage

* Program to implement Newton's Divided Difference formula.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C :tangle 5.c :results output :wrap src text
#include <stdio.h>

int main() {
  int n, i, j = 1;
  double xp, yp, f1, f2 = 0;

  printf("Enter the number of inputs: ");
  scanf("%d", &n);

  double x[n], y[n];

  printf("Enter input values:\n");
  for (i = 0; i < n; i++) {
    printf("x%d=", i);
    scanf("%lf", x + i);
    printf("y%d=", i);
    scanf("%lf", y + i);
  }

  yp = y[0];

  printf("Enter interpolation point x: ");
  scanf("%lf", &xp);

  do {
    for (i = 0; i < n - 1; i++) 
      y[i] = ((y[i + 1] - y[i]) / (x[i + j] - x[i]));

    f1 = 1;

    for (i = 0; i < j; i++) {
      f1 *= (xp - x[i]);
    }

    f2 += (y[0] * f1);
    j++;
  } while ((--n) > 1);

  yp += f2;
  
  printf("Interpolated value for %lf is %lf.", xp, yp);
}
#+end_src

#+LATEX: \clearpage

#+begin_src text
Enter the number of inputs: 4
Enter input values:
x0=3
y0=9
x1=5
y1=12 
x2=9
y2=666
x3=15
y3=10245
Enter interpolation point x: 999
Interpolated value for 999.000000 is 9525764925.000002.
#+end_src

* Program for solving numerical integration by trapezoidal method.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C  :tangle 6.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>

#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)

double trapezoidal_integral(double a, double b, double n) {
  double h = (b - a) / n;

  double s = (f(a) + f(b)) / 2;

  // Add the other heights

  for (int i = 1; i < n; i++)
    s += f(a + i * h);

  return s * h;
}

int main() {
  printf("The area under the curve f(x) = 352x^3 - 64x^2 + 198x - 36 from x=3 "
         "to x=4 is %lf",
         trapezoidal_integral(3, 4, 10000));
}
#+end_src

* Program for solving numerical integration by Simpson's 1/3 rule.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C  :tangle 7.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>

#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)

double simpsons_1_3(double a, double b, double n) {
  double h = (b - a) / n;

  double s = f(a) + f(b);

  // Add the other heights

  for (int i = 1; i < n; i += 2)
    s += 4 * f(a + i * h);

  for (int i = 2; i < n; i += 2)
    s += 2 * f(a + i * h);

  return (h / 3) * s;
}

int main() {
  printf("The area under the curve f(x) = 352x^3 - 64x^2 + 198x - 36 from x=3 "
         "to x=4 is %lf",
         simpsons_1_3(3, 4, 10000));
}
#+end_src

#+LATEX: \clearpage
* Program for solving numerical integration by Simpson's 3/8 rule.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C  :tangle 8.c :results output :exports both :wrap src text
#include <math.h>
#include <stdio.h>

#define f(x) ((352 * x * x * x) - (64 * x * x) + (198 * x) - 36)

double simpsons_1_3(double a, double b, double n) {
  double h = (b - a) / n;

  double s = f(a) + f(b);

  // Add the other heights

  for (int i = 1; i < n; i++)
    s += 3 * f(a + i * h);

  for (int i = 3; i < n; i += 3)
    s -= f(a + i * h);

  return (h / 8) * 3 * s;
}

int main() {
  printf("The area under the curve f(x) = 352x^3 - 64x^2 + 198x - 36 from x=3 "
         "to x=4 is %lf",
         simpsons_1_3(3, 4, 10000));
}
#+end_src

* Program for finding inverse of linear equations using Gauss Jordan method.
#+ATTR_LATEX: :options frame=single,breaklines=true
#+begin_src C  :tangle 9.c :results output :exports both :wrap src text
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

double **inverse(double **matrix, int order) {
  double **inverse = calloc(order, sizeof(double *));
  for (int i = 0; i < order; i++) {
    inverse[i] = calloc(2 * order, sizeof(double));
    inverse[i][order + i] = 1;
    memcpy(inverse[i], matrix[i], order * sizeof(double));
  }

  for (int i = 0; i < order; i++) {
    for (int j = 0; j < order; j++) {
      if (i == j)
        continue;
      double r = inverse[j][i] / inverse[i][i];

      for (int k = 0; k < order * 2; k++)
        inverse[j][k] -= r * inverse[i][k];
    }
  }

  for (int i = 0; i < order; i++) {
    for (int j = 0; j < order; j++)
      inverse[i][j + order] /= inverse[i][i];
  }

  return inverse;
}

int main() {
  const int ORDER = 3;

  double **matrix = malloc(ORDER * sizeof(double *));
  for (int i = 0; i < ORDER; i++)
    matrix[i] = malloc(ORDER * sizeof(double));

  matrix[0][0] = 92, matrix[0][1] = 4.5, matrix[0][2] = 61; // 92x + 4.5y + 61z = 0
  matrix[1][0] = -2, matrix[1][1] = 0, matrix[1][2] = 92387; // -2x + 92387z = 0
  matrix[2][0] = -2, matrix[2][1] = 0, matrix[2][2] = -23; // -2x - 23z = 0

  double **inv = inverse(matrix, ORDER);

  for (int i = 0; i < ORDER; i++) {
    for (int j = 0; j < ORDER; j++)
      printf("%lf ", inv[i][j + ORDER]);
    printf("\n");
    free(inv[i]);
    free(matrix[i]);
  }
  free(inv);
  free(matrix);
}
#+end_src