bkm.c

/* SPDX-License-Identifier: MIT */
/*
 * bkm.c
 *
 * Initial implementations of lg2 and pow2 using long doubles.
 * Written to test the implemented algorithms in general,
 * before tailoring them for fxp's
 *
 * Also generates the tables of pre-calculated values for
 * long and int fxp implementations of BKM-based lg2 and pow2
 *
 * By Raul Saavedra, Bonn Germany
 */

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include "fxp_extern.h"

const long double FXP_ZERO_LD = 1.0E-124L;

#define DASHES "=========================================\n"

static const long double A_2[] = {
    1.0000000000000000000000000000000000000000000000000000000000000000000000000000L,
    0.5849625007211561814537389439478165087598144076924810604557526545410982276485L,
    0.3219280948873623478703194294893901758648313930245806120547563958159347765589L,
    0.1699250014423123629074778878956330175196288153849621209115053090821964552970L,
    0.0874628412503394082540660108104043540112672823448206881266090643866965081686L,
    0.0443941193584534376531019906736094674630459333742491317685543002674288465967L,
    0.0223678130284545082671320837460849094932677948156179815932199216587899627785L,
    0.0112272554232541203378805844158839407281095943600297940811823651462712311786L,
    0.0056245491938781069198591026740666017211096815383520359072957784732489771013L,
    0.0028150156070540381547362547502839489729507927389771959487826944878598909400L,
    0.0014081943928083889066101665016890524233311715793462235597709051792834906001L,
    0.0007042690112466432585379340422201964456668872087249334581924550139514213168L,
    0.0003521774803010272377989609925281744988670304302127133979341729842842377649L,
    0.0001760994864425060348637509459678580940163670081839283659942864068257522373L,
    0.0000880524301221769086378699983597183301490534085738474534831071719854721939L,
    0.0000440268868273167176441087067175806394819146645511899503059774914593663365L,
    0.0000220136113603404964890728830697555571275493801909791504158295359319433723L,
    0.0000110068476674814423006223021573490183469930819844945565597452748333526464L,
    0.0000055034343306486037230640321058826431606183125807276574241540303833251704L,
    0.0000027517197895612831123023958331509538486493412831626219340570294203116559L,
    0.0000013758605508411382010566802834037147561973553922354232704569052932922954L,
    0.0000006879304394358496786728937442939160483304056131990916985043387874690617L,
    0.0000003439652607217645360118314743718005315334062644619363447395987584138324L,
    0.0000001719826406118446361936972479533123619972434705828085978955697643547921L,
    0.0000000859913228686632156462565208266682841603921494181830811515318381744650L,
    0.0000000429956620750168703982940244684787907148132725669106053076409624949917L,
    0.0000000214978311976797556164155504126645192380395989504741781512309853438587L,
    0.0000000107489156388827085092095702361647949603617203979413516082280717515504L,
    0.0000000053744578294520620044408178949217773318785601260677517784797554422804L,
    0.0000000026872289172287079490026152352638891824761667284401180026908031182361L,
    0.0000000013436144592400232123622589569799954658536700992739887706412976115422L,
    0.0000000006718072297764289157920422846078078155859484240808550018085324187007L,
    0.0000000003359036149273187853169587152657145221968468364663464125722491530858L,
    0.0000000001679518074734354745159899223037458278711244127245990591908996412262L,
    0.0000000000839759037391617577226571237484864917411614198675604731728132152582L,
    0.0000000000419879518701918839775296677020135040214077417929807824842667285938L,
    0.0000000000209939759352486932678195559552767641474249812845414125580747434389L,
    0.0000000000104969879676625344536740142096218372850561859495065136990936290929L,
    0.0000000000052484939838408141817781356260462777942148580518406975851213868092L,
    0.0000000000026242469919227938296243586262369156865545638305682553644113887909L,
    0.0000000000013121234959619935994960031017850191710121890821178731821983105443L,
    0.0000000000006560617479811459709189576337295395590603644549624717910616347038L,
    0.0000000000003280308739906102782522178545328259781415615142931952662153623493L,
    0.0000000000001640154369953144623242936888032768768777422997704541618141646683L,
    0.0000000000000820077184976595619616930350508356401599552034612281802599177300L,
    0.0000000000000410038592488303636807330652208397742314215159774270270147020117L,
    0.0000000000000205019296244153275153381695384157073687186580546938331088730952L,
    0.0000000000000102509648122077001764119940017243502120046885379813510430378661L,
    0.0000000000000051254824061038591928917243090559919209628584150482483994782302L,
    0.0000000000000025627412030519318726172939815845367496027046030028595094737777L,
    0.0000000000000012813706015259665053515049475574143952543145124550608158430592L,
    0.0000000000000006406853007629833949364669629701200556369782295210193569318434L,
    0.0000000000000003203426503814917330334121037829290364330169106716787999052925L,
    0.0000000000000001601713251907458754080007074659337446341494733882570243497196L,
    0.0000000000000000800856625953729399268240176265844257044861248416330071223615L,
    0.0000000000000000400428312976864705191179247866966320469710511619971334577509L,
    0.0000000000000000200214156488432353984854413866994246781519154793320684126179L,
    0.0000000000000000100107078244216177339743404416874899847406043033792202127070L,
    0.0000000000000000050053539122108088756700751579281894640362199287591340285355L,
    0.0000000000000000025026769561054044400057638132352058574658089256646014899499L,
    0.0000000000000000012513384780527022205455634651853807110362316427807660551208L,
    0.0000000000000000006256692390263511104084521222346348012116229213309001913762L,
    0.0000000000000000003128346195131755552381436585278035120438976487697544916191L,
    0.0000000000000000001564173097565877776275512286165232838833090480508502328437L,
    0.0000000000000000000782086548782938888158954641464170239072244145219054734086L,
    0.0000000000000000000391043274391469444084776945327473574450334092075712154016L,
    0.0000000000000000000195521637195734722043713378812583900953755962557525252782L,
    0.0000000000000000000097760818597867361022187915943503728909029699365320287407L,
    0.0000000000000000000048880409298933680511176764606054809062553340323879609794L,
    0.0000000000000000000024440204649466840255609083961603140683286362962192177597L,
    0.0000000000000000000012220102324733420127809717395445504379645613448652614939L,
    0.0000000000000000000006110051162366710063906152551383735699323415812152114058L,
    0.0000000000000000000003055025581183355031953399739107113727036860315024588989L,
    0.0000000000000000000001527512790591677515976780735407368332862218276873443537L,
    0.0000000000000000000000763756395295838757988410584167137033767056170417508383L,
    0.0000000000000000000000381878197647919378994210346199431733717514843471513618L,
    0.0000000000000000000000190939098823959689497106436628681671067254111334889005L,
    0.0000000000000000000000095469549411979844748553534196582286585751228071408728L,
    0.0000000000000000000000047734774705989922374276846068851506055906657137209047L,
    0.0000000000000000000000023867387352994961187138442777065843718711089344045782L,
    0.0000000000000000000000011933693676497480593569226324192944532044984865894525L,
    0.0000000000000000000000005966846838248740296784614396011477934194852481410926L,
    0.0000000000000000000000002983423419124370148392307506484490384140516252814304L,
    0.0000000000000000000000001491711709562185074196153830361933046331030629430117L,
    0.0000000000000000000000000745855854781092537098076934460888486730708440475045L,
    0.0000000000000000000000000372927927390546268549038472050424734256652501673274L,
    0.0000000000000000000000000186463963695273134274519237230207489851150821191330L,
    0.0000000000000000000000000093231981847636567137259618916352525606281553180093L,
    0.0000000000000000000000000046615990923818283568629809533488457973317312233323L,
    0.0000000000000000000000000023307995461909141784314904785572277779202790023236L,
    0.0000000000000000000000000011653997730954570892157452397493151087737428485431L,
    0.0000000000000000000000000005826998865477285446078726199923328593402722606924L,
    0.0000000000000000000000000002913499432738642723039363100255852559084863397344L,
    0.0000000000000000000000000001456749716369321361519681550201473345138307215067L,
    0.0000000000000000000000000000728374858184660680759840775119123438968122488047L,
    0.0000000000000000000000000000364187429092330340379920387564158411083803465567L,
    0.0000000000000000000000000000182093714546165170189960193783228378441837282509L,
    0.0000000000000000000000000000091046857273082585094980096891901482445902524441L,
    0.0000000000000000000000000000045523428636541292547490048446022564529197237262L,
    0.0000000000000000000000000000022761714318270646273745024223029238091160103901L
};

