HP Prime: Minimum Distance Between a Point and a Line

HP Prime: Minimum Distance Between a Point and a Line

Introduction

We have a line in the form of y = m * x + b, where m is the slope of the line and b is the y-intercept of the line, and a separate point (px, py). The task is to find the minimum distance, or the shortest distance, between the point and the line. The separate point is not required to be on the line. The line and point are in two-dimensional space.

If the point (px, py) is not on the line, then theoretically, there are an infinite amount of distances between the point and the line. However, to get the shortest distance, draw a path that is “directly straight” to the line. This is achieved by choosing a line that connects the (px, py) that is a line that is orthogonal (perpendicular) to the line y = m * x + b.

The line drawn is of the form y = -1/m * x + b1. The slope of the orthogonal line is -1/m. Assuming that m ≠ 0, the y-intercept of the orthogonal line is b1 = y1 + x1 / m.

The next step is to find where the two lines intersect, which is done by solving the following system for x and y:

y = m * x + b

y = -m / n + b1

Label the intersection point (x1, y1). The minimum distance will be calculated as follows:

dist = √( (x1 – px)^2 + (y1 – py)^2 ) = abs( (x1 – px) + (y1 – py)*i)

If m = 0, the line is in the form of y = b. The orthogonal line is x = px, and the distance is simply abs( (y1 – py)*i ).

HP Prime Code: PTLNDIST

EXPORT PTLNDIST()

BEGIN

// 2024-07-21 EWS

// radian

HAngle:=0;

LOCAL px,py,m,b;

INPUT({m,b,px,py},

"Point-Line Distance (px, py), y=mx+b",

{"m:","b:","px:","py:"},

{"m: slope"," b: y-intercept",

"point x","point y"});

LOCAL y0;

y0:=m*px+b;

LOCAL b1,mt,x1,y1,dist,str;

IF m≠0 THEN

b1:=py+px/m;

mt:=[[−m,1],[1/m,1]]^-1*[[b],[b1]];

x1:=mt[1,1];

y1:=mt[2,1];

dist:=ABS((x1-px)+(y1-py)*√(-1));

ELSE

x1:=px;

y1:=b;

dist:=ABS((y1-py)*√(-1));

END;

// print results

PRINT();

PRINT("Results:");

PRINT("Intersect point:");

PRINT("x: "+STRING(x1));

PRINT("y: "+STRING(y1));

PRINT("");

IF m≠0 THEN

str:="Y="+STRING(-1/m)+"*X+"+STRING(b1);

ELSE

str:="X="+STRING(x1);

END;

PRINT("Orthogonal Line:");

PRINT(str);

PRINT("");

PRINT("Minimum Distance:");

PRINT(dist);

RETURN {x1,y1,str,dist};

END;

Note:

√(-1) represents the imaginary number ⅈ ( [ Shift ], [ 2 ] ).

Inputs:

* The slope of the y-intercept of the line y = m * x + b (no vertical lines, but m can be zero)

* The point (px, py)

Outputs:

* The line that runs through point (px, yx) that is orthogonal to y = m * x + b. The slope and y-intercept of the orthogonal line, which the line will be stated in a string

* The intersection point of the two lines.

* The distance between (px, yx) and the intersection point. (dist)

Examples

Example 1:

Inputs: Line: y = 5 x – 2, Point: (-1, -5)

m = 5

b = -2

px = -1

py = -5

Results:

Intersect point:

x = -0.615384615386

y = -5.07692307692

Orthogonal Line:

Y = -0.2 * X – 5.2

Minimum distance:

0.392232270274

Example 2:

Inputs: Line: y = 6, Point: (3, -9)

m = 0

b = 6

px = 3

py = -9

Results:

Intersect point:

x = 0.764705882353

y = 4.05882352941

Orthogonal Line:

X = 3

Minimum distance:

15

Example 3:

Inputs: Line: y = 4 x + 1, Point: (5, 3)

m = 4

b = 1

px = 5

py = 3

Results:

Intersect point:

x = 0.764705882353

y = 4.05882352941

Orthogonal Line:

Y = -0.25 * X + 4.25

Minimum distance:

4.36564125066

Source

