Python Numbers: A Deep Dive Into Numeric Data Types

# Creating integers
age = 25                  # Direct assignment
population = 7_800_000_000  # Using underscores for readability (Python 3.6+)
hex_value = 0xFF          # Hexadecimal (255)
octal_value = 0o100       # Octal (64)
binary_value = 0b1010     # Binary (10)

# Converting other types to integers
count = int(25.7)         # Converting from float (truncates to 25)
digit_str = int("42")     # Converting from string (42)
bool_val = int(True)      # Boolean to int (1)
empty_int = int()         # Creates zero

print(f"Age: {age}")                    # => 25
print(f"Population: {population}")      # => 7800000000
print(f"Hex to decimal: {hex_value}")   # => 255
print(f"Binary to decimal: {binary_value}")  # => 10

# Python's unlimited integer size
massive_number = 2**1000  # 2 raised to the power of 1000
print(f"Digits in 2^1000: {len(str(massive_number))}")  # => 302 digits
                    

# Creating floating point numbers
temperature = 98.6           # Direct decimal assignment
precise = float(28)          # Converting from integer (28.0)
scientific = 2.5e-3          # Scientific notation (0.0025)
negative_float = -3.14       # Negative floating point

# Various ways to create floats
from_division = 10 / 3       # Division of integers yields float (3.3333...)
from_string = float("42.5")  # Converting from string
nan_value = float("nan")     # Not a Number (NaN)
infinity = float("inf")      # Positive infinity

print(f"Temperature: {temperature}")                 # => 98.6
print(f"Scientific notation: {scientific}")          # => 0.0025
print(f"Division result: {from_division:.4f}")       # => 3.3333 (formatted)
print(f"Is infinity > 1000000?: {infinity > 1000000}")  # => True

# Understanding floating point precision
print(f"0.1 + 0.2 = {0.1 + 0.2}")                   # => 0.30000000000000004
print(f"Is 0.1 + 0.2 == 0.3? {0.1 + 0.2 == 0.3}")   # => False
print(f"Rounded: {round(0.1 + 0.2, 1)}")            # => 0.3
                    

# Very small number
tiny = 0.0000000123
print(tiny)                    # => 1.23e-08

# Very large number
huge = 6700000000.0
print(huge)                    # => 6.7e+09

# Explicit scientific notation in code
planck_length = 1.616255e-35   # Physics constant
avogadro = 6.02214076e23       # Chemistry constant

# Formatting scientific notation for display
print(f"{avogadro:.2e}")       # => 6.02e+23 (formatted for readability)
                    

# Creating complex numbers
z1 = 3 + 4j              # Direct assignment
z2 = complex(1, 2)       # Using constructor: 1+2j

# Accessing parts of complex numbers
print(f"z1 = {z1}")                 # => (3+4j)
print(f"Real part: {z1.real}")      # => 3.0
print(f"Imaginary part: {z1.imag}") # => 4.0

# Basic operations with complex numbers
z3 = z1 + z2             # Addition: (4+6j)
z4 = z1 * z2             # Multiplication: (3+4j)*(1+2j) = (3-8+4j+8j²) = (3-8+4j-8) = (-5+4j)
modulus = abs(z1)        # Magnitude: sqrt(3² + 4²) = sqrt(25) = 5.0

print(f"z1 + z2 = {z3}")
print(f"z1 * z2 = {z4}")
print(f"Modulus of z1: {modulus}")  # Distance from origin in complex plane

# Practical application: Polar conversion
import math
r = abs(z1)
theta = math.atan2(z1.imag, z1.real)
print(f"Polar form: {r}∠{theta:.2f} radians")
                    

# Order of operations examples
result1 = 2 + 3 * 4       # Multiplication happens first
print(f"2 + 3 * 4 = {result1}")         # 14, not 20

result2 = (2 + 3) * 4     # Parentheses take precedence
print(f"(2 + 3) * 4 = {result2}")       # 20

# More complex example
result3 = 2 ** 3 * 4 + 5  # Exponent, then multiplication, then addition
print(f"2 ** 3 * 4 + 5 = {result3}")    # 37

# When in doubt, use parentheses for clarity
complex_calculation = ((10 + 5) * 3 - 4) / 2
print(f"((10 + 5) * 3 - 4) / 2 = {complex_calculation}")  # 20.5
                    

from decimal import Decimal, getcontext

# The floating point problem
print(f"0.1 + 0.2 = {0.1 + 0.2}")                  # 0.30000000000000004
print(f"Is 0.1 + 0.2 == 0.3? {0.1 + 0.2 == 0.3}")  # False