/*
 * Calculate the log2 of a long double number
 * (Implementing and testing the logarithm algorithm first
 * with floating points, before implementing it for fxp's)
 *
 * Based on the general algorithm to calculate binary logarithms
 * explained by Clay Turner in IEEE Signal Processing Magazine,
 * Sep/2010. His short article has D. E. Knuth's "The Art of
 * Computer Programming Vol 2: Seminumerical Algorithms",
 * 2nd ed., 1981 (pages 441 - 446) as the only reference.
 * Additional useful resources: CORDIC and/or BKM algorithms
 */
long double my_log2(long double x)
{
        if (x < 0.0L) return FXP_UNDEF_LD;
        if (x <= FXP_ZERO_LD) return FXP_NINF_LD;
        long double z = x;
        int c = 0; // characteristic
        while (z >= 2.0L) {
                c++;
                z /= 2;
        }
        while (z < 1.0L) {
                c--;
                z *= 2;
        }
        // Here we have already calculated the log characteristic c, and
        // we have z satisfying: 1 <= z < 2, so we can use it to calculate
        // the mantissa
        //printf("log c: %d\n", c);
        long double m = 0.0L; // mantissa
        long double b = 0.5L;
        int nb = 64; // desired number of mantissa bits to process
        while (nb > 0) {
                z = z * z;
                if (z >= 2.0L) {
                        z /= 2;
                        m += b;
                }
                b = b/2;
                nb--;
        }
        // Here we have already calculated the mantissa down to
        // nb bits of precision
        //printf("log m: %Lf\n", m);
        long double full_log = ((long double) c) + m;
        //printf("log(%Lf) = c + m = %Lf\n", x, full_log);
        return full_log;
}

