Get multiplier and divisors for an input frequency for PLL using Continued Fractions

Program.cs

using System;
using System.Collections.Generic;
using System.Numerics;
using Fractions;

static Fraction FractionCoefsToRational(IReadOnlyList<BigInteger> approx, int index = 0, int cutoff = -1)
{
    if (cutoff > approx.Count || cutoff < 0) cutoff = approx.Count;
    
    var ret = Fraction.Zero;
    if (index >= cutoff)
        return ret;

    ret += approx[index++];
    
    if (index < cutoff)
    {
        ret += 1 / FractionCoefsToRational(approx, index, cutoff);            
    }

    return ret.Reduce();
}

static BigInteger[] GetFractionCoefs(Fraction f)
{
    var ret = new Queue<BigInteger>();
    
    while (f.Numerator != BigInteger.Zero)
    {
        var r = f.Numerator % f.Denominator;
        ret.Enqueue((f.Numerator - r) / f.Denominator);

        f = new Fraction(f.Denominator, r);
    }

    return ret.ToArray();
}

static void GetFactors(Fraction freq, out BigInteger m, out BigInteger n, out BigInteger k, out int tryCount, int spreadDirection)
{
    spreadDirection = Math.Sign(spreadDirection);
    
    tryCount = 1;
    
    m = 1;
    n = 1; 
    k = 1;

    Fraction referenceFraction = default;
    var referenceInitialized = false;

    while (tryCount <= 50)
    {
        var coefs = GetFractionCoefs(freq);

        var newFraction = FractionCoefsToRational(coefs);

        if (!referenceInitialized)
        {
            referenceFraction = newFraction;
            referenceInitialized = true;
        }

        var cutoff = coefs.Length;
        
        while (newFraction.Numerator > 512)
        {
            newFraction = FractionCoefsToRational(coefs, 0, --cutoff);
        }

        n = BigInteger.Zero;
        
        const int newNLimit = 512;

        for (var i = 2; i < newNLimit; i++)
        {
            var newN = new BigInteger(i);

            if (newFraction.Denominator % newN != BigInteger.Zero ||
                newFraction.Denominator / newN >= newNLimit) continue;
            
            n = newN;
            break;
        }
        
        if (n == BigInteger.Zero)
        {
            var spreadCount = spreadDirection * tryCount;
            
            if (tryCount > 10) spreadCount = spreadDirection * 10 * (tryCount - 10);
            if (tryCount > 20) spreadCount = spreadDirection * 100 * (tryCount - 20);
            if (tryCount > 30) spreadCount = spreadDirection * 1000 * (tryCount - 30);
            if (tryCount > 40) spreadCount = spreadDirection * 10000 * (tryCount - 40);

            freq = referenceFraction + new Fraction(spreadCount, 1000000);
            
            tryCount++;
            continue;
        }
        
        m = newFraction.Numerator;
        k = newFraction.Denominator / n;
        
        break;
    }
}

// Preparation

// Fixed seed for repeatability
var rand = new Random(25543);

var listOfFrequencies = new List<Fraction>();

var randCount = 500_000;

for (var i = 0; i < randCount; i++)
{
    // Generate a random frequency from 3 MHz to 15 MHz
    var freq = new Fraction(rand.Next(3_000_000, 15_000_000), 1_000_000);
    freq /= 50; // 50 is the input frequency
    
    listOfFrequencies.Add(freq);
}

Console.WriteLine("--------");

var outlierCount = 0;
var maxDrift = Fraction.Zero;

foreach (var freq in listOfFrequencies)
{
    GetFactors(freq, out var m, out var n, out var k, out var tryCount, -1);

    var originalFreq = freq * 50;
    var newFrequency = new Fraction(m, n * k) * 50;

    var drift = ((newFrequency - originalFreq) / originalFreq * 100 ).Abs();

    if (drift > maxDrift) maxDrift = drift;

    if (drift > new Fraction(3, 100))
    {
        outlierCount++;
    }
    
    // Console.WriteLine(
    //     "original frequency: {0:0.000000} MHz, resulting frequency: {1:0.000000} MHz, drift: {2:0.00} %, direction: {3}, tries: {4}", 
    //     originalFreq.ToDecimal(), 
    //     newFrequency.ToDecimal(), 
    //     drift.ToDecimal(),
    //     1,
    //     tryCount
    // );
    // Console.WriteLine("frac: {3}, m: {0}, n: {1}, k: {2},", m, n, k, new Fraction(m, n * k));
    //
    // Console.WriteLine("--------");   
}

Console.WriteLine("Outliers: {0}, Percentage of outliers: {1:0.00} %", outlierCount, (double) outlierCount / listOfFrequencies.Count * 100);
Console.WriteLine("Maximum drift: {0:0.00} %", maxDrift.ToDecimal());