# The decimal solution
a = Decimal('0.1')
b = Decimal('0.2')
c = a + b
print(f"Decimal 0.1 + 0.2 = {c}")                 # 0.3
print(f"Is Decimal 0.1 + 0.2 == 0.3? {c == Decimal('0.3')}")  # True

# Setting precision (defaults to 28 places)
getcontext().prec = 50
result = Decimal(1) / Decimal(3)
print(f"1/3 with 50 digits of precision: {result}")

# Financial calculation example
def calculate_compound_interest(principal, rate, years, compounds_per_year=1):
    """Calculate compound interest with precise decimal arithmetic."""
    p = Decimal(str(principal))
    r = Decimal(str(rate))
    n = Decimal(str(compounds_per_year))
    t = Decimal(str(years))
    
    # Formula: P(1 + r/n)^(nt)
    amount = p * (1 + r/n) ** (n * t)
    return amount

# Calculate 5% interest compounded monthly for 10 years on $1000
investment = calculate_compound_interest(1000, 0.05, 10, 12)
print(f"Investment after 10 years: ${investment:.2f}")  # $1647.01

# Contrast with float calculation
float_result = 1000 * (1 + 0.05/12) ** (12 * 10)
print(f"Float calculation: ${float_result:.2f}")        # $1647.01
print(f"Difference: ${float_result - float(investment):.15f}")  # Very small difference
                    

from fractions import Fraction

# Creating fractions
half = Fraction(1, 2)
third = Fraction(1, 3)
from_decimal = Fraction('0.25')
from_string = Fraction('3/4')

print(f"Half: {half}")              # 1/2
print(f"Third: {third}")            # 1/3
print(f"From decimal: {from_decimal}")  # 1/4
print(f"From string: {from_string}")   # 3/4

# Arithmetic with fractions
sum_frac = half + third
print(f"1/2 + 1/3 = {sum_frac}")     # 5/6

product = half * from_string
print(f"1/2 * 3/4 = {product}")      # 3/8

# Converting between forms
print(f"0.25 as a fraction: {Fraction('0.25')}")  # 1/4
print(f"1/3 as a float: {float(third)}")         # 0.3333333333333333
print(f"1/3 as a decimal: {Decimal(third.numerator) / Decimal(third.denominator)}")

# Accessing parts of a fraction
print(f"Numerator of 3/4: {from_string.numerator}")    # 3
print(f"Denominator of 3/4: {from_string.denominator}")  # 4

# Example: Calculating probabilities
def coin_toss_probability(heads, total_tosses):
    """Return the probability of getting exactly 'heads' in 'total_tosses'."""
    from math import comb
    
    # The probability as a Fraction
    probability = Fraction(comb(total_tosses, heads), 2**total_tosses)
    return probability

# Probability of exactly 3 heads in 5 tosses
prob = coin_toss_probability(3, 5)
print(f"Probability of 3 heads in 5 tosses: {prob}")  # 10/32 (simplifies to 5/16)
print(f"As a decimal: {float(prob):.4f}")            # 0.3125
                    

# Division by zero
try:
    result = 10 / 0
except ZeroDivisionError as e:
    print(f"Error: {e}")  # Error: division by zero

# Better approach: check before dividing
divisor = 0
if divisor != 0:
    result = 10 / divisor
else:
    result = float('inf')  # or some other appropriate value
    print("Used infinity as result for division by zero")

# Integer division surprises
print(f"5 / 2 = {5 / 2}")     # 2.5 (normal division)
print(f"5 // 2 = {5 // 2}")   # 2 (floor division)

# Negative number division
print(f"-5 / 2 = {-5 / 2}")     # -2.5 (normal division)
print(f"-5 // 2 = {-5 // 2}")   # -3 (floor division, rounds down)
                    

# Incorrect way to compare floats
x = 0.1 + 0.2
y = 0.3
print(f"Is x == y? {x == y}")  # False! (0.30000000000000004 != 0.3)

# Better approach 1: Use round for display purposes
print(f"Rounded x: {round(x, 1)}")  # 0.3

# Better approach 2: Check if difference is very small
epsilon = 1e-10  # A very small number
print(f"Is x almost equal to y? {abs(x - y) < epsilon}")  # True

# Best approach 3: Use isclose from math module (Python 3.5+)
import math
print(f"Is x close to y? {math.isclose(x, y)}")  # True
                    