Tremblay, Christopher. Mathematics for Game Developers. Thomson Course Technology. Boston, MA. 2004. ISBN 1-59200-038-X.

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

HP Prime Update: Official Firmware 15157

 HP has released an official firmware for the HP Prime:  Firmware 15157.  When you start the HP Connectivity Kit, the program should prompt you to download the firmware.

More information and downloads from TI-Planet (the page is in French):

https://tiplanet.org/forum/viewtopic.php?p=276332#p276332

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Casio fx-CG 50 and Python: Gaussian Quadrature

 Casio fx-CG 50 and Python:  Gaussian Quadrature

The Method

Gaussian quadrature is a method of numerical integration of a function f(x) over an interval [a, b] where a < b. Quadrature works similar to the Trapezoid Rule or Simpson's Rule, but instead of equally spaced points of measure, points are determined by the roots of a certain polynomial. Polynomials used include the Legendre, Jacobi, Chebyshev, Hermite, and Laguerre. Each point, xi, will be assigned a corresponding certain weight, wi, designed to make the approximation the most accurate as possible.

In general, ∫( f(x) dx, x = a to x = b) ≈ Σ( wi * f(xi) from i = 1 to n)

For today’s blog, the Legendre polynomial is used and the domain of integration is restricted to the interval [ -1, 1 ]. Hence:

∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n)

Steps:

1. Determine the number of points (order) is needed.

2. Find the roots of the Legendre polynomial, P_n(x) = 0. This determines the points xi.

3. Calculate the weights as follows:

wi = 2 / ((1 – xi^2) * (P’_n(xi))^2)

where P’_n(xi) is the numerical derivative of the Legendre polynomial at point x = xi.

4. Approximate ∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n).

Legendre Polynomials

There are several ways to determine the Legendre polynomials of order n, these are just a few:

P_n(x) = 1 / (2^n * n!) d^n/dx^n (x^2 – 1)^n

P_n(x) = Σ( comb(n, k), * comb(n + k, k) * comb(x -1, 2)^k for k = 0 to n)

where comb is the combination method

Recursive method:

(n + 1) * P_n(x) = (2*n + 1) * x * P_n-1(x) – n * P_n-2(x)

where P_0(x) = 1 and P_1(x) = x.

A table of Legendre polynomials of orders 0 to 10 is found here (scroll down a bit):

https://en.wikipedia.org/wiki/Legendre_polynomials

An Example: Order 3

Let’s determine the points xi and the weights wi for order 3 (n = 3).

The Legendre polynomial of order 3 is:

P_3(x) = 1 / 2 * (5 * x^3 – 3 * x) = 5 / 2 * x^3 – 3 / 2 * x

And it’s derivative:

P’_3(x) = 15 / 2 * x^2 – 3 / 2

Solve P_3(x) = 0 yields:

5 / 2 * x^3 – 3 / 2 * x = 0

x * (5 / 2 * x^2 – 3 / 2) = 0

x= 0 or 5 / 2 * x^2 – 3 / 2 = 0:

5 / 2 * x^2 – 3 / 2 = 0

5 / 2 * x^2 = 3 / 2

x^2 = 3 / 5

x = ±√(3 / 5)

The roots are, from lowest value to highest value are: x1 = - √(3 / 5), x2 = 0, x3 = √(3 / 5)

Determining the associated weights yields: w1 = 5 / 9 , w2 = 8 / 9, w3 = 5 / 9

Our points and weights are:

x1 = -√(3 / 5), w1 = 5 / 9

x2 = 0, w2 = 8 / 9

x3 = √(3 / 5), w3 = 5 / 9

Casio fx-CG 50 Program: QUADRAT0

(The last character in the file name is a zero)

“3PT QUADRATURE”

“-1 TO 1”

“NO QUOTES”

“F(X)”? → fn1

∫(fn1, A, B) → D

{-√(3 _| 5), 0, √(3 _| 5) } → List 1

{ 5 _| 9, 8 _| 9, 5 _| 9 } → List 2

0 → S

For 1 → I To 3

List 1[ I ] → X

List 2[ I ] → W

S + W × fn1 → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

_| is the fraction symbol, pressed by the [ [] / [] ] key.

