# Crystals, Oscillators & ceramic resonators.  All freq in MHz
# TODO: find way to specify tolerance: currently just using stated frequency
#       and can count number of zeros...
# TODO: Modify to allow M or K to be decimal point to specify other freq ranges

# Oscillators (cans)
OSC     20.0000
OSC     28.3220
OSC     42.0000
OSC     36.0000
OSC     36.00000
OSC     50.000
OSC     60.0

# Bare crystals
XTAL    .032768
XTAL    .0384
XTAL    1.8432
XTAL    3.58    (tiny)
XTAL    4.00    (tiny)
XTAL    4.5
XTAL    7.3728
XTAL    6.000
XTAL    6.5536
XTAL    8.0000
XTAL    14.31818
XTAL    20.000
XTAL    27.115
XTAL    27.125
XTAL    30.000
XTAL    31
# UNKNOWNS:
# BOMAF C-6 1W14T | MEW 31 PU47931-2
# A set marked 301-502 on side, top:  (14,4,3,2,-1,-5,-8,-9,-11,-12,-13)

# Ceramic resonators
RES     4.000
RES     8.000
RES     7.37        # ? "737 Cm 219"
RES     10.7
RES     16.93

##</code><code>##

actualFreq (error) = startFreq / divisor "div =" factors

##</code><code>##

$ ./find_xtal.pl 175k k5


    175,000 (      0) =  42,000,000 /  240  div = 2^4 * 3 * 5     
                      =  27,125,000 /  155  div = 5 * 31          
    174,935 (     65) =  27,115,000 /  155  div = 5 * 31          
    174,927 (     73) =  60,000,000 /  343  div = 7^3             
    175,141 (    141) =  31,000,000 /  177  div = 3 * 59          
    174,827 (    173) =  28,322,000 /  162  div = 2 * 3^4         
    174,825 (    175) =  50,000,000 /  286  div = 2 * 11 * 13     
    174,757 (    243) =  36,000,000 /  206  div = 2 * 103         
    174,611 (    389) =  14,318,180 /   82  div = 2 * 41          
    175,409 (    409) =  10,700,000 /   61  div = 61              
    175,438 (    438) =  20,000,000 /  114  div = 2 * 3 * 19      
                      =  30,000,000 /  171  div = 3^2 * 19        
                      =  60,000,000 /  342  div = 2 * 3^2 * 19    
                      =  50,000,000 /  285  div = 3 * 5 * 19      
    174,536 (    464) =  16,930,000 /   97  div = 97              
    175,476 (    476) =   7,370,000 /   42  div = 2 * 3 * 7       


##</code><code>##

#!/usr/bin/perl
my $usage = <<EOHDR;

    find_xtal.pl <DesiredFreq> [<Tolerance>]

    Find the crystal (and integral divisor(s)) to come within <Tolerance> Hz of 
    <DesiredFreq>.  Specify frequency as #.# or #m# for MHz, or #k# for kHz.
    If you don't specify <Tolerance>, 5% is assumed.

EOHDR
#
# 20130216 MCMason: Added tolerance, sort by error (asc) and funkiness within error.
#                   Find *all* divs within TOL range
# 20120429 MCMason: original version
#
#   TODO:   Trim report: If there are multiple ways to get same frequency, keep only
#           the top N of best one(s) (based on funkiness)
#   TODO:   Improve funkiness calculation: e.g. try to get something proportional
#           to # dividers & difficulty.  Perhaps something like $e * (log($f)/log(2))
#
use strict;
use warnings;
use autodie;

my $dbg_funkiness=0;

my $Freq = shift or die $usage;
$Freq = txt_2_MHz($Freq);
my $Tol = shift;
if (defined $Tol) {
    $Tol = txt_2_MHz($Tol);
}
else {
    $Tol = $Freq * .05;
}

my @crystals = read_crystals();
my @results;

sub compute {
    my ($Fx, $d) = @_;
    return undef unless $d;

    my $f = int($Fx / $d);
    my $err = abs($Freq - $f);
    return undef if $err > $Tol;

    my $ar = [ factorize($d) ];
    my $funk = factor_funkiness(@$ar);
    return [ $Fx, $d, $f, $err, $ar, $funk ];
}