# String to number conversions
try:
    num = int("3.14")  # This fails!
except ValueError as e:
    print(f"Error: {e}")  # Error: invalid literal for int() with base 10: '3.14'

# Correct approach
num = int(float("3.14"))  # 3

# Unexpected string concatenation
value = 5
print("The value is " + str(value))  # Need explicit conversion for concatenation

# Safer approach: use f-strings
print(f"The value is {value}")  # No explicit conversion needed

# Implicit type conversion in expressions
result = 5 + 3.14  # Integer and float mixed
print(f"5 + 3.14 = {result}")  # 8.14 (converts to float)
print(f"Type of result: {type(result)}")  # <class 'float'>
                    

# Boolean value coercion
true_value = True
print(f"True as an integer: {int(true_value)}")  # 1
print(f"True + True: {true_value + true_value}")  # 2

# None handling
try:
    result = 5 + None  # This is an error
except TypeError as e:
    print(f"Error: {e}")  # Error: unsupported operand type(s) for +: 'int' and 'NoneType'

# Better approach: check for None
value = None
if value is not None:
    result = 5 + value
else:
    result = 5
    print("Used default value because input was None")
                    

from decimal import Decimal, getcontext

# Set precision for financial calculations
getcontext().prec = 28

class FinancialCalculator:
    """A calculator for various financial operations with high precision."""
    
    @staticmethod
    def compound_interest(principal, rate, time, compounds_per_year=1):
        """
        Calculate compound interest.
        
        Args:
            principal: Initial investment amount
            rate: Annual interest rate (decimal, e.g. 0.05 for 5%)
            time: Time in years
            compounds_per_year: Number of compounds per year
            
        Returns:
            Final amount after interest
        """
        p = Decimal(str(principal))
        r = Decimal(str(rate))
        t = Decimal(str(time))
        n = Decimal(str(compounds_per_year))
        
        # Formula: P(1 + r/n)^(nt)
        amount = p * (1 + r/n) ** (n * t)
        return amount
    
    @staticmethod
    def mortgage_payment(principal, rate, years):
        """
        Calculate monthly mortgage payment.
        
        Args:
            principal: Loan amount
            rate: Annual interest rate (decimal)
            years: Loan term in years
            
        Returns:
            Monthly payment amount
        """
        p = Decimal(str(principal))
        r = Decimal(str(rate)) / 12  # Monthly rate
        n = Decimal(str(years)) * 12  # Number of payments
        
        # Formula: P * (r * (1 + r)^n) / ((1 + r)^n - 1)
        payment = p * (r * (1 + r) ** n) / ((1 + r) ** n - 1)
        return payment
    
    @staticmethod
    def present_value(future_value, rate, time):
        """
        Calculate present value of a future amount.
        
        Args:
            future_value: Value at the end of the period
            rate: Annual discount rate (decimal)
            time: Time in years
            
        Returns:
            Present value
        """
        fv = Decimal(str(future_value))
        r = Decimal(str(rate))
        t = Decimal(str(time))
        
        # Formula: FV / (1 + r)^t
        pv = fv / (1 + r) ** t
        return pv

# Example usage
calculator = FinancialCalculator()

# Compound interest example
investment = 10000
rate = 0.05  # 5%
years = 10
compounds = 12  # monthly compounding

final_amount = calculator.compound_interest(investment, rate, years, compounds)
print(f"Investment of ${investment} after {years} years at {rate*100}% interest: ${final_amount:.2f}")

# Mortgage calculation example
home_price = 300000
down_payment = 60000
loan_amount = home_price - down_payment
mortgage_rate = 0.045  # 4.5%
loan_term = 30  # years

monthly_payment = calculator.mortgage_payment(loan_amount, mortgage_rate, loan_term)
print(f"Monthly mortgage payment on ${loan_amount} at {mortgage_rate*100}% for {loan_term} years: ${monthly_payment:.2f}")

# Present value example
retirement_goal = 1000000
expected_return = 0.07  # 7%
time_to_retirement = 25  # years

needed_now = calculator.present_value(retirement_goal, expected_return, time_to_retirement)
print(f"To have ${retirement_goal} in {time_to_retirement} years, you need ${needed_now:.2f} now with {expected_return*100}% returns.")
                    

import statistics
import math