If you have a monochrome calculator, such as the fx-9750G series or the fx-9860G series, leave out the color commands (Black, Blue, Red) in front of the Locate commands.

The function is stored in function memory, slot 1, which is found in the OPTN menu.

Variables:

D = actual integral, using the integral command of the calculators

S = approximation determined the quadrature

Examples

Radians angle assumed	D (actual)	S (approximate)
∫( sin x dx, x = -1 to x = 1)	0	0
∫( e^(.5 * x) dx, x = -1 to x = 1)	2.084381222	2.084380222
∫( ln(x + 4) / 2 dx, x = -1 to x – 1)	1.365058335	1.365060354

7^th Order and Expanding the Range

The 7^th order Legendre polynomial is:

P_7 = 1 / 16 * (429 * x^7 – 693 * x^5 + 315 * x^3 – 35 * x)

The associated roots (xi) and weights (wi) are:

Point, xi =	Weight, wi =
-0.949107912342759	0.129484966168870
-0.741531185599394	0.279705391489277
-0.405845151377397	0.381830050505119
0	0.417959183673469
0.405845151377397	0.381830050505119
0.741531185599394	0.279705391489277
0.949107912342759	0.129484966168870

These values are from the table presented in “Gauss-Kronrod quadrature formula” (see Source). I verified the roots with Wolfram Alpha, which you can see here:

https://www.wolframalpha.com/input?i2d=true&i=Divide%5B1%2C16%5D+*++%5C%2840%29429+*+Power%5Bx%2C7%5D+%E2%80%93+693+*+Power%5Bx%2C5%5D+%2B+315+*+Power%5Bx%2C3%5D+%E2%80%93+35+*+x%5C%2841%29%3D0

We are not limited to the interval [-1, 1]. To use any interval [a, b], we can scale the points as so:

xi ‘ = (b – a) / 2 * xi + (b + a) / 2

wi’ = (b – a) / 2 * wi

Casio fx-CG50 Program: GAUSS7

This the 7-point Gaussian quadrature, the calculator basic version.

“7PT GAUSS RULE”

“NO QUOTES”

“F(X)”? → fn1

“A”? → A

“B”? → B

∫(fn1, A, B) → D

{ -0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759 } → List 1

{ 0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870 } → List 2

0 → S

For 1 → I To 7

(B – A) ÷ 2 × List 1[ I ] + (B + A) ÷ 2 → X

List 2[ I ] × (B – A) ÷ 2 × fn1 + S → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

Python Code: gauss7.py

The only module needed is the math module, so it should work with all calculators. This script was created with the fx-CG 50.

from math import *

# 7 point Gauss Integral Approximation

def f(x):

return (x-1/2)*(x+1/4)

# xi and wi

xi=[-0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759]

wi=[0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870]

# sum

s=0

# limits

a=eval(input("a? "))

b=eval(input("b? "))

# integral approximation

for i in range(7):

  x=(b-a)/2*xi[i]+(b+a)/2

  w=(b-a)/2*wi[i]

  s+=w*f(x)

print("Approx. Integral:\n",s)

Examples

Radians mode	D (actual)	S (Casio Basic)	S (Python)
∫( 1 / √(x + 3) dx, 1, 5)	1.656854249	1.658042205	1.656854249504868
∫( .03 * e^(-x^2) dx, -1, 1)	0.04480944797	0.04484387349	0.04480944866632411
∫( e^(-x) dx, -5, 0)	147.4131591	147.5268171	147.4131590824608
∫( .3 * ( cos x – 0.05) dx, 0, π / 2)	0.2764380551	0.2767068146	0.2764380551025117
∫( (x – 1 / 2) * (x + 1 / 4) dx, -1, 1)	0.4166666667	0.4168576653	0.41666666666867568

Notice a difference between the Basic and Python?

Observations

* The program calls for approximations of xi and wi that are carried out 15 digits. Since the roots of the Legendre polynomials are most likely to be irrational, the more digits used, the better (generally).

* Most calculators will store up to 10 to 13 digits internally.

* By contrast, Python will store many more decimal places, at least 15 to 17.