for my $Fx (@crystals) {
    # 'perfect' divisor
    my $div = $Fx / $Freq;

    # Surrounding integral divisors
    my $d = int $div;
    while (my $ar = compute($Fx, $d)) {
        last if ! defined $ar;
        my ($xFx, $xd, $xf, $xerr, $xar, $xfunk) = @$ar;
        print "Fx:$xFx, d:$xd, f:$xf, err:$xerr, funk:$xfunk  (div:$div)\n" if $dbg_funkiness;
        push @results, $ar;
        --$d;
    }
    $d = int $div+1;
    while (my $ar = compute($Fx, $d)) {
        last if ! defined $ar;
        my ($xFx, $xd, $xf, $xerr, $xar, $xfunk) = @$ar;
        print "Fx:$xFx, d:$xd, f:$xf, err:$xerr, funk:$xfunk  (div:$div)\n" if $dbg_funkiness;
        push @results, $ar;
        ++$d;
    }
}

print "\n\n";
@results = sort { $$a[3] <=> $$b[3] or $$a[5] <=> $$b[5] } @results;

my $prev_freq = -1;
for my $ar (@results) {

    my $funk = $dbg_funkiness ? " \tfunkiness=$$ar[5]" : "";

    if ($prev_freq ne $$ar[2]) { 
        printf "%11.11s (%7s) = %11s / %4s  div = %-16s$funk\n",
                commify($$ar[2]), commify($$ar[3]), commify($$ar[0]), commify($$ar[1]),
                join(" * ", @{$$ar[4]});
        $prev_freq = $$ar[2];
    }
    else {
        printf "                      = %11s / %4s  div = %-16s$funk\n",
                commify($$ar[0]), commify($$ar[1]),
                join(" * ", factorize($$ar[1]));
    }
}

sub txt_2_MHz {
    my $t = shift;
    my $orig = $t;

    # Multiplier (default is MHz for /\d+\.\d*/)
    my $mult = 1000000;

    if ($t =~ /m/i) {
        #  1M8432 => 1,843,000
        $t =~ s/m/./;
        $mult = 1000000;
    }
    elsif ($t =~ /k/i) {
        # 32k768  =>    32,768
        $t =~ s/k/./;
        $mult = 1000;
    }

    if ($t =~ /\d+(\.\d*)?|\.\d+/) {
        return $t * $mult;
    }

    die "txt_2_MHz: Unexpected input '$orig' => '$t', can't determine frequency.";
}

sub read_crystals {
    open my $FH, '<', 'crystals.txt';
    my %freqs;
    while (<$FH>) {
        next if /^\s*#/;
        next if /^\s*$/;
        if (/([.0-9]+)/) {
            $freqs{$1*1_000_000}=0;
        }
    }
    return sort { $a <=> $b } keys %freqs;
}

sub commify {
    my $s = shift;
    $s =~ s/(\d)(\d\d\d)$/$1,$2/;
    $s =~ s/(\d)(\d\d\d,)/$1,$2/g;
    return $s;
}

sub factor_funkiness {
    my @factors = @_;
    my $funkiness = 0;
    print "funkiness(",join(", ", @factors),")\n" if $dbg_funkiness;
    for my $t (@factors) {
        my ($f, $e) = split /\^/, $t;
        $e=1 if ! defined $e;

        if ($f == 2) {
            $funkiness += .05*$e; # no change
        }
        elsif ($f == 3) {
            $funkiness += 0.5 * $e;
        }
        elsif ($f == 5) {
            $funkiness += 0.6 * $e;
        }
        elsif ($f < 32) {
            $funkiness += 1 * $e;
        }
        else {
            $funkiness += 2 * $e;
        }
        print "\t$t => $funkiness\n" if $dbg_funkiness;
    }

    print "\tfinal == $funkiness\n" if $dbg_funkiness;
    return $funkiness;
}

sub factorize {
    my $num = shift;
    my @factors = ();

    my $factor=2;
    my $exp=0;
    while ($num%$factor == 0) {
        ++$exp;
        $num /= $factor;
    }
    if ($exp > 1) { push @factors, "$factor^$exp" }
    elsif ($exp)  { push @factors, $factor }

    $factor=3;
    while ($factor*$factor <= $num) {
        $exp = 0;
        while ($num%$factor == 0) {
            ++$exp;
            $num /= $factor;
        }
        if ($exp > 1) { push @factors, "$factor^$exp" }
        elsif ($exp)  { push @factors, $factor }
        $factor += 2;
    }
    push @factors, $num if $num>1;

    return @factors;
}