Blog Entry © Tuesday, October 21, 2025, by James Pate Williams, Jr., Solving Low Density Subset Sum Problems Using the LLL-Lattice Reduction Algorithm

// Algorithm found in the "Handbook of
// Applied Cryptography" (c) 1997 by
// Alfred J. Menezes, Paul C van Oorschot,
// and Scott Vanstone 3.105 Algorithm
// Chapter 3 pages 120 - 121

#pragma once
class LLL_Lattice
{

private:

	static double Scalar(
		int n,
		std::vector<double> u,
		std::vector<double> v);
	static void RED(
		int k, int l, int n,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& h,
		std::vector<std::vector<double>>& mu);
	static void SWAP(
		int k, int k1, int kmax, int n,
		double m, std::vector<double>& B,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& bs,
		std::vector<std::vector<double>>& h,
		std::vector<std::vector<double>>& mu);

public:

	static bool LLL(
		int n,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& h);
};

#include "pch.h"
#include "LLL_Lattice.h"

double LLL_Lattice::Scalar(
    int n,
    std::vector<double> u,
    std::vector<double> v)
{
    // Calculate the scalar product of two vectors [1..n]
    double sum = 0.0;

    for (int i = 1; i <= n; i++) sum += u[i] * v[i];
    return sum;
}

void LLL_Lattice::RED(
    int k, int l, int n,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& h,
    std::vector<std::vector<double>>& mu)
{
    int i, q = (int)(0.5 + mu[k][l]);

    if (fabs(mu[k][l]) > 0.5)
    {
        for (i = 1; i <= n; i++)
        {
            b[k][i] -= q * b[l][i];
            h[i][k] -= q * h[i][l];
        }

        mu[k][l] -= q;

        for (i = 1; i <= l - 1; i++)
        {
            mu[k][i] -= q * mu[l][i];
        }
    }
}

void LLL_Lattice::SWAP(
    int k, int k1, int kmax, int n,
    double m, std::vector<double>& B,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& bs,
    std::vector<std::vector<double>>& h,
    std::vector<std::vector<double>>& mu)
{
    double C, t;
    int i, j;
    std::vector<double> c(n + 1);

    for (j = 1; j <= n; j++)
    {
        t = b[k][j];
        b[k][j] = b[k1][j];
        b[k1][j] = t;
        t = h[j][k];
        h[j][k] = h[j][k1];
        h[j][k1] = t;
    }

    if (k > 2)
    {
        for (j = 1; j <= k - 2; j++)
        {
            t = mu[k][j];
            mu[k][j] = mu[k1][j];
            mu[k1][j] = t;
        }
    }

    C = B[k] + m * m * B[k1];
    mu[k][k1] = m * B[k1] / C;

    for (i = 1; i <= n; i++)
    {
        c[i] = bs[k1][i];
    }

    for (i = 1; i <= n; i++)
    {
        bs[k1][i] = bs[k][i] + m * c[i];
    }
    
    for (i = 1; i <= n; i++)
    {
        bs[k][i] = -m * bs[k][i] + B[k] * c[i] / C;
    }

    B[k] *= B[k1] / C;
    B[k1] = C;

    for (i = k + 1; i <= kmax; i++)
    {
        t = mu[i][k];
        mu[i][k] = mu[i][k1] - m * t;
        mu[i][k1] = t + mu[k][k1] * mu[i][k];
    }
}