* The approximation uses a lot of calculations. Round off errors can affect the final results, even when the round off is small.

* I recommend the Quadrature method if you are using Python or any other platform that allows for a lot of decimal places.

* Thankfully, most calculators have the numeric integration function.

Sources

“Guass-Kronrod quadrature formula” Wikipedia. Last Edited December 27, 2023. Accessed July 3, 2024. https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula

“Gaussian quadrature” Wikipedia. Last Edited July 7, 2024. Accessed July 7, 2024. https://en.wikipedia.org/wiki/Gaussian_quadrature

This is a long one, folks. Thank you and have a great day. To all the students, may your year of studies be full of success, discovery, and joy.

Until next time,

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

TI 30Xa Algorithm: The Sigmoid Function and its Integral

 TI 30Xa Algorithm: The Sigmoid Function and its Integral

The Sigmoid Function

The sigmoid function is defined as:

S = 1 / (1 + e^(-x))

If we multiply the last equation by e^(x) / e^(x), we get:

S = e^(x) / (e^(x) + 1) = e^(x) / (1 + e^(x))

The sigmoid function takes functions from the real numbers and maps them to the interval (0,1). If x → -∞, S → 0. If x → +∞, S → 1.

The inverse of the sigmoid function is the logit function, from which we can derive:

S = 1 / (1 + e^(-x))

1 / S = 1 + e^(-x)

1 / S – 1 = e^(-x)

ln( 1 / S – 1 ) = -x

ln( (1 – S) / S ) = -x (Note: 1 / S – 1 = 1 / S – S / S = (1 – S) / S )

-ln( (1 – S) / S ) = x

ln( S / (1 – S)) = x (Note: For any x, -ln(x) = ln(x^(-1)) = ln(1 / x))

The rest of the blog will focus on the sigmoid function.

The Integral of the Sigmoid Function

Funding the area of the curve under the sigmoid function is pretty straight forward.

∫ 1 / (1 + e^(-x)) dx

= ∫ e^(x) / (e^(x) + 1) dx

Let u = e^(x) + 1. The du = e^(x) dx and:

∫ du / u

= ln | u | + C

= ln (e^(x) +1) + C (since e^(x) + 1 > 0 for all real x)

The definite integral (area) can be calculated as:

x = b

∫ e^(x) / (e^(x) + 1) dx

x = a

= ln (1 + e^b) – ln (1 + e^a)

What if we are given the area, A (capital A), from negative infinity to a value x? Here we are finding the lower tail area.

t = x

∫ e^(t) / (e^(t) + 1) dx = A

t = -∞

ln(1 + e^(x)) – ln(1 + e^(-∞)) = A

Note that e^(-∞) → 0 as t → -∞

We can estimate that:

ln(1 + e^(x)) – ln( 1 ) = A

ln(1 + e^(x)) – 0 = A

1 + e^(x) = e^(A)

e^(x) = e^(A) – 1

x = ln(e^(A) – 1)

TI-30Xa Algorithms

Sigmoid function:

S = 1 / (1 + e^(-x)) = e^(x) / (e^(x) + 1)

Algorithm:

[ ( ] x [ +/- ] [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ 1/x ] [ = ]

Example:

S(1.5):

[ ( ] 1.5 [ +/- ] [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ 1/x ] [ = ]

Result: 0.817574476

Integral of the sigmoid function:

x = b

∫ e^(x) / (e^(x) + 1) dx

x = a

Algorithm:

[ ( ] b [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ]

[ - ] [ ( ] a [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ] [ = ]

Example:

x = 3

∫ e^(x) / (e^(x) + 1) dx

x = 0

(a = 0, b = 3)

[ ( ] 3 [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ]

[ - ] [ ( ] 0 [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ] [ = ]

Result: 2.355440171

Find x given lower tail area:

x = ln(e^(A) – 1)

Algorithm:

[ ( ] A [ 2^nd ] [ LN ] {e^x} [ - ] 1 [ ) ] [ LN ] [ = ]

Example:

Area: A = 0.5

[ ( ] 0.5 [ 2^nd ] [ LN ] {e^x} [ - ] 1 [ ) ] [ LN ] [ = ]

Result: -0.4327521296

Until next time,

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.