/*
 * Implementation of the BKM L-Mode algorithm for log2 calculation
 * using long doubles:
 * https://en.wikipedia.org/wiki/BKM_algorithm
 */
long double my_log2_bkm(long double xin, int nbits)
{
        if (xin < 0.0L) return FXP_UNDEF_LD;
        if (xin <= FXP_ZERO_LD) return FXP_NINF_LD;
        long double argument = xin;
        int c = 0; // characteristic
        while (argument >= 2.0L) {
                c++;
                argument /= 2;
        }
        while (argument < 1.0L) {
                c--;
                argument *= 2;
        }
        // Here we have already calculated the log characteristic c, and
        // we have z satisfying: 1 <= z < 2
        long double x = 1.0L;
        long double y = 0.0L;
        long double s = 0.5L;
        // Notice ss starting with 0.5. Skipping the 1.0 because we know
        // zz < 2, so that first test if (zz < z) when ss == 1.0 would be
        // false for sure
        long double z;
        for (int k = 1; k < nbits; k++) {
                // Notice that z is for sure in [1, 2), and yy
                // will remain in the range of log(z), and that is:
                // [0, 1), but zz can in fact get bigger than 2
                // repeatedly in these iterations, even if not by much.
                z = x + x * s;
                printf("\tlog2 bkm iteration %d, z:%0.5LE\n", k, z);
                if (z <= argument) {
                        x = z;
                        y += A_2[k];
                        printf("\t\tUpdating y:%0.5LE\n", y);
                }
                s *= 0.5;
        }
        long double full_log2 = ((long double) c) + y;
        printf("log2(%.5LE) = c (%d)+ m (%.5LE) = %.5LE\n", xin, c, y, full_log2);
        return full_log2;
}

/*
 * Implementation of the BKM E-Mode algorithm for 2^n calculation
 * using long doubles
 */