class TemperatureConverter:
    """A utility for temperature conversions and analysis."""
    
    @staticmethod
    def celsius_to_fahrenheit(celsius):
        """Convert Celsius to Fahrenheit."""
        return (celsius * 9/5) + 32
    
    @staticmethod
    def fahrenheit_to_celsius(fahrenheit):
        """Convert Fahrenheit to Celsius."""
        return (fahrenheit - 32) * 5/9
    
    @staticmethod
    def celsius_to_kelvin(celsius):
        """Convert Celsius to Kelvin."""
        return celsius + 273.15
    
    @staticmethod
    def kelvin_to_celsius(kelvin):
        """Convert Kelvin to Celsius."""
        return kelvin - 273.15
    
    @staticmethod
    def fahrenheit_to_kelvin(fahrenheit):
        """Convert Fahrenheit to Kelvin."""
        celsius = TemperatureConverter.fahrenheit_to_celsius(fahrenheit)
        return TemperatureConverter.celsius_to_kelvin(celsius)
    
    @staticmethod
    def kelvin_to_fahrenheit(kelvin):
        """Convert Kelvin to Fahrenheit."""
        celsius = TemperatureConverter.kelvin_to_celsius(kelvin)
        return TemperatureConverter.celsius_to_fahrenheit(celsius)
    
    @classmethod
    def analyze_temperatures(cls, temperatures, unit='C'):
        """
        Analyze a list of temperatures.
        
        Args:
            temperatures: List of temperature values
            unit: The unit of the input temperatures ('C', 'F', or 'K')
            
        Returns:
            Dictionary with statistics and all temperature representations
        """
        # Convert all temperatures to all units
        all_units = []
        
        for temp in temperatures:
            if unit == 'C':
                celsius = temp
                fahrenheit = cls.celsius_to_fahrenheit(temp)
                kelvin = cls.celsius_to_kelvin(temp)
            elif unit == 'F':
                fahrenheit = temp
                celsius = cls.fahrenheit_to_celsius(temp)
                kelvin = cls.celsius_to_kelvin(celsius)
            elif unit == 'K':
                kelvin = temp
                celsius = cls.kelvin_to_celsius(temp)
                fahrenheit = cls.celsius_to_fahrenheit(celsius)
            else:
                raise ValueError("Unit must be 'C', 'F', or 'K'")
                
            all_units.append({
                'C': celsius,
                'F': fahrenheit,
                'K': kelvin
            })
        
        # Calculate statistics on the original unit
        stats = {
            'mean': statistics.mean(temperatures),
            'median': statistics.median(temperatures),
            'stdev': statistics.stdev(temperatures) if len(temperatures) > 1 else 0,
            'min': min(temperatures),
            'max': max(temperatures),
            'range': max(temperatures) - min(temperatures)
        }
        
        return {
            'original_unit': unit,
            'statistics': stats,
            'all_conversions': all_units
        }

# Example usage
converter = TemperatureConverter()

# Simple conversion examples
freezing_c = 0
freezing_f = converter.celsius_to_fahrenheit(freezing_c)
freezing_k = converter.celsius_to_kelvin(freezing_c)

print(f"Water freezes at {freezing_c}°C, {freezing_f}°F, or {freezing_k}K")

# Temperature analysis example
daily_highs_f = [78.3, 82.1, 77.6, 83.4, 80.2, 81.9, 76.5]
analysis = converter.analyze_temperatures(daily_highs_f, unit='F')

print("\nDaily Temperature Analysis:")
print(f"Original unit: {analysis['original_unit']}")
print("\nStatistics:")
for key, value in analysis['statistics'].items():
    print(f"  {key}: {value:.2f}°{analysis['original_unit']}")

print("\nConversions (first 3 days):")
for i, day in enumerate(analysis['all_conversions'][:3], 1):
    print(f"  Day {i}: {day['C']:.2f}°C, {day['F']:.2f}°F, {day['K']:.2f}K")
                    

import math
from collections import Counter