bool LLL_Lattice::LLL(
    int n,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& h)
{
    // Lattice reduction algorithm
 
    double m;
    std::vector<double> B(n + 1);
    std::vector<double> bv(n + 1);
    std::vector<double> bw(n + 1);
    std::vector<std::vector<double>> bs(n + 1,
        std::vector<double>(n + 1));
    std::vector<std::vector<double>> mu(n + 1,
        std::vector<double>(n + 1));
    int i, j, k, k1, kmax = 1, l;

    for (i = 1; i <= n; i++)
        bv[i] = bs[1][i] = b[1][i];

    B[1] = Scalar(n, bv, bv);

    for (i = 1; i <= n; i++)
    {
        for (j = 1; j <= n; j++)
        {
            h[i][j] = 0.0;
        }
        
        h[i][i] = 1.0;
    }

    for (k = 2; k <= n; k++)
    {
        if (k <= kmax)
            goto Label3;

        kmax = k;
        for (i = 1; i <= n; i++)
        {
            bs[k][i] = b[k][i];
        }

        for (j = 1; j <= k - 1; j++)
        {
            for (l = 1; l <= n; l++)
            {
                bv[l] = b[k][l];
                bw[l] = bs[j][l];
            }

            mu[k][j] = Scalar(n, bv, bw) / B[j];

            for (i = 1; i <= n; i++)
            {
                bs[k][i] -= mu[k][j] * bs[j][i];
            }
        }

        for (i = 1; i <= n; i++)
        {
            bv[i] = bs[k][i];
        }

        B[k] = Scalar(n, bv, bv);

        if (B[k] == 0.0)
            return false;

    Label3:

        k1 = k - 1;
        m = mu[k][k1];
        RED(k, k1, n, b, h, mu);

        if (B[k] < (0.75 - m * m) * B[k1])
        {
            SWAP(k, k1, kmax, n, m, B, b, bs, h, mu);
            k = max(2, k1);
            goto Label3;
        }

        for (l = k - 2; l >= 1; l--)
        {
            RED(k, l, n, b, h, mu);
        }
    }

    return true;
}

Blog Entry © Sunday, October 19, 2025, by James Pate Williams, Jr. LIPCalculator (Large Integer Package Calculator)

Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.

Blog Entry © Thursday, August 7, 2025, z-Component of Cylindrical Velocity Based on an Equation in Classical Electrodynamics Second Edition © 1975 by John David Jackson

Blog Entry © Monday August 4, 2025, by James Pate Williams, Jr. Bessel Functions of the First and Second Kind, Their Derivatives and Zeros

Blog Entry © Saturday, August 2, 2025, by James Pate Williams, Jr. Orthonormal Transverse Magnetic and Transverse Electric Fields Win32 Desktop C/C++ Application in the Release Configuration

Blog Entry © Friday, August 1, 2025, by James Pate Williams, Jr. Numerically Solving a Two-Dimensional Elliptic Partial Differential Equation (PDE) Boundary Value Problem (BVP) and a New MP3

Blog Entry © Thursday, July 31, 2025, Numerically Solving Four Second Order Linear Ordinary Differential Equation Boundary Value Problems by James Pate Williams, Jr.

Blog Entry (c) Wednesday, July 30, 2025, by James Pate Williams, Jr. Four Root Finding Algorithms and C Source

function exp(-x) - sin(0.5 * pi * x)
The interval found by bisection is:
a = +0.4435735341  b = +0.4435735341
bisection # iterations = 50
The root by regula falsi:
x = +0.4435735341  f(x) = -0.0000000000
regula falsi # iterations = 8
The root by the secant method:
x = +0.4435735341  f(x) = +0.0000000000
secant method # iterations = 9
The root by Newton's Method:
x = +0.4435735341  f(x) = -0.0000000000
Newton's Method # iterations = 6

D:\roots\x64\Release\roots.exe (process 28152) exited with code 0 (0x0).
Press any key to close this window . . .
#include <math.h>
#include <stdio.h>

typedef double real;

static real f(real x)
{
	double pi2 = 2.0 * atan(1.0);
	return(exp(-x) - sin(pi2 * x));
}

static real g(real x)
{
	double pi2 = 2.0 * atan(1.0);
	return(-exp(-x) - pi2 * cos(pi2 * x));
}

static int bisection(
	real(*f)(real), real* a, real* b, real xtol, int* flag)
{
	real error, fa, fm, xm;
	int n = 0;

	fa = (*f)(*a);
	if (fa * (*f)(*b) > 0.0)
	{
		*flag = -1;
		return(n);
	}
	error = fabsl(*b - *a);
	while (error > xtol)
	{
		n++;
		error *= 0.5;
		if (error <= xtol)
		{
			*flag = 0;
			return(n);
		}
		xm = 0.5 * (*a + *b);
		if (xm + error == xm)
		{
			*flag = 1;
			return(n);
		}
		fm = (*f)(xm);
		if (fa * fm > 0.0)
		{
			*a = xm;
			fa = fm;
		}
		else
			*b = xm;
	}
	*flag = 2;
	return(n);
}