long double my_pow2_bkm(long double n, int nbits)
{
        // Notice that 2^n == 2^(whole(n) + frac(n)) ==
        // 2^whole(n) * 2^frac(n) = pow2(w) * pow2(f)
        long long w = truncl(n);
        long double frac = n - w;
        // Calculate pow2(w), and prepare argument for BKM
        long double pow2w;
        long double argument;
        if (n >= 0) {
                pow2w = (long double) (1l << w);
                argument = frac;    // argument >= 0
        } else {
                pow2w = 1.0 / ((long double) (1l << (-w + 1)));
                argument = 1 + frac; // Notice argument in (0, 1]
        }
        printf("n:%.18Lf,  w:%lld,  pow2w:%18Lf,  argument:%.18Lf\n", \
                    n, w, pow2w, argument);
        // The BKM algorithm in E-Mode (for exponential)
        // calculates pow2(argument), argument in [0, 1) here
        long double x = 1.0L, y = 0.0L, s = 1.0L;
        for (int k = 0; k < nbits; k++) {
                long double const  z = y + A_2[k];
                printf("\tpow2 bkm iteration %d, z:%.18Lf, x:%.18Lf\n", \
                            k, z, x);
                if (z <= argument) {
                        y = z;
                        x = x + x*s;
                        printf("\t\tUpdating y:%.18Lf x:%.18Lf\n", y, x);
                }
                s *= 0.5L;
        }
        // Here x == 2^argument == pow2(f)
        // Calculate the full 2^n = pow2w * pow2f
        long double full_pow2 = pow2w * x;
        printf("pow2(%.18Lf) = %.18Lf\n", n, full_pow2);
        return full_pow2;
}

int main(void)
{
        printf("\n%sGenerate table of log values for BKM.c\n%s", DASHES, DASHES);
        const int nfbits = 63;
        printf("\nValues for the BKM array, long hex format, %d frac bits:\n", nfbits);
        int n = sizeof(A_2) / sizeof(A_2[0]);
        int j = 1;
        unsigned long mult = 1llu << nfbits;
        printf("\nstatic const unsigned long FXP_BKM_LOGS_L[] = {\n");
        for (int i=0; i < n; i++) {
                long double numld = A_2[i];
                unsigned long num = truncl( numld * mult );
                //printf("%d:\t0x%llX \t(%LE)\n", i, num, numld);
                printf("0x%lX", num);
                if (num == 0L) break;
                printf(", ");
                if (((i + 1) % 4) == 0) printf("\n");
        }
        printf("\n};\n");

        unsigned long msbits, lsbits;
        printf("\nSplit values for the BKM arrays using only ints:\n");
        printf("\nFirst array:\n");
        printf("\nstatic const ulongy FXP_BKM_LOGS_NEW[] = {\n");
        for (int i=0; i < n; i++) {
                long double numld = A_2[i];
                unsigned long num = truncl( numld * mult );
                msbits = (num & 0xFFFFFFFF00000000llu) >> 32;
                printf("{ 0x%08lX, ", msbits);
                //if (((i + 1) % 4) == 0) printf("\n");
                lsbits = (num & 0x00000000FFFFFFFFllu);
                printf("0x%08lX }", lsbits);
                if (i<n-1) printf(", ");
		printf("\n");
                if ((msbits == 0) && (lsbits == 0)) break;
        }
        printf("};\n");

        //long double x1 = my_log2_bkm(1.9999999999, 31);
        //long double x2 = my_log2_bkm(0.9999999999, 31);

        printf("\n");
        long double p = 1.0L - 1.0L/(1u << 30);
        long double x3 = my_pow2_bkm(p, 60);
        printf("Expected 2^(1 - 1/2^30):     %.18Lf\n\n", powl(2.0L, p));
        //long double x4 = my_pow2_bkm(1.5, 31);
        //long double x5 = my_pow2_bkm(2.0, 31);
        //long double x6 = my_pow2_bkm(2.5, 31);


/*
        printf("\nPrecision of long double pow2() calculations using BKM:\n");
        long double invln2 = 1 / logl(2);
        for (int bkmd = 29; bkmd <= 31; bkmd++) {
                printf("\nChecking BKM depth %d\n", bkmd);
                // Making sure it works as expected even beyond 32 bits
                long double nums[] = {0.0, 0.25, 0.5, 1.0, 2.0, 3.1, 31.1, 40.1, \
                                    -0.25, -0.5, -1.0, -2.0, -3.1, -31.1, -40.1};
                int n = sizeof(nums) / sizeof(nums[0]);
                for (int i=0; i < n; i++) {
                        long double tgt = powl(2.0, nums[i]);
                        long double mypow2 = my_pow2_bkm(nums[i], bkmd);
                        printf("\n\tref: pow(2, %.3LE) = %.20LE\n",
                                nums[i], tgt);
                        printf("\tmy_pow2(%.3LE)     = %.20LE (- ref: %1.2LE)\n",
                                nums[i], mypow2, (mypow2 - tgt));
                }
        }
*/

}