class StatisticsEngine:
    """Implementation of basic statistical functions."""
    
    @staticmethod
    def mean(data):
        """Calculate arithmetic mean."""
        if not data:
            return None
        return sum(data) / len(data)
    
    @staticmethod
    def median(data):
        """Calculate median value."""
        if not data:
            return None
            
        sorted_data = sorted(data)
        n = len(sorted_data)
        
        if n % 2 == 0:
            # Even number of items
            middle1 = sorted_data[n // 2 - 1]
            middle2 = sorted_data[n // 2]
            return (middle1 + middle2) / 2
        else:
            # Odd number of items
            return sorted_data[n // 2]
    
    @staticmethod
    def mode(data):
        """Find the most common value(s)."""
        if not data:
            return None
            
        # Count occurrences of each value
        counts = Counter(data)
        max_count = max(counts.values())
        
        # Return all values that occur max_count times
        return [value for value, count in counts.items() if count == max_count]
    
    @staticmethod
    def variance(data, sample=True):
        """
        Calculate variance.
        
        Args:
            data: List of numeric values
            sample: If True, calculates sample variance (n-1 denominator),
                   otherwise population variance (n denominator)
        """
        if not data or len(data) < (2 if sample else 1):
            return None
            
        # Calculate mean
        data_mean = StatisticsEngine.mean(data)
        
        # Sum of squared differences from the mean
        squared_diff_sum = sum((x - data_mean) ** 2 for x in data)
        
        # Divide by n-1 for sample variance, or n for population variance
        denominator = len(data) - 1 if sample else len(data)
        return squared_diff_sum / denominator
    
    @staticmethod
    def std_dev(data, sample=True):
        """Calculate standard deviation."""
        variance = StatisticsEngine.variance(data, sample)
        if variance is None:
            return None
        return math.sqrt(variance)
    
    @staticmethod
    def covariance(data_x, data_y, sample=True):
        """Calculate covariance between two datasets."""
        if len(data_x) != len(data_y) or len(data_x) < (2 if sample else 1):
            return None
        
        mean_x = StatisticsEngine.mean(data_x)
        mean_y = StatisticsEngine.mean(data_y)
        
        # Sum of products of differences from means
        sum_product_diff = sum((x - mean_x) * (y - mean_y) for x, y in zip(data_x, data_y))
        
        # Divide by n-1 for sample covariance, or n for population covariance
        denominator = len(data_x) - 1 if sample else len(data_x)
        return sum_product_diff / denominator
    
    @staticmethod
    def correlation(data_x, data_y):
        """Calculate Pearson correlation coefficient."""
        if len(data_x) != len(data_y) or len(data_x) < 2:
            return None
        
        # Calculate covariance
        cov = StatisticsEngine.covariance(data_x, data_y)
        
        # Calculate standard deviations
        std_x = StatisticsEngine.std_dev(data_x)
        std_y = StatisticsEngine.std_dev(data_y)
        
        # Correlation is covariance divided by product of standard deviations
        if std_x == 0 or std_y == 0:
            return 0  # No correlation if either dataset has no variation
        
        return cov / (std_x * std_y)
    
    @staticmethod
    def summarize(data):
        """Generate a complete statistical summary."""
        if not data:
            return {"error": "Empty dataset"}
        
        return {
            "count": len(data),
            "mean": StatisticsEngine.mean(data),
            "median": StatisticsEngine.median(data),
            "mode": StatisticsEngine.mode(data),
            "range": max(data) - min(data),
            "variance": StatisticsEngine.variance(data),
            "std_dev": StatisticsEngine.std_dev(data),
            "min": min(data),
            "max": max(data)
        }

# Example usage
stats = StatisticsEngine()

# Test data
heights = [165, 170, 168, 172, 175, 173, 168, 170, 171]
weights = [62, 65, 63, 70, 72, 68, 63, 64, 66]

# Basic statistics
print(f"Mean height: {stats.mean(heights):.2f} cm")
print(f"Median height: {stats.median(heights):.2f} cm")
print(f"Mode height: {stats.mode(heights)} cm")
print(f"Standard deviation: {stats.std_dev(heights):.2f} cm")

# Relationship between height and weight
correlation = stats.correlation(heights, weights)
print(f"Correlation between height and weight: {correlation:.4f}")

# Complete summary
summary = stats.summarize(heights)
print("\nComplete Height Statistics:")
for key, value in summary.items():
    if isinstance(value, float):
        print(f"  {key}: {value:.2f}")
    else:
        print(f"  {key}: {value}")
                    

import cmath  # The complex math module

# Creating complex numbers
z1 = 3 + 4j
z2 = complex(2, -1)  # 2-1j

# Basic properties
print(f"z1 = {z1}")
print(f"Real part: {z1.real}")
print(f"Imaginary part: {z1.imag}")
print(f"Conjugate: {z1.conjugate()}")  # 3-4j

# Complex arithmetic
print(f"z1 + z2 = {z1 + z2}")         # 5+3j
print(f"z1 - z2 = {z1 - z2}")         # 1+5j
print(f"z1 * z2 = {z1 * z2}")         # 10+5j
print(f"z1 / z2 = {z1 / z2}")         # 0.4+2.2j

# Polar form
r = abs(z1)                          # Magnitude (modulus)
theta = cmath.phase(z1)              # Phase angle in radians
print(f"Polar form: {r}∠{theta:.2f} radians")
print(f"Polar form: {r}∠{theta * 180/cmath.pi:.2f} degrees")

# Converting between rectangular and polar forms
# Polar to rectangular: r * (cos θ + j sin θ)
r, theta = 5, cmath.pi/4
z_rect = r * (math.cos(theta) + 1j * math.sin(theta))
print(f"Polar to rectangular: {z_rect}")

# Or using cmath.rect
z_rect2 = cmath.rect(r, theta)
print(f"Using cmath.rect: {z_rect2}")

# Complex functions
print(f"Square root of -1: {cmath.sqrt(-1)}")           # 1j
print(f"Natural log of j: {cmath.log(1j)}")             # 1.5707963267948966j
print(f"e^(jπ): {cmath.exp(1j * cmath.pi)}")           # (-1+1.2246467991473532e-16j) ≈ -1
print(f"cos(j): {cmath.cos(1j)}")                      # 1.5430806348152437

# Practical application: Impedance in electrical circuits
class ComplexImpedance:
    """Calculate impedance in AC circuits."""
    
    @staticmethod
    def resistor(resistance):
        """Resistor impedance: purely real."""
        return complex(resistance, 0)
    
    @staticmethod
    def capacitor(capacitance, frequency):
        """Capacitor impedance: negative imaginary component."""
        if capacitance == 0 or frequency == 0:
            return complex(0, float('inf'))
        reactance = -1 / (2 * math.pi * frequency * capacitance)
        return complex(0, reactance)
    
    @staticmethod
    def inductor(inductance, frequency):
        """Inductor impedance: positive imaginary component."""
        reactance = 2 * math.pi * frequency * inductance
        return complex(0, reactance)
    
    @staticmethod
    def series(*impedances):
        """Impedances in series are added."""
        total = complex(0, 0)
        for z in impedances:
            total += z
        return total
    
    @staticmethod
    def parallel(*impedances):
        """Impedances in parallel: reciprocal of sum of reciprocals."""
        reciprocal_sum = complex(0, 0)
        for z in impedances:
            if z == complex(0, 0):
                return complex(0, 0)  # Short circuit
            reciprocal_sum += 1 / z
        return 1 / reciprocal_sum

# Example: RLC circuit analysis
R = 100  # 100 ohms
L = 0.1  # 0.1 henry
C = 1e-6  # 1 microfarad
f = 1000  # 1000 Hz

z_resistor = ComplexImpedance.resistor(R)
z_inductor = ComplexImpedance.inductor(L, f)
z_capacitor = ComplexImpedance.capacitor(C, f)

# Series RLC circuit
z_series = ComplexImpedance.series(z_resistor, z_inductor, z_capacitor)
print(f"\nRLC Series Circuit at {f} Hz:")
print(f"Resistor: {z_resistor} ohms")
print(f"Inductor: {z_inductor} ohms")
print(f"Capacitor: {z_capacitor} ohms")
print(f"Total impedance: {z_series} ohms")
print(f"Magnitude: {abs(z_series):.2f} ohms")
print(f"Phase angle: {cmath.phase(z_series) * 180/cmath.pi:.2f} degrees")
                    

import decimal
from decimal import Decimal
from fractions import Fraction
import math

# Setting precision for decimal calculations
decimal.getcontext().prec = 50  # 50 digits of precision

# Calculating pi to high precision
def calculate_pi_archimedes(iterations):
    """
    Calculate π using Archimedes' method.
    Converges slowly but demonstrates arbitrary precision.
    """
    # Start with a hexagon inscribed in a circle
    sides = 6
    # Initial approximation based on hexagon
    a = Decimal(1)  # Length of one side
    s = Decimal(sides * a / 2)  # Perimeter / 2
    r = Decimal(1)  # Radius of circle
    
    # Iteratively double the number of sides
    for i in range(iterations):
        sides *= 2
        # Update the side length based on geometry
        a = Decimal(Decimal(2) - Decimal(Decimal(4) - a ** Decimal(2)).sqrt()) / Decimal(2)
        s = sides * a / 2
    
    return s

# Calculate pi to 30 decimal places
pi_approx = calculate_pi_archimedes(20)
print(f"π calculated to 50 digits: {pi_approx}")
print(f"Python's math.pi: {math.pi}")
print(f"Difference: {abs(pi_approx - Decimal(math.pi))}")

# Working with extremely large numbers
factorial_100 = math.factorial(100)
print(f"100! has {len(str(factorial_100))} digits")
print(f"First 50 digits of 100!: {str(factorial_100)[:50]}...")

# Computing with fractions for exact results
def compute_sum_of_fractions(n):
    """Calculate sum of 1/i for i from 1 to n as an exact fraction."""
    result = Fraction(0, 1)
    for i in range(1, n + 1):
        result += Fraction(1, i)
    return result

# The harmonic series
harmonic_10 = compute_sum_of_fractions(10)
print(f"Sum of 1/i for i from 1 to 10: {harmonic_10}")
print(f"As decimal: {float(harmonic_10):.10f}")

# Calculating precise exponential function
def precise_exp(x, terms=50):
    """
    Calculate e^x with arbitrary precision using the Taylor series.
    """
    result = Decimal(1)  # First term is 1
    term = Decimal(1)
    
    for i in range(1, terms):
        # Each term is previous term * x / i
        term *= Decimal(x) / Decimal(i)
        result += term
        
        # Stop if term becomes too small to affect result
        if term < Decimal('1e-50'):
            break
    
    return result

# Calculate e to high precision
e_approx = precise_exp(1, 100)
print(f"e calculated to 50 digits: {e_approx}")
print(f"Python's math.e: {math.e}")

# Infinite precision square root
def precise_sqrt(n, iterations=20):
    """Calculate square root using Newton's method with arbitrary precision."""
    # Initial guess
    x = Decimal(n) / 2
    
    for i in range(iterations):
        # Newton's method: x' = (x + n/x) / 2
        x = (x + Decimal(n) / x) / 2
    
    return x

# Calculate square root of 2 to high precision
sqrt2 = precise_sqrt(2, 30)
print(f"√2 calculated to 50 digits: {sqrt2}")
print(f"Python's math.sqrt(2): {math.sqrt(2)}")
                    

import math

def trapezoidal_rule(func, a, b, n):
    """
    Numerically integrate a function using the trapezoidal rule.
    
    Args:
        func: Function to integrate
        a, b: Integration limits
        n: Number of trapezoids (intervals)
    
    Returns:
        Numerical approximation of the integral
    """
    h = (b - a) / n  # Width of each trapezoid
    result = 0.5 * (func(a) + func(b))  # End points
    
    for i in range(1, n):
        x = a + i * h
        result += func(x)
    
    return result * h

def simpson_rule(func, a, b, n):
    """
    Numerically integrate a function using Simpson's rule.
    
    Args:
        func: Function to integrate
        a, b: Integration limits
        n: Number of intervals (must be even)
    
    Returns:
        Numerical approximation of the integral
    """
    if n % 2 != 0:
        n += 1  # Ensure n is even
    
    h = (b - a) / n
    result = func(a) + func(b)  # End points
    
    for i in range(1, n):
        x = a + i * h
        weight = 4 if i % 2 == 1 else 2
        result += weight * func(x)
    
    return result * h / 3

def numerical_derivative(func, x, h=1e-5):
    """
    Calculate the numerical derivative of a function at a point.
    
    Args:
        func: Function to differentiate
        x: Point at which to calculate derivative
        h: Step size for finite difference
    
    Returns:
        Numerical approximation of the derivative
    """
    # Central difference formula: f'(x) ≈ [f(x+h) - f(x-h)] / (2h)
    return (func(x + h) - func(x - h)) / (2 * h)

def second_derivative(func, x, h=1e-4):
    """
    Calculate the numerical second derivative of a function at a point.
    
    Args:
        func: Function to differentiate
        x: Point at which to calculate derivative
        h: Step size for finite difference
    
    Returns:
        Numerical approximation of the second derivative
    """
    # Second derivative formula: f''(x) ≈ [f(x+h) - 2f(x) + f(x-h)] / h²
    return (func(x + h) - 2 * func(x) + func(x - h)) / (h**2)

# Test functions
def f1(x):
    """Simple polynomial: x^3"""
    return x**3

def f2(x):
    """Gaussian function: e^(-x²)"""
    return math.exp(-x**2)

def f3(x):
    """Sine function: sin(x)"""
    return math.sin(x)

# Integration examples
print("Integration Examples:")

# Integrate x³ from 0 to 1 (analytical result: 1/4 = 0.25)
result1_trap = trapezoidal_rule(f1, 0, 1, 100)
result1_simp = simpson_rule(f1, 0, 1, 100)
print(f"∫₀¹ x³ dx (Trapezoidal): {result1_trap:.8f}")
print(f"∫₀¹ x³ dx (Simpson): {result1_simp:.8f}")
print(f"Analytical result: 0.25")

# Integrate e^(-x²) from -1 to 1
result2_trap = trapezoidal_rule(f2, -1, 1, 100)
result2_simp = simpson_rule(f2, -1, 1, 100)
print(f"∫₋₁¹ e^(-x²) dx (Trapezoidal): {result2_trap:.8f}")
print(f"∫₋₁¹ e^(-x²) dx (Simpson): {result2_simp:.8f}")

# Integrate sin(x) from 0 to π (analytical result: 2)
result3_trap = trapezoidal_rule(f3, 0, math.pi, 100)
result3_simp = simpson_rule(f3, 0, math.pi, 100)
print(f"∫₀ᵗ sin(x) dx (Trapezoidal): {result3_trap:.8f}")
print(f"∫₀ᵗ sin(x) dx (Simpson): {result3_simp:.8f}")
print(f"Analytical result: 2")

# Differentiation examples
print("\nDifferentiation Examples:")

# Derivative of x³ at x=2 (analytical result: 3x² = 12)
deriv1 = numerical_derivative(f1, 2)
print(f"d/dx(x³) at x=2: {deriv1:.8f}")
print(f"Analytical result: 12")

# Derivative of e^(-x²) at x=0 (analytical result: -2xe^(-x²) = 0)
deriv2 = numerical_derivative(f2, 0)
print(f"d/dx(e^(-x²)) at x=0: {deriv2:.8f}")
print(f"Analytical result: 0")

# Second derivative of sin(x) at x=π/4 (analytical result: -sin(x) = -1/√2)
second_deriv = second_derivative(f3, math.pi/4)
print(f"d²/dx²(sin(x)) at x=π/4: {second_deriv:.8f}")
print(f"Analytical result: {-math.sin(math.pi/4):.8f}")

# Practical application: Numerical solution to a physical problem
def projectile_motion(angle_degrees, initial_velocity, x):
    """
    Calculate the height of a projectile at a given horizontal distance.
    
    Args:
        angle_degrees: Launch angle in degrees
        initial_velocity: Initial velocity in m/s
        x: Horizontal distance in meters
    
    Returns:
        Height of projectile at distance x
    """
    angle_rad = math.radians(angle_degrees)
    g = 9.81  # Gravitational acceleration in m/s²
    
    # Projectile motion equation: y = x*tan(θ) - (g*x²)/(2*v₀²*cos²(θ))
    term1 = x * math.tan(angle_rad)
    term2 = (g * x**2) / (2 * initial_velocity**2 * math.cos(angle_rad)**2)
    return term1 - term2

def find_maximum_range(angle_degrees, initial_velocity):
    """
    Calculate the maximum range of a projectile.
    
    Args:
        angle_degrees: Launch angle in degrees
        initial_velocity: Initial velocity in m/s
    
    Returns:
        Maximum horizontal distance
    """
    angle_rad = math.radians(angle_degrees)
    g = 9.81
    
    # Range formula: R = (v₀²*sin(2θ))/g
    return (initial_velocity**2 * math.sin(2 * angle_rad)) / g

def find_maximum_height(angle_degrees, initial_velocity):
    """
    Calculate the maximum height reached by a projectile.
    
    Args:
        angle_degrees: Launch angle in degrees
        initial_velocity: Initial velocity in m/s
    
    Returns:
        Maximum height
    """
    angle_rad = math.radians(angle_degrees)
    g = 9.81
    
    # Maximum height formula: h = (v₀²*sin²(θ))/(2g)
    return (initial_velocity**2 * math.sin(angle_rad)**2) / (2 * g)

# Projectile motion example
angle = 45  # degrees
velocity = 20  # m/s

max_range = find_maximum_range(angle, velocity)
max_height = find_maximum_height(angle, velocity)

print("\nProjectile Motion Analysis:")
print(f"Launch angle: {angle}°")
print(f"Initial velocity: {velocity} m/s")
print(f"Maximum range: {max_range:.2f} m")
print(f"Maximum height: {max_height:.2f} m")

# Generate trajectory data points
print("\nTrajectory points (x, y):")
for i in range(0, 11):
    x = i * max_range / 10
    y = projectile_motion(angle, velocity, x)
    print(f"({x:.2f}, {y:.2f})")