static int regula(
	real (*f)(real), real a, real b, real xtol,
	real ftol, int ntol, real* w, int* flag)
{
	int n = 0;
	real fa, fb, fw, signfa, prvsfw;

	fa = f(a);
	if (fa >= 0.0) signfa = +1.0; else signfa = -1.0;
	fb = f(b);
	if (signfa * fb >= 0.0)
	{
		*flag = -1;
		return n;
	}
	*w = a;
	fw = fa;
	for (n = 0; n <= ntol; n++)
	{
		if (fabs(a - b) <= xtol)
		{
			*flag = 0;
			return n;
		}
		if (fabs(fw) <= ftol)
		{
			*flag = 1;
			return n;
		}
		*w = (fa * b - fb * a) / (fa - fb);
		if (fw >= 0.0) prvsfw = +1.0; else prvsfw = -1.0;
		fw = f(*w);
		if (signfa * fw > 0.0)
		{
			a = *w;
			fa = fw;
			if (fw * prvsfw > 0.0) fb = 0.5 * fb;
		}
		else
		{
			b = *w;
			fb = fw;
			if (fw * prvsfw > 0.0) fa = 0.5 * fa;
		}
	}
	*flag = 2;
	return n;
}

static int secantMethod(
	real(*f)(real), real xtol, real ftol,
	int ntol, real xm1, real x0, real* w)
{
	real fm1 = f(xm1), f0 = f(x0);
	real df = fabs(fm1 - f0), f1 = 0.0;
	real dx = fabs(xm1 - x0), x1 = 0.0;
	int n = 0;

	while (n < ntol && df > ftol && dx > xtol)
	{
		x1 = (f0 * xm1 - fm1 * x0) / (f0 - fm1);
		f1 = f(x1);
		df = fabs(f(x1) - f0);
		dx = fabs(x1 - x0);
		fm1 = f0;
		f0 = f1;
		xm1 = x0;
		x0 = x1;
		n++;
	}

	*w = x1;
	return n;
}

static int NewtonsMethod(
	real(*f)(real), real(*g)(real),
	real xtol, real ftol, int ntol,
	real* w)
{
	// f is the function
	// g is the function's derivative
	// xtol is root's tolerance
	// ftol is the function's tolerance
	// ntol is the maximum # of iterations
	real f1 = 0.0, g1 = 0.0, x0 = *w, x1 = 0.0;
	real f0 = f(x0), g0 = g(x0);
	real deltaX = DBL_MAX, deltaF = DBL_MAX;
	int n = 0;

	while (n < ntol && deltaX > xtol && deltaF > ftol)
	{
		x1 = x0 - f0 / g0;
		f1 = f(x1);
		g1 = g(x1);
		deltaX = fabs(x1 - x0);
		deltaF = fabs(f1 - f0);
		f0 = f1;
		g0 = g1;
		x0 = x1;
		n++;
	}

	*w = x1;
	return n;
}

int main(void)
{
	int flag = 0, ntol = 0;
	real a0 = 0, b0 = 0, ftol = 0, w1 = 0, w2 = 1.0;
	real a1 = 0, b1 = 0, w3 = 0, xtol = 0;

	a0 = 0.0;
	b0 = 1.0;
	ntol = 128;
	ftol = 1.0e-15;
	xtol = 1.0e-15;
	int its1 = bisection(f, &a0, &b0, xtol, &flag);
	a1 = 0.0;
	b1 = 1.0;
	int its2 = regula(f, a1, b1, xtol, ftol, ntol, &w1, &flag);
	int its3 = secantMethod(f, ftol, xtol, ntol, 0.0, 1.0, &w2);
	int its4 = NewtonsMethod(f, g, xtol, ftol, ntol, &w3);
	printf("function exp(-x) - sin(0.5 * pi * x)\n");
	printf("The interval found by bisection is:\n");
	printf("a = %+13.10lf  b = %+13.10lf\n", a0, b0);
	printf("bisection # iterations = %ld\n", its1);
	printf("The root by regula falsi:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w1, f(w1));
	printf("regula falsi # iterations = %ld\n", its2);
	printf("The root by the secant method:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w2, f(w2));
	printf("secant method # iterations = %ld\n", its3);
	printf("The root by Newton's Method:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w2, f(w3));
	printf("Newton's Method # iterations = %ld\n", its4);
	return(0);
}