0_Fundamentals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ★ Fundamentals ★"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0.1 Horner’s method or  Nested multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import traceback\n",
    "import math\n",
    "import numpy as np\n",
    "import unittest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "testcase1 (__main__.nest_unittest) ... ok\n",
      "testcase2 (__main__.nest_unittest) ... ok\n",
      "testcase3 (__main__.nest_unittest) ... ok\n",
      "\n",
      "----------------------------------------------------------------------\n",
      "Ran 3 tests in 0.009s\n",
      "\n",
      "OK\n"
     ]
    }
   ],
   "source": [
    "def nest(degree, coefficients, x = 0, base_points = None) -> float:\n",
    "    \"\"\"\n",
    "    Evaluates polynomial from nested form using Horner’s Method\n",
    "    \n",
    "    Examples:\n",
    "        P(x) = 3 * x^2 + 5 * x − 1 and evaluate P(x = 1)\n",
    "        Use nest(2, [3, 5, -1], 1) to get the value of above polynomial\n",
    "    \n",
    "    Arguments:\n",
    "        degree (int): degree of polynomial\n",
    "        coefficients (list): list of coefficients \n",
    "        x (float): x-coordinate x at which to evaluate (Default : 0)\n",
    "        base_points (list): list of base points (Default : None)\n",
    "        \n",
    "    Return:\n",
    "         value y of polynomial at x\n",
    "         \n",
    "    Raises:\n",
    "        ValueError:\n",
    "            coefficients is null\n",
    "            degree is negative\n",
    "            degree is not equal len(base_points)\n",
    "            degree len(coefficients) is not equal len(base_points) + 1\n",
    "    \"\"\"\n",
    "    try:\n",
    "        if base_points is None :\n",
    "            base_points = np.zeros(degree).tolist()\n",
    "        \n",
    "        if(degree < 0):\n",
    "            raise ValueError('degree is negative')\n",
    "        \n",
    "        if(degree != len(base_points)):\n",
    "            raise ValueError('degree is not consistent with base points')\n",
    "        \n",
    "        if coefficients is None :\n",
    "            raise ValueError('coefficients is null')\n",
    "        \n",
    "        if(degree + 1 != len(coefficients)):\n",
    "            raise ValueError\n",
    "            \n",
    "        if(len(coefficients) != len(base_points) + 1):\n",
    "            raise ValueError \n",
    "        \n",
    "        # Check whether coefficients is type of ndarray\n",
    "        if(type(coefficients).__module__ == np.__name__):\n",
    "            coefficients = coefficients.tolist()\n",
    "        \n",
    "        # Check whether base_points is type of ndarray\n",
    "        if(type(base_points).__module__ == np.__name__):\n",
    "            base_points = base_points.tolist()\n",
    "            \n",
    "        y = coefficients.pop(0)\n",
    "        \n",
    "        for i in range(degree):\n",
    "            y = y * (x - base_points[i]) + coefficients[i]\n",
    "            \n",
    "        return y\n",
    "    except ValueError as e:\n",
    "        print('Exception : ValueError {0}'.format(str(e)))\n",
    "        traceback.print_exception()\n",
    "\n",
    "# unittest nest(...) \n",
    "class nest_unittest(unittest.TestCase):\n",
    "    # P(x = 0) = 3\n",
    "    def testcase1(self):\n",
    "        self.assertEqual(nest(0, [3], 0), 3)\n",
    "    # P(x = 2.1) = x\n",
    "    def testcase2(self):\n",
    "        self.assertEqual(nest(1, [1, 0], 2.1), 2.1)\n",
    "    # P(x = 4) = 3 * x^2 + 5 * x + 7\n",
    "    def testcase3(self):\n",
    "        self.assertEqual(nest(2, [3, 5, 7], 4), 75)\n",
    "\n",
    "unittest.main(argv=['ignored', '--verbose'], exit=False);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example : $P(x) = 2x^{4} + 3x^{3} - 3x^{2} + 5x - 1$\n",
    "Evaluate $P(0.5)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.25"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nest(4, [2, 3, -3, 5, -1], 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example \n",
    "\n",
    "$P(x) = x^6 - 2x^5 + 3x^4 - 4x^3 + 5x^2 - 6x + 7$ at $x = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nest(6, [1, -2, 3, -4, 5, -6, 7], 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.1 Computer Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate $P(x) = 1 + x + ... + x^{50}$ at x = 1.00001. <br/>\n",
    "Find the error of the computation by comparing with the equivalent expression $Q(x) = (x^{51} - 1) / (x - 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(1.00001) = 51.01275208274999\n",
      "Q(1.00001) = 51.01275208274523\n",
      "error      = 00.00000000000476\n"
     ]
    }
   ],
   "source": [
    "x = 1.00001\n",
    "p = nest(50, np.ones(51, dtype=np.int), x) \n",
    "q = (x ** 51 - 1) / (x - 1)\n",
    "print('P(1.00001) = {0}'.format(p))\n",
    "print('Q(1.00001) = {0}'.format(q))\n",
    "print('error      = {0:017.14f}'.format(np.abs(p - q)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Evaluate $P(x) = 1 - x + x^2 - x^3 + ... + x^{98} - x^{99}$ at x = 1.00001. <br/>\n",
    "Find a simpler, equivalent expression, and use it to estimate the error of the nested multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(1.00001) = -0.0005002450796476321\n",
      "Q(1.00001) = -0.0005002450796474608\n",
      "error      =  0.0000000000000001713\n"
     ]
    }
   ],
   "source": [
    "x = 1.00001\n",
    "p = nest(99, [-1 if i % 2 == 0 else 1 for i in range(100)], x)\n",
    "q = (1 - x ** 100) / (1 + x) # q(x) = (1 - x^100) / (1 + x) = p(x)\n",
    "print('P(1.00001) = {0}'.format(p))\n",
    "print('Q(1.00001) = {0}'.format(q))\n",
    "print('error      = {0:22.19f}'.format(np.abs(p - q)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0.4 Loss of significance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.4 Computer Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Calculate the expressions that follow in double precision arithmetic $x = 10^{-1},\\cdots,10^{-14}$. Then, using an alternative form of the expression that doesn’t suffer from subtracting nearly equal numbers, repeat the calculation and make a table of results. Report the number of correct digits in the original expression for each $x$. <br/>\n",
    "## $(a)\\frac{1 - secx}{tan^{2}x}$\n",
    "## $(b)\\frac{1 - (1 - x)^{3}}{x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x = 0.10000000000000) = -0.49874791371141\n",
      "P(x = 0.01000000000000) = -0.49998749979096\n",
      "P(x = 0.00100000000000) = -0.49999987501429\n",
      "P(x = 0.00010000000000) = -0.49999999362793\n",
      "P(x = 0.00001000000000) = -0.50000004133685\n",
      "P(x = 0.00000100000000) = -0.50004445029084\n",
      "P(x = 0.00000010000000) = -0.51070259132757\n",
      "P(x = 0.00000001000000) = 0.00000000000000\n",
      "P(x = 0.00000000100000) = 0.00000000000000\n",
      "P(x = 0.00000000010000) = 0.00000000000000\n",
      "P(x = 0.00000000001000) = 0.00000000000000\n",
      "P(x = 0.00000000000100) = 0.00000000000000\n",
      "P(x = 0.00000000000010) = 0.00000000000000\n",
      "P(x = 0.00000000000001) = 0.00000000000000\n"
     ]
    }
   ],
   "source": [
    "Px_a = lambda x : (1 - np.cos(x) ** -1) / np.tan(x) ** 2\n",
    "\n",
    "for i in range(1, 15):\n",
    "    x = 10 ** (-i)\n",
    "    print( 'P(x = {0:>.14f}) = {1:>.14f}'.format(x, Px_a(x)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative form :  \n",
    "## $\\frac{1 - secx}{tan^{2}x} \\Rightarrow \\frac{cosx \\cdot (1 - \\frac{1}{cosx})}{cosx \\cdot (\\frac{sinx}{cosx})^2} \\Rightarrow \\frac{cosx \\cdot (1 - \\frac{1}{cosx})}{\\frac{sin^2{x}}{cosx}} \\Rightarrow \\frac{cos^{2}x \\cdot (1 - \\frac{1}{cosx})}{sin^2{x}} \\Rightarrow \\frac{-\\cos{x} \\cdot (1 - cosx)}{1 - cos^{2}x} \\Rightarrow \\frac{-cosx}{(1 + cosx)} \\Rightarrow \\frac{-1}{(1 + secx)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x = 0.10000000000000) = -0.49874791371143\n",
      "P(x = 0.01000000000000) = -0.49998749979166\n",
      "P(x = 0.00100000000000) = -0.49999987499998\n",
      "P(x = 0.00010000000000) = -0.49999999875000\n",
      "P(x = 0.00001000000000) = -0.49999999998750\n",
      "P(x = 0.00000100000000) = -0.49999999999987\n",
      "P(x = 0.00000010000000) = -0.50000000000000\n",
      "P(x = 0.00000001000000) = -0.50000000000000\n",
      "P(x = 0.00000000100000) = -0.50000000000000\n",
      "P(x = 0.00000000010000) = -0.50000000000000\n",
      "P(x = 0.00000000001000) = -0.50000000000000\n",
      "P(x = 0.00000000000100) = -0.50000000000000\n",
      "P(x = 0.00000000000010) = -0.50000000000000\n",
      "P(x = 0.00000000000001) = -0.50000000000000\n"
     ]
    }
   ],
   "source": [
    "Px_a = lambda x : -1 / (1 + np.cos(x) ** -1)\n",
    "\n",
    "for i in range(1, 15):\n",
    "    x = 10 ** (-i)\n",
    "    print( 'P(x = {0:>.14f}) = {1:>.14f}'.format(x, Px_a(x)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x = 0.10000000000000) = 2.71000000000000\n",
      "P(x = 0.01000000000000) = 2.97010000000000\n",
      "P(x = 0.00100000000000) = 2.99700100000000\n",
      "P(x = 0.00010000000000) = 2.99970001000016\n",
      "P(x = 0.00001000000000) = 2.99997000008378\n",
      "P(x = 0.00000100000000) = 2.99999700004161\n",
      "P(x = 0.00000010000000) = 2.99999969866072\n",
      "P(x = 0.00000001000000) = 2.99999998176759\n",
      "P(x = 0.00000000100000) = 2.99999991515421\n",
      "P(x = 0.00000000010000) = 3.00000024822111\n",
      "P(x = 0.00000000001000) = 3.00000024822111\n",
      "P(x = 0.00000000000100) = 2.99993363483964\n",
      "P(x = 0.00000000000010) = 3.00093283556180\n",
      "P(x = 0.00000000000001) = 2.99760216648792\n"
     ]
    }
   ],
   "source": [
    "Px_b = lambda x : (1 - (1 - x) ** 3) / x\n",
    "\n",
    "for i in range(1, 15):\n",
    "    x = 10 ** (-i)\n",
    "    print( 'P(x = {0:>.14f}) = {1:>.14f}'.format(x, Px_b(x)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative form :  \n",
    "## $\\frac{1 - (1 - x)^{3} }{x} \\Rightarrow \\frac{1 - (1 - 3x + 3x^2 - x^3)}{x} \\Rightarrow x^2 - 3x + 3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x = 0.10000000000000) = 2.71000000000000\n",
      "P(x = 0.01000000000000) = 2.97010000000000\n",
      "P(x = 0.00100000000000) = 2.99700100000000\n",
      "P(x = 0.00010000000000) = 2.99970001000000\n",
      "P(x = 0.00001000000000) = 2.99997000010000\n",
      "P(x = 0.00000100000000) = 2.99999700000100\n",
      "P(x = 0.00000010000000) = 2.99999970000001\n",
      "P(x = 0.00000001000000) = 2.99999997000000\n",
      "P(x = 0.00000000100000) = 2.99999999700000\n",
      "P(x = 0.00000000010000) = 2.99999999970000\n",
      "P(x = 0.00000000001000) = 2.99999999997000\n",
      "P(x = 0.00000000000100) = 2.99999999999700\n",
      "P(x = 0.00000000000010) = 2.99999999999970\n",
      "P(x = 0.00000000000001) = 2.99999999999997\n"
     ]
    }
   ],
   "source": [
    "Px_b = lambda x : x ** 2 - 3 * x + 3\n",
    "\n",
    "for i in range(1, 15):\n",
    "    x = 10 ** (-i)\n",
    "    print( 'P(x = {0:>.14f}) = {1:>.14f}'.format(x, Px_b(x)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "2. Find the smallest value of $p$ for which the expression calculated in double precision arithmetic at $x = 10^{-p}$ has no correct significant digits. <br/>\n",
    "\n",
    "## $(a)\\frac{tanx - x}{x^3}$\n",
    "## $(b)\\frac{e^{x} + cos{x} - sin{x} - 2}{x^3}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x = 1.000000e-01) = 0.33467208545054\n",
      "P(x = 1.000000e-02) = 0.33334666720702\n",
      "P(x = 1.000000e-03) = 0.33333346673159\n",
      "P(x = 1.000000e-04) = 0.33333333651891\n",
      "P(x = 1.000000e-05) = 0.33333287573298\n",
      "P(x = 1.000000e-06) = 0.33330746474457\n",
      "P(x = 1.000000e-07) = 0.33087224502121\n",
      "P(x = 1.000000e-08) = 0.00000000000000\n",
      "P(x = 1.000000e-09) = 0.00000000000000\n",
      "P(x = 1.000000e-10) = 0.00000000000000\n",
      "The smallest value of p is 8\n",
      "\n",
      "P(x = 1.000000e-01) = 0.34166670684543\n",
      "P(x = 1.000000e-02) = 0.33416666678221\n",
      "P(x = 1.000000e-03) = 0.33341684968491\n",
      "P(x = 1.000000e-04) = 0.33351099659740\n",
      "P(x = 1.000000e-05) = 0.44408920985006\n",
      "P(x = 1.000000e-06) = 0.00000000000000\n",
      "P(x = 1.000000e-07) = -222044.60492503116257\n",
      "P(x = 1.000000e-08) = 0.00000000000000\n",
      "P(x = 1.000000e-09) = 0.00000000000000\n",
      "P(x = 1.000000e-10) = 0.00000000000000\n",
      "The smallest value of p is 6\n"
     ]
    }
   ],
   "source": [
    "def Px2(x, Px):\n",
    "    while True:\n",
    "        yield x, Px(x)\n",
    "        x /= 10\n",
    "\n",
    "Px2a_generator = Px2(0.1, lambda x : (np.tan(x) - x) / np.power(x, 3))\n",
    "\n",
    "for i in range(10):\n",
    "    print( 'P(x = {0:>e}) = {1:>.14f}'.format(*next(Px2a_generator)))\n",
    "print( 'The smallest value of p is 8' )\n",
    "print()\n",
    "\n",
    "Px2b_generator = Px2(0.1, lambda x : (np.exp(x) + np.cos(x) - np.sin(x) - 2) / np.power(x, 3))\n",
    "\n",
    "for i in range(10):\n",
    "    print( 'P(x = {0:>e}) = {1:>.14f}'.format(*next(Px2b_generator)))\n",
    "print( 'The smallest value of p is 6' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "3. Evaluate the quantity $ a + \\sqrt{a^2 + b^2} $ to four correct signiﬁcant digits, where $ a = -12345678987654321 $ and $ b = 123 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ a + \\sqrt{a^2 + b^2} \\Rightarrow \\frac{(a - \\sqrt{a^2 + b^2})(a + \\sqrt{a^2 + b^2})}{(a - \\sqrt{a^2 + b^2})} \\Rightarrow \\frac{a^2 - (a^2 + b^2)}{(a - \\sqrt{a^2 + b^2})} \\Rightarrow \\frac{-b^2}{(a - \\sqrt{a^2 + b^2})}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.127245e-13\n"
     ]
    }
   ],
   "source": [
    "Px3 = lambda a, b : - np.power(b, 2) / (a - np.sqrt(a * a + b * b))\n",
    "print('{0:e}'.format(Px3(-12345678987654321.0, 123.0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "4. Evaluate the quantity $ \\sqrt{c^2 + d} - c $ to four correct signiﬁcant digits, where $c = 246886422468$ and $d = 13579$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ \\sqrt{c^2 + d} - c \\Rightarrow \\frac{(\\sqrt{c^2 + d} - c)(\\sqrt{c^2 + d} + c)}{(\\sqrt{c^2 + d} + c)} \\Rightarrow \\frac{(c^2 + d) - c^2}{\\sqrt{c^2 + d} + c} \\Rightarrow \\frac{d}{\\sqrt{c^2 + d} + c}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.750050e-08\n"
     ]
    }
   ],
   "source": [
    "Px4 = lambda c, d : d / (np.sqrt(c * c + d) + c)\n",
    "print('{0:e}'.format(Px4(246886422468.0, 13579.0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "5. Consider a right triangle whose legs are of length 3344556600 and 1.2222222. How much longer is the hypotenuse than the longer leg ? Give your answer with at least four correct digits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$a = 3344556600$, $b = 1.2222222$, the hypotenuse $c$ is that $c^2 = a^2 + b^2$ <br/>\n",
    "\n",
    "And $c - a \\Rightarrow \\sqrt{a^2 + b^2} - a \\Rightarrow \\frac{(\\sqrt{a^2 + b^2} - a)(\\sqrt{a^2 + b^2} + a)}{(\\sqrt{a^2 + b^2} + a)} \\Rightarrow \\frac{b^2}{\\sqrt{a^2 + b^2} + a}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.233221e-10\n"
     ]
    }
   ],
   "source": [
    "diff = lambda a, b : np.power(b, 2) / (np.sqrt(a * a + b * b) + a)\n",
    "print('{0:e}'.format(diff(3344556600.0, 1.2222222)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## MIT License\n",
    "\n",
    "Copyright (c) Jim00000\n",
    "\n",
    "Permission is hereby granted, free of charge, to any person obtaining a copy\n",
    "of this software and associated documentation files (the \"Software\"), to deal\n",
    "in the Software without restriction, including without limitation the rights\n",
    "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
    "copies of the Software, and to permit persons to whom the Software is\n",
    "furnished to do so, subject to the following conditions:\n",
    "\n",
    "The above copyright notice and this permission notice shall be included in all\n",
    "copies or substantial portions of the Software."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}