text:
#pragma once

#include "../common.h"
#include "../concepts.h"
#include <boost/container/flat_set.hpp>
#include <boost/outcome/result.hpp>
#include <boost/outcome/try.hpp>
#include <fmt/format.h>
#include <gsl/gsl-lite.hpp>
#include <string>
#include <string_view>
#include <yaml-cpp/yaml.h>

namespace angonoka::validation {
using boost::container::flat_set;
namespace bo = boost::outcome_v2;
using result = bo::result<void, std::string>;
using namespace fmt::literals;

template <typename T>
concept Check
    = std::is_invocable_v<T, const YAML::Node&, gsl::czstring>;
template <typename T>
concept Attribute = std::is_convertible_v<
    decltype(std::declval<T>().name),
    std::string_view> && Check<T>;
template <typename T>
concept AttrOrStr = String<T> || Attribute<T>;

/**
    YAML scalar.

    Example:

    hello: "world"
           ^
           | Scalar

    scalar()

    Means that the value has to be a scalar and
    not a map or a sequence, etc.

    @return Check function
*/
constexpr Check auto scalar()
{
    return
        [](const YAML::Node& node, std::string_view scope) -> result {
            if (!node || node.IsNull())
                return R"("{}" can't be empty.)"_format(scope);
            if (!node.IsScalar())
                return R"("{}" has invalid type.)"_format(scope);
            return bo::success();
        };
}

namespace detail {
    /**
        Helper class for required and optional YAML fields.

        @var name   Parameter's name
        @var check  Function to apply to the field
    */
    template <Check T> struct functor {
        gsl::czstring name;
        T check;

        constexpr functor(gsl::czstring name, T check)
            : name{name}
            , check{check}
        {
        }
        explicit constexpr functor(gsl::czstring name)
            : functor{name, scalar()}
        {
        }
    };

    /**
        Join YAML path parts.

        @param a Root path
        @param b New path

        @return Concatenation of parts separated by a ".".
    */
    template <String T1, String T2> std::string join(T1&& a, T2&& b)
    {
        if (std::empty(a)) return std::forward<T2>(b);
        return "{}.{}"_format(
            std::forward<T1>(a),
            std::forward<T2>(b));
    }
} // namespace detail

/**
    Requred YAML field.

    Example:

    required("hello")

    Means that the field "hello" is required and
    it has to be a scalar value.

    @var name   Parameter's name
    @var check  Function to apply to the field
*/
template <Check T> struct required : detail::functor<T> {
    using detail::functor<T>::functor;
    result
    operator()(const YAML::Node& node, std::string_view scope) const
    {
        if (const auto n = node[this->name])
            return this->check(n, detail::join(scope, this->name));
        return R"("{}" is missing a "{}" attribute.)"_format(
            scope,
            this->name);
    }
};

template <Check T> required(gsl::czstring, T) -> required<T>;
required(gsl::czstring)->required<decltype(scalar())>;

namespace detail {
    /**
        Extract the map attribute's name.

        If an attribute is a string literal, pass the argument as is.

        @param attr Either an attribute or a string literal

        @return Attribute name
    */
    constexpr auto attr_name(gsl::czstring attr) { return attr; }
    constexpr auto attr_name(Attribute auto&& attr)
    {
        return attr.name;
    }

    /**
        Extract or construct an attribute check function.

        If an attrubte is a string literal, construct
        the required attrubte with the string literal as it's
        name.

        @param attr Either an attribute or a string literal

        @return Check function
    */
    constexpr auto attr_check(gsl::czstring attr)
    {
        return required(attr);
    }
    constexpr auto attr_check(Attribute auto&& attr) { return attr; }
} // namespace detail

/**
    Optional YAML field.

    Example:

    optional("hello")

    Means that the field "hello" is optional and if
    present, it has to be a scalar value.

    @var name  Parameter's name
    @var check Function to apply to the field
*/
template <Check T> struct optional : detail::functor<T> {
    using detail::functor<T>::functor;

    result
    operator()(const YAML::Node& node, std::string_view scope) const
    {
        if (const auto n = node[this->name])
            return this->check(n, detail::join(scope, this->name));
        return bo::success();
    }
};

template <Check T> optional(gsl::czstring, T) -> optional<T>;
optional(gsl::czstring)->optional<decltype(scalar())>;

/**
    YAML array.

    Validates each value of the array with the provided function.

    Example:

    sequence(scalar())

    Means that the value has to be a sequence (array) of
    scalar values.

    @param check Function to apply to each item

    @return Check function
*/
constexpr Check auto sequence(Check auto check)
{
    return [=](const YAML::Node& node,
               std::string_view scope) -> result {
        if (!node || !node.IsSequence()) {
            return R"("{}" is expected to be a sequence.)"_format(
                scope);
        }
        for (gsl::index i{0}; i < std::size(node); ++i) {
            BOOST_OUTCOME_TRY(
                check(node[i], "{}[{}]"_format(scope, i)));
        }
        return bo::success();
    };
}

constexpr Check auto sequence() { return sequence(scalar()); }

/**
    YAML map.

    Matches specified parameters exactly, no extra fields permitted.

    Example:

    attributes("first", optional("second"))

    Means that the value has to be a map with a required field
    "first", which has to be a scalar and an optional field "second"
    which also has to be a scalar.

    @param attrs Sequence of optional or required parameters

    @return Check function
*/
constexpr Check auto attributes(AttrOrStr auto... attrs)
{
    return [=](const YAML::Node& node,
               std::string_view scope = {}) -> result {
        if (!node || node.IsScalar() || node.IsSequence())
            return R"("{}" is expected to be a map.)"_format(scope);
        flat_set<std::string_view> unique_fields;
        for (auto&& n : node) {
            const auto& attr_name = n.first.Scalar();
            if (attr_name.empty())
                return R"(Empty attribute in "{}".)"_format(scope);
            if (!unique_fields.emplace(attr_name).second) {
                return R"(Duplicate attribute "{}" in "{}".)"_format(
                    attr_name,
                    scope);
            }
            if (!((attr_name == detail::attr_name(attrs)) || ...)) {
                return R"(Unexpected attribute "{}" in "{}".)"_format(
                    attr_name,
                    scope);
            }
        }
        result r = bo::success();
        ((r = detail::attr_check(attrs)(node, scope)) && ...);
        return r;
    };
}

/**
    YAML map.

    Validates each value of the map with the provided function.
    Used when the number of map fields may vary.

    Example:

    foo:
      bar1: 1
      bar2: 2
      bar3: 3

    values(scalar())

    @param check Function to apply to each value

    @return Check function
*/
constexpr Check auto values(Check auto check)
{
    return [=](const YAML::Node& node,
               std::string_view scope) -> result {
        if (!node || !node.IsMap())
            return R"("{}" is expected to be a map.)"_format(scope);
        for (auto&& n : node) {
            BOOST_OUTCOME_TRY(check(
                n.second,
                detail::join(scope, n.first.Scalar())));
        }
        return bo::success();
    };
}

/**
    Match at least one of the validators.

    Example:

    required("example", any_of(scalar(), attributes("foo", "bar")))

    Means the value has to either be a singular scalar value or
    a map with 2 fields "foo" and "bar.

    example: "hello"

    or

    example:
      foo: 1
      bar: 2

    @param check Functions to match

    @return Check function
*/
constexpr Check auto any_of(Check auto... checks)
{
    return [=](const YAML::Node& node,
               std::string_view scope) -> result {
        result r = bo::success();
        ((r = checks(node, scope)) || ...);
        return r;
    };
}
} // namespace angonoka::validation

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <omp.h>
#include <math.h>
#include "parmt_config.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "parmt_mtsearch.h"
#include "parmt_utils.h"
#include "compearth.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

#ifndef MAX
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
#endif

static int getNpMax(const int iobs, const struct parmtData_struct data,
                    const struct localMT_struct mtloc);
/*
static int setDataOnG(const int iobs, const int iloc, const int npmax,
                       const struct parmtData_struct data,
                       double *__restrict__ G);
*/
int parmt_computeL1StationMagnitude64f(const int ldm, const int nmt,
                                       const int npts, const int blockSize,
                                       const double *__restrict__ Gxx,
                                       const double *__restrict__ Gyy,
                                       const double *__restrict__ Gzz,
                                       const double *__restrict__ Gxy,
                                       const double *__restrict__ Gxz,
                                       const double *__restrict__ Gyz,
                                       const double *__restrict__ mts,
                                       const double *__restrict__ d,
                                       double *__restrict__ mags);

int parmt_locSearchL164f(const MPI_Comm locComm,
                         const int iobs, const int blockSize,
                         const int nlags, const bool lwantLags,
                         const bool lrescale,
                         struct localMT_struct mtloc,
                         struct parmtData_struct *data,
                         const double *__restrict__ CeInv,
                         double *__restrict__ phi,
                         double *__restrict__ var,
                         int *__restrict__ lags)
{
    const char *fcnm = "parmt_locSearchL164f\0";
    double *G, *d, *phiLoc, *phiWork, *varLoc, *varWork, xrescaleNum;
    int *lagLoc, *lagWork, ierr, iloc, jndx, jloc, k,
        myloc, npmax, nlocProcs, npts;
    const int master = 0;
    ierr = 0;
    myloc =-1;
    if (mtloc.myid == master)
    {
        MPI_Comm_rank(locComm, &myloc);
        MPI_Comm_size(locComm, &nlocProcs);
    }
    MPI_Bcast(&myloc,     1, MPI_INT, master, mtloc.comm);
    MPI_Bcast(&nlocProcs, 1, MPI_INT, master, mtloc.comm);
    if (iobs < 0 || iobs >= data->nobs)
    {
        printf("%s: Invalid observation number %d\n", fcnm, iobs);
        return -1;
    }
    if (mtloc.ldm <= 6)
    {
        printf("%s: Error leading dimension must be at least 6\n", fcnm);
        return -1;
    }
    if (mtloc.mts == NULL || phi == NULL)
    {
        if (mtloc.mts == NULL)
        {
            fprintf(stderr, "%s: Error mts is NULL\n", __func__);
        }
        if (phi == NULL){fprintf(stderr, "%s: Error phi is NULL\n", __func__);}
        return -1;
    }
    npmax = getNpMax(iobs, *data, mtloc);
    if (npmax < 1)
    {
        printf("%s: No data points\n", fcnm);
        return -1;
    }
    // Set space and copy data
    lagLoc = NULL;
    lagWork = NULL;
    varWork = NULL;
    phiWork = NULL;
    G = memory_calloc64f(6*npmax);
    d = memory_calloc64f(npmax);
    varLoc = array_set64f(npmax, DBL_MAX, &ierr);
    phiLoc = memory_calloc64f(mtloc.nmt);
    if (nlocProcs > 1)
    {
        if (mtloc.myid == master)
        {
            phiWork = memory_calloc64f(data->nlocs*mtloc.nmtAll);
            varWork = memory_calloc64f(npmax);
            if (lwantLags)
            {
                lagWork = memory_calloc32i(data->nlocs*mtloc.nmtAll);
            }
        }
        else
        {
            phiWork = memory_calloc64f(1);
            varWork = memory_calloc64f(1);
            lagWork = memory_calloc32i(1);
        }
    }
    else
    {
        phiWork = phi;
        varWork = var;
        lagWork = lags;
    }
    if (lwantLags)
    {
        lagLoc = memory_calloc32i(mtloc.nmt);
    }
    else
    {
        lagLoc = memory_calloc32i(1);
    }
    npts = data->data[iobs].npts;
    cblas_dcopy(npts, data->data[iobs].data, 1, d, 1);
    for (jloc=0; jloc<data->nlocs; jloc=jloc+nlocProcs)
    {
        iloc = jloc + myloc;
        if (iloc >= data->nlocs){goto NEXT_LOCATION;}
        k = iobs*data->nlocs + iloc;
        // Get the Green's functions onto the matrix G
        ierr = parmt_utils_setDataOnG(iobs, iloc, npmax, *data, G);
        if (ierr != 0)
        {
            printf("%s: Error setting Greens functions\n", fcnm);
            return -1;
        }
        // Tally the objective function
        ierr = parmt_mtSearchL164f(mtloc.ldm, mtloc.nmt,
                                   npts, blockSize,
                                   nlags, lwantLags,
                                   lrescale,
                                   &G[0*npmax], &G[1*npmax], &G[2*npmax],
                                   &G[3*npmax], &G[4*npmax], &G[5*npmax],
                                   CeInv, mtloc.mts, d,
                                   phiLoc, varLoc, lagLoc);
        if (ierr != 0)
        {
            printf("%s: Error calling XC64f search %d %d\n",
                   fcnm, myloc, mtloc.myid);
            return -1;
        }
        // Get the results on the master
        jndx = 0;
        if (mtloc.myid == master){jndx = iloc*mtloc.nmtAll;}
        if (mtloc.commSize == 1)
        {
            array_copy64f_work(mtloc.nmt, phiLoc, &phiWork[jndx]);
            if (lwantLags)
            {
                array_copy32i_work(mtloc.nmt, lagLoc, &lagWork[jndx]);
            }
        }
        else
        {
            MPI_Gatherv(phiLoc, mtloc.nmt, MPI_DOUBLE,
                        &phiWork[jndx], mtloc.nmtProc, mtloc.offset,
                        MPI_DOUBLE, master, mtloc.comm);
            if (lwantLags)
            {
                MPI_Gatherv(lagLoc, mtloc.nmt, MPI_INT,
                            &lagWork[jndx], mtloc.nmtProc, mtloc.offset,
                            MPI_INT, master, mtloc.comm);
            }
        }
NEXT_LOCATION:;
    }
    // Reduce the variance onto varWork
    MPI_Reduce(varLoc, varWork, npts, MPI_DOUBLE, MPI_MIN,
               master, mtloc.comm);
    // Have the location masters reduce their result onto the master
    if (mtloc.myid == master)
    {
        if (nlocProcs > 1)
        {
            MPI_Reduce(phiWork, phi, data->nlocs*mtloc.nmtAll,
                       MPI_DOUBLE, MPI_SUM, master, locComm);
            MPI_Reduce(varWork, var, npts,
                       MPI_DOUBLE, MPI_MIN, master, locComm);
            if (lwantLags)
            {
                MPI_Reduce(lagWork, lags, data->nlocs*mtloc.nmtAll,
                           MPI_INT, MPI_SUM, master, locComm);
            }
        }
    }
    cblas_dscal(npts, 1.0/(double) mtloc.nmtAll, varLoc, 1);
    // Free memory
    memory_free32i(&lagLoc);
    memory_free64f(&phiLoc);
    memory_free64f(&varLoc);
    memory_free64f(&d);
    memory_free64f(&G);
    if (nlocProcs > 1)
    {
        memory_free64f(&phiWork);
        memory_free64f(&varWork);
        memory_free32i(&lagWork);
    }
    phiWork = NULL;
    lagWork = NULL;
    // Block until everyone is done
    if (mtloc.myid == master){MPI_Barrier(locComm);}
    MPI_Barrier(mtloc.comm);
    return 0;
}

int parmt_locSearchXC64f(const MPI_Comm locComm,
                         const int iobs, const int blockSize,
                         const int nlags, const bool lwantLags,
                         struct localMT_struct mtloc,
                         struct parmtData_struct *data,
                         double *__restrict__ phi, int *__restrict__ lags)
{
    const char *fcnm = "parmt_locSearchXC64f\0";
    double *G, *d, *phiLoc, *phiWork;
    int *lagLoc, *lagWork, ierr, iloc, jndx, jloc, k,
        myloc, npmax, nprocs, npts;
    const int master = 0;
    ierr = 0;
    myloc =-1;
    if (mtloc.myid == master)
    {
        MPI_Comm_rank(locComm, &myloc);
        MPI_Comm_size(locComm, &nprocs);
    }
    MPI_Bcast(&myloc, 1, MPI_INT, master, mtloc.comm);
    MPI_Bcast(&nprocs, 1, MPI_INT, master, mtloc.comm);
    if (iobs < 0 || iobs >= data->nobs)
    {
        printf("%s: Invalid observation number %d\n", fcnm, iobs);
        return -1;
    }
    if (mtloc.ldm <= 6)
    {
        printf("%s: Error leading dimension must be at least 6\n", fcnm);
        return -1;
    }
    if (mtloc.mts == NULL || phi == NULL)
    {
        if (mtloc.mts == NULL)
        {
            fprintf(stderr, "%s: Error mts is NULL\n", __func__);
        }
        if (phi == NULL){fprintf(stderr, "%s: Error phi is NULL\n", __func__);}
        return -1;
    }
    npmax = getNpMax(iobs, *data, mtloc);
    if (npmax < 1)
    {
        printf("%s: No data points\n", fcnm);
        return -1;
    }
    // Set space and copy data
    G = memory_calloc64f(6*npmax);
    d = memory_calloc64f(npmax);
    phiLoc = memory_calloc64f(mtloc.nmt);
    if (nprocs > 1)
    {
        if (mtloc.myid == master)
        {
            phiWork = memory_calloc64f(data->nlocs*mtloc.nmtAll);
            if (lwantLags)
            {
                lagWork = memory_calloc32i(data->nlocs*mtloc.nmtAll);
            }
        }
        else
        {
             phiWork = memory_calloc64f(1);
             lagWork = memory_calloc32i(1);
        }
    }
    else
    {
        phiWork = phi;
        lagWork = lags;
    }
    if (lwantLags)
    {
        lagLoc = memory_calloc32i(mtloc.nmt);
    }
    else
    {
        lagLoc = memory_calloc32i(1);
    }
    npts = data->data[iobs].npts;
    cblas_dcopy(npts, data->data[iobs].data, 1, d, 1);
    for (jloc=0; jloc<data->nlocs; jloc=jloc+nprocs)
    {
        iloc = jloc + myloc;
        if (iloc >= data->nlocs){goto NEXT_LOCATION;}
        k = iobs*data->nlocs + iloc;
        // Get the Green's functions onto the matrix G
        ierr = parmt_utils_setDataOnG(iobs, iloc, npmax, *data, G);
        if (ierr != 0)
        {
            printf("%s: Error setting Greens functions\n", fcnm);
            return -1;
        }
        // Tally the objective function
        ierr = parmt_mtSearchXC64f(mtloc.ldm, mtloc.nmt,
                                   npts, blockSize,
                                   nlags, lwantLags,
                                   &G[0*npmax], &G[1*npmax], &G[2*npmax],
                                   &G[3*npmax], &G[4*npmax], &G[5*npmax],
                                   mtloc.mts, d, phiLoc, lagLoc);
        if (ierr != 0)
        {
            printf("%s: Error calling XC64f search %d %d\n",
                   fcnm, myloc, mtloc.myid);
            return -1;
        }
        // Get the results on the master
        jndx = 0;
        if (mtloc.myid == master){jndx = iloc*mtloc.nmtAll;}
        MPI_Gatherv(phiLoc, mtloc.nmt, MPI_DOUBLE,
                    &phiWork[jndx], mtloc.nmtProc, mtloc.offset,
                    MPI_DOUBLE, master, mtloc.comm);
        if (lwantLags)
        {
            MPI_Gatherv(lagLoc, mtloc.nmt, MPI_INT,
                        &lagWork[jndx], mtloc.nmtProc, mtloc.offset,
                        MPI_INT, master, mtloc.comm);
        }
NEXT_LOCATION:;
        MPI_Barrier(mtloc.comm);
    }
    // Have the location masters reduce their result onto the master
    if (mtloc.myid == master && nprocs > 1)
    {
        MPI_Reduce(phiWork, phi, data->nlocs*mtloc.nmtAll,
                   MPI_DOUBLE, MPI_SUM, master, locComm);
        if (lwantLags)
        {
            MPI_Reduce(lagWork, lags, data->nlocs*mtloc.nmtAll,
                       MPI_INT, MPI_SUM, master, locComm);
        }
    }
    // Free memory
    memory_free32i(&lagLoc);
    memory_free64f(&phiLoc);
    memory_free64f(&d);
    memory_free64f(&G);
    if (nprocs > 1)
    {
        memory_free64f(&phiWork);
        memory_free32i(&lagWork);
    }
    phiWork = NULL;
    lagWork = NULL;
    return 0;
}

static int getNpMax(const int iobs, const struct parmtData_struct data,
                    const struct localMT_struct mtloc)
{
    const char *fcnm = "getNpMax\0";
    int iloc, k, npmax, npts;
    const int master = 0;
    // Get the workspace size
    npmax = 0;
    npts = data.data[iobs].header.npts;
    if (mtloc.myid == master)
    {
        for (iloc=0; iloc<data.nlocs; iloc++)
        {
            k = iobs*data.nlocs + iloc;
            npmax = MAX(npmax, data.data[iobs].header.npts);
            if (data.data[iobs].header.npts != data.sacGxx[k].header.npts ||
                data.data[iobs].header.npts != data.sacGyy[k].header.npts ||
                data.data[iobs].header.npts != data.sacGzz[k].header.npts ||
                data.data[iobs].header.npts != data.sacGxy[k].header.npts ||
                data.data[iobs].header.npts != data.sacGxz[k].header.npts ||
                data.data[iobs].header.npts != data.sacGyz[k].header.npts)
            {
                printf("%s: Inconsistent sizes\n", fcnm);
                npmax = 0;
                break;
            }
        }
    }
    MPI_Bcast(&npmax, 1, MPI_INTEGER, master, mtloc.comm);
    if (npmax < 1)
    {
        printf("%s: Error no data points\n", fcnm);
        return -1;
    }
    return npmax;
}

/*
static int setDataOnG(const int iobs, const int iloc, const int npmax,
                       const struct parmtData_struct data,
                       double *__restrict__ G)
{
    const char *fcnm = "setDataOnG\0";
    int k;
    if (iobs < 0 || iobs >= data.nobs)
    {
        printf("%s: Error invalid observation\n", fcnm);
        return -1;
    }
    if (iloc < 0 || iloc >= data.nlocs)
    {
        printf("%s: Error invalid location\n", fcnm);
        return -1;
    }
    k = iobs*data.nlocs + iloc;
    if (npmax > data.sacGxx[k].npts || npmax > data.sacGyy[k].npts ||
        npmax > data.sacGzz[k].npts || npmax > data.sacGxy[k].npts ||
        npmax > data.sacGxz[k].npts || npmax > data.sacGyz[k].npts)
    {
        printf("%s: npmax is too small\n", fcnm);
        return -1;
    }
    if (data.sacGxx[k].data == NULL || data.sacGyy[k].data == NULL ||
        data.sacGzz[k].data == NULL || data.sacGxy[k].data == NULL ||
        data.sacGxz[k].data == NULL || data.sacGyz[k].data == NULL)
    {
        printf("%s: Greens function is null\n", fcnm);
        return -1;
    }
    // Get the Green's functions onto the matrix G
    array_zeros64f_work(6*npmax, G);
    cblas_dcopy(data.sacGxx[k].header.npts,
                data.sacGxx[k].data, 1, &G[0*npmax], 1);
    cblas_dcopy(data.sacGyy[k].header.npts,
                data.sacGyy[k].data, 1, &G[1*npmax], 1);
    cblas_dcopy(data.sacGzz[k].header.npts,
                data.sacGzz[k].data, 1, &G[2*npmax], 1);
    cblas_dcopy(data.sacGxy[k].header.npts,
                data.sacGxy[k].data, 1, &G[3*npmax], 1);
    cblas_dcopy(data.sacGxz[k].header.npts,
                data.sacGxz[k].data, 1, &G[4*npmax], 1);
    cblas_dcopy(data.sacGyz[k].header.npts,
                data.sacGyz[k].data, 1, &G[5*npmax], 1);
    return 0;
}
*/

ext: c


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef vbitset_deeaeda9_1fc9_43e4_9026_c5eeb93c6ebb_h
#define vbitset_deeaeda9_1fc9_43e4_9026_c5eeb93c6ebb_h

#include <intrin.h>
#include <gslib/config.h>
#include <gslib/std.h>

__gslib_begin__

inline int bit_count(uint n)
{
    uint t = n - ((n >> 1) & 033333333333) - ((n >> 2) & 011111111111);
    return ((t + (t >> 3)) & 030707070707) % 63;
}

inline uint bitelim_first_nonzero(uint n) { return n & (n - 1); }
inline uint bitmask_first_nonzero(uint n) { return n ^ (n - 1); }

inline int bitscan_forward(uint n)
{
    uint p;
    return _BitScanForward((unsigned long*)&p, n) ? (int)p : -1;
}

inline int bitscan_reverse(uint n)
{
    uint p;
    return _BitScanReverse((unsigned long*)&p, n) ? (int)p : -1;
}

class vbitset
{
public:
    typedef vector<uint> container;
    static const uint knot_size = sizeof(uint) * 8;

protected:
    container           _container;
    int                 _size;

public:
    vbitset() { _size = 0; }
    int knot_count() const { return (int)_container.size(); }
    uint knot_at(int i) const
    {
        assert(i < knot_count());
        return _container.at(i);
    }
    void set_knot(int p, uint k)
    {
        if(p >= knot_count()) {
            assert(!"unexpected.");
            return;
        }
        _container.at(p) = k;
    }
    uint& ref_knot_at(int i)
    {
        assert(i < knot_count());
        return _container.at(i);
    }
    int capacity() const { return knot_size * knot_count(); }
    int size() const { return _size; }
    int length() const { return _size; }
    void set(int pos) { reset(pos, true); }
    void reset(int pos, bool b = false)
    {
        int cc = cycle_count(pos);
        int cidx = cycle_index(pos);
        ensure_container_capacity(pos, cc, cidx);
        uint& n = _container.at(cc);
        uint mask = (1 << cidx);
        b ? (n |= mask) : (n &= ~mask);
    }
    void reset(bool b)
    {
        int cc = cycle_count(_size);
        int cidx = cycle_index(_size);
        int m = b ? -1 : 0;
        memset(&_container.front(), m, cc * sizeof(uint));
        uint mask = (1 << (cidx + 1)) - 1;
        uint& n = _container.at(cc);
        b ? (n |= mask) : (n &= ~mask);
    }
    bool test(int pos)
    {
        if(pos >= _size)
            return false;
        int cc = cycle_count(pos);
        int cidx = cycle_index(pos);
        uint n = knot_at(cc);
        uint mask = (1 << cidx);
        return (n & mask) != 0;
    }
    void resize(int s)
    {
        if(s == 0) {
            _container.clear();
            _size = 0;
            return;
        }
        if(s == _size)
            return;
        else if(s < _size) {
            ensure_container_capacity(s - 1);
            return;
        }
        else {
            int cc = cycle_count(s - 1);
            _container.resize(cc + 1);
            _size = s;
        }
    }
    void rotate_left(int bias)
    {
        if(!bias)
            return;
        if(bias < 0)
            return rotate_right(-bias);
        bias %= _size;
        if(!bias)
            return;
        int cc = cycle_count(bias);
        int cm = cycle_index(bias);
        rotate_left(cc, cm);
    }
    void rotate_right(int bias)
    {
        if(!bias)
            return;
        if(bias < 0)
            return rotate_left(-bias);
        bias %= _size;
        if(!bias)
            return;
        int cc = cycle_count(bias);
        int cm = cycle_index(bias);
        rotate_right(cc, cm);
    }
    // todo: following
    void rotate_left(int cc, int cm)
    {
        assert(cc >= 0 && cm >= 0);
        if(!cm) {
            uint* ptr = new uint[cc];
            int crst = (int)_container.size() - cc;
            memcpy_s(ptr, sizeof(uint) * cc, &_container.front(), _container.capacity());
            memcpy_s(&_container.front(), sizeof(uint) * crst, &_container.at(cc), sizeof(uint) * crst);
            memcpy_s(&_container.front(), sizeof(uint) * cc, ptr, sizeof(uint) * cc);
            delete [] ptr;
        }
        else {

        }
    }
    void rotate_right(int cc, int cm)
    {
        assert(cc >= 0 && cm >= 0);
    }
    void shift_left(int bias);
    void shift_right(int bias);
    void shift_left_1(int bias);
    void shift_right_1(int bias);

public:
    static int cycle_count(int p) { return p / knot_size; }
    static int cycle_index(int p) { return p % knot_size; }

protected:
    void ensure_container_capacity(int pos)
    {
        if(pos < _size)
            return;
        int cc = cycle_count(pos);
        int cidx = cycle_index(pos);
        ensure_container_capacity(pos, cc, cidx);
    }
    void ensure_container_capacity(int pos, int cc, int cidx)
    {
        assert(cc * knot_size + cidx == pos);
        if(pos < _size)
            return;
        int kc = knot_count();
        assert(kc <= cc);
        if(kc == cc) {
            _container.push_back(0);
            _size = pos + 1;
        }
        else {
            cc ++;
            _container.resize(cc);
            int delta = cc - kc;
            memset(&_container.at(kc), 0, delta * sizeof(uint));
            _size = pos + 1;
        }
    }
};

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
 *  Created on: Oct. 19, 2019
 *      Author: Govind
 */

/*
 * standard includes */
	#include <stdio.h>
	#include <stdlib.h>

/*
 *  includes for GSL components
 *  	- use double precision
 */
	#include <gsl/gsl_vector_double.h>
	#include <gsl/gsl_matrix_double.h>
	#include <gsl/gsl_rng.h>
	#include <gsl/gsl_randist.h>

/*
 * FUNCTIONS
 */

/*
 * simple Fibonacci sequence generator function, using recursion
 * */
size_t fib(size_t k)
{
	if (k==0)
	{
		return 0;
	}
	else if (k==1)
	{
		return 1;
	}
	else /* k >= 2 */
	{
		return fib(k-1) + fib(k-2);
	}
}

gsl_matrix *embt_mm(const gsl_matrix *U, const gsl_matrix *V, size_t N)
{
	gsl_matrix *W = gsl_matrix_alloc(N,N);
	double dp;
	double uk,vk;

	for (size_t i=0; i != N; ++i)
	{
		for (size_t j=0; j != N; ++j)
		{
			/* compute element (i,j) of W */
			dp = 0;
			for (size_t k=0; k != N; ++k)
			{
				uk = gsl_matrix_get(U,i,k);
				vk = gsl_matrix_get(V,k,j);
				dp += uk*vk;
			}
			gsl_matrix_set(W,i,j,dp);
		}
	}
	return W;
}

gsl_vector gsl_print_vector(const gsl_vector V, size_t N)
{
	  /* Print the vector*/
	   for (size_t i = 0; i < N; i++)
	   {
	       printf ("v_%d = %g\n", i, gsl_vector_get (V, i));
	   }
}
gsl_matrix embt_print_matrix(const gsl_matrix V, size_t N)
{
	printf("Matrix m\n");
	  for (size_t i=0;i!=N;i++)
	    {
	      for (size_t j=0;j!=N;j++)
		{
		  printf("%f ",gsl_matrix_get(V,i,j));
		}
	      printf("\n");
	    }
	  printf("\n");
}
int main()
{
	/*
	 * INITIALIZE PARAMETERS
	 */
		/* vectors parameters */
		size_t		N=10; /* index type, vector sizes */
		gsl_vector *a = gsl_vector_alloc(N); /* allocate vector from heap of size N */
		gsl_vector *b = gsl_vector_alloc(N); /* allocate vector from heap of size N */
		gsl_vector *c = gsl_vector_calloc(N); /* allocate vector of size N but initialize entries to zero */

		/* random number generator parameters */
		const gsl_rng_type *T;
		gsl_rng *r; /* handle for our random number generator */

		/* matrix parameters */
		gsl_matrix *A = gsl_matrix_alloc(N,N);
		gsl_matrix *B = gsl_matrix_alloc(N,N);
		gsl_matrix *C = gsl_matrix_calloc(N,N);

	/*
	 * SET UP RANDOM NUMBER GENERATION
	 */
		gsl_rng_env_setup();
		T = gsl_rng_default;
		r = gsl_rng_alloc(T);

	/*
	 *  VECTOR OPERATIONS
	 */
		/* set the vector elements */
		for (size_t i = 0; i != N; ++i)
		{
			gsl_vector_set(a, i, fib(i)); /* set element i of vector a to Fibonacci number i */
			gsl_vector_set(b, i, gsl_ran_flat(r,-1.0,+1.0)); /* set element of vector b to random number */
		}

		/* c = a + b */
		gsl_vector_add(c, a); /* c += a */
		gsl_vector_add(c, b); /* c += b */

		/* print results */
		for (size_t i = 0; i != N; ++i)
		{
			printf("i=%d, a(i)=%f, b(i)=%f, c(i)=%f\n", i,
			gsl_vector_get(a, i),
			gsl_vector_get(b, i),
			gsl_vector_get(c, i));
		}

	/*
	 *	MATRIX OPERATIONS - your homework!! :)
	 */
		/* fill A with first N*N Fibonacci numbers, starting with row 1 (cols 1-10), then row 2, etc. */
		for (size_t i=0; i != N; ++i)
		{
			for (size_t j = 0; j != N; ++j)
			{
				gsl_matrix_set(A, i, j, (double) fib(j+i*N));
			}
		}

		/* fill B with N*N random numbers, uniformly distributed over the interval (-100, 100) */
		for (size_t i=0; i != N; ++i)
		{
			for (size_t j = 0; j != N; ++j)
			{
				gsl_matrix_set(B, i, j, gsl_ran_flat(r,-100.0,+100.0));
			}
		}

		/* make C the product of A and B */
		C = embt_mm(A,B,N);

		/* print the results */
		embt_print_vector(a,N);
		embt_print_vector(b,N);
		embt_print_vector(c,N);
		embt_print_matrix(A,N);
		embt_print_matrix(B,N);
		embt_print_matrix(C,N);

	/* de-allocate memory */
	gsl_vector_free(a);
	gsl_vector_free(b);
	gsl_vector_free(c);
	gsl_matrix_free(A);
	gsl_matrix_free(B);
	gsl_matrix_free(C);
	return EXIT_SUCCESS;
}

max_stars_repo_path: C program/DSP/GSL./Matrix_mul.c

max_issues_repo_path: C program/DSP/GSL./Matrix_mul.c

max_forks_repo_path: C program/DSP/GSL./Matrix_mul.c


text:
/*
 * Gsl_fcmp: test gsl_fcmp (float comparison)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_fcmp.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_fcmp.c $(pkg-config --libs gsl)
 */

#include <gsl/gsl_sys.h>
#include <stdio.h>

void main(void) {

    double n1, n2, eps;

    n1 = 1.0;
    n2 = 1.1;

    eps = 0.2;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    eps = 0.1;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    eps = 0.05;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 0.01;
    n2 = 0.0;
    eps = 0.1;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 0.01;
    n2 = 0.000001;
    eps = 0.1;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 1.01;
    n2 = 1.000001;
    eps = 0.1;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 0.51;
    n2 = 0.500001;
    eps = 0.1;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 0.51;
    n2 = 0.500001;
    eps = 0.01;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

    n1 = 0.51;
    n2 = 0.5;
    eps = 0.02;

    printf("%f ~ %f (+/- %f) : %d\n", n1, n2, eps, gsl_fcmp(n1, n2, eps));

}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
/* fit/gsl_fit.h
 *
 * Copyright (C) 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_FIT_H__
#define __GSL_FIT_H__

#include <stdlib.h>
#include <gsl/gsl_math.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

int gsl_fit_linear (const double * x, const size_t xstride,
                    const double * y, const size_t ystride,
                    const size_t n,
                    double * c0, double * c1,
                    double * cov00, double * cov01, double * cov11,
                    double * sumsq);

int gsl_fit_wlinear (const double * x, const size_t xstride,
                     const double * w, const size_t wstride,
                     const double * y, const size_t ystride,
                     const size_t n,
                     double * c0, double * c1,
                     double * cov00, double * cov01, double * cov11,
                     double * chisq);

int
gsl_fit_linear_est (const double x,
                    const double c0, const double c1,
                    const double c00, const double c01, const double c11,
                    double *y, double *y_err);

int gsl_fit_mul (const double * x, const size_t xstride,
                 const double * y, const size_t ystride,
                 const size_t n,
                 double * c1,
                 double * cov11,
                 double * sumsq);

int gsl_fit_wmul (const double * x, const size_t xstride,
                  const double * w, const size_t wstride,
                  const double * y, const size_t ystride,
                  const size_t n,
                  double * c1,
                  double * cov11,
                  double * sumsq);

int
gsl_fit_mul_est (const double x,
                 const double c1,
                 const double c11,
                 double *y, double *y_err);

/* choose better names!! */

int gsl_fit_poly (const double * x,
                  const double * w,
                  const double * y,
                  size_t n,
                  double * c, size_t m,
                  double * chisq);

int gsl_fit_fns (const double * A,
                 const double * w,
                 const double * y,
                 size_t n,
                 double * c, size_t m,
                 double * chisq);

int gsl_fit_linear_nd (double * m, double * y, double * w);

__END_DECLS

#endif /* __GSL_FIT_H__ */


text:
/* fit/linear.c
 *
 * Copyright (C) 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fit.h>

/* Fit the data (x_i, y_i) to the linear relationship

   Y = c0 + c1 x

   returning,

   c0, c1  --  coefficients
   cov00, cov01, cov11  --  variance-covariance matrix of c0 and c1,
   sumsq   --   sum of squares of residuals

   This fit can be used in the case where the errors for the data are
   uknown, but assumed equal for all points. The resulting
   variance-covariance matrix estimates the error in the coefficients
   from the observed variance of the points around the best fit line.
*/

int
gsl_fit_linear (const double *x, const size_t xstride,
                const double *y, const size_t ystride,
                const size_t n,
                double *c0, double *c1,
                double *cov_00, double *cov_01, double *cov_11, double *sumsq)
{
  double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;

  size_t i;

  for (i = 0; i < n; i++)
    {
      m_x += (x[i * xstride] - m_x) / (i + 1.0);
      m_y += (y[i * ystride] - m_y) / (i + 1.0);
    }

  for (i = 0; i < n; i++)
    {
      const double dx = x[i * xstride] - m_x;
      const double dy = y[i * ystride] - m_y;

      m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
      m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
    }

  /* In terms of y = a + b x */

  {
    double s2 = 0, d2 = 0;
    double b = m_dxdy / m_dx2;
    double a = m_y - m_x * b;

    *c0 = a;
    *c1 = b;

    /* Compute chi^2 = \sum (y_i - (a + b * x_i))^2 */

    for (i = 0; i < n; i++)
      {
        const double dx = x[i * xstride] - m_x;
        const double dy = y[i * ystride] - m_y;
        const double d = dy - b * dx;
        d2 += d * d;
      }

    s2 = d2 / (n - 2.0);        /* chisq per degree of freedom */

    *cov_00 = s2 * (1.0 / n) * (1 + m_x * m_x / m_dx2);
    *cov_11 = s2 * 1.0 / (n * m_dx2);

    *cov_01 = s2 * (-m_x) / (n * m_dx2);

    *sumsq = d2;
  }

  return GSL_SUCCESS;
}

/* Fit the weighted data (x_i, w_i, y_i) to the linear relationship

   Y = c0 + c1 x

   returning,

   c0, c1  --  coefficients
   s0, s1  --  the standard deviations of c0 and c1,
   r       --  the correlation coefficient between c0 and c1,
   chisq   --  weighted sum of squares of residuals */

int
gsl_fit_wlinear (const double *x, const size_t xstride,
                 const double *w, const size_t wstride,
                 const double *y, const size_t ystride,
                 const size_t n,
                 double *c0, double *c1,
                 double *cov_00, double *cov_01, double *cov_11,
                 double *chisq)
{

  /* compute the weighted means and weighted deviations from the means */

  /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */

  double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;

  size_t i;

  for (i = 0; i < n; i++)
    {
      const double wi = w[i * wstride];

      if (wi > 0)
        {
          W += wi;
          wm_x += (x[i * xstride] - wm_x) * (wi / W);
          wm_y += (y[i * ystride] - wm_y) * (wi / W);
        }
    }

  W = 0;                        /* reset the total weight */

  for (i = 0; i < n; i++)
    {
      const double wi = w[i * wstride];

      if (wi > 0)
        {
          const double dx = x[i * xstride] - wm_x;
          const double dy = y[i * ystride] - wm_y;

          W += wi;
          wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
          wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
        }
    }

  /* In terms of y = a + b x */

  {
    double d2 = 0;
    double b = wm_dxdy / wm_dx2;
    double a = wm_y - wm_x * b;

    *c0 = a;
    *c1 = b;

    *cov_00 = (1 / W) * (1 + wm_x * wm_x / wm_dx2);
    *cov_11 = 1 / (W * wm_dx2);

    *cov_01 = -wm_x / (W * wm_dx2);

    /* Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2 */

    for (i = 0; i < n; i++)
      {
        const double wi = w[i * wstride];

        if (wi > 0)
          {
            const double dx = x[i * xstride] - wm_x;
            const double dy = y[i * ystride] - wm_y;
            const double d = dy - b * dx;
            d2 += wi * d * d;
          }
      }

    *chisq = d2;
  }

  return GSL_SUCCESS;
}

int
gsl_fit_linear_est (const double x,
                    const double c0, const double c1,
                    const double c00, const double c01, const double c11,
                    double *y, double *y_err)
{
  *y = c0 + c1 * x;
  *y_err = sqrt (c00 + x * (2 * c01 + c11 * x));
  return GSL_SUCCESS;
}

int
gsl_fit_mul (const double *x, const size_t xstride,
             const double *y, const size_t ystride,
             const size_t n,
             double *c1, double *cov_11, double *sumsq)
{
  double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;

  size_t i;

  for (i = 0; i < n; i++)
    {
      m_x += (x[i * xstride] - m_x) / (i + 1.0);
      m_y += (y[i * ystride] - m_y) / (i + 1.0);
    }

  for (i = 0; i < n; i++)
    {
      const double dx = x[i * xstride] - m_x;
      const double dy = y[i * ystride] - m_y;

      m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
      m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
    }

  /* In terms of y =  b x */

  {
    double s2 = 0, d2 = 0;
    double b = (m_x * m_y + m_dxdy) / (m_x * m_x + m_dx2);

    *c1 = b;

    /* Compute chi^2 = \sum (y_i -  b * x_i)^2 */

    for (i = 0; i < n; i++)
      {
        const double dx = x[i * xstride] - m_x;
        const double dy = y[i * ystride] - m_y;
        const double d = (m_y - b * m_x) + dy - b * dx;
        d2 += d * d;
      }

    s2 = d2 / (n - 1.0);        /* chisq per degree of freedom */

    *cov_11 = s2 * 1.0 / (n * (m_x * m_x + m_dx2));

    *sumsq = d2;
  }

  return GSL_SUCCESS;
}

int
gsl_fit_wmul (const double *x, const size_t xstride,
              const double *w, const size_t wstride,
              const double *y, const size_t ystride,
              const size_t n,
              double *c1, double *cov_11, double *chisq)
{

  /* compute the weighted means and weighted deviations from the means */

  /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */

  double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;

  size_t i;

  for (i = 0; i < n; i++)
    {
      const double wi = w[i * wstride];

      if (wi > 0)
        {
          W += wi;
          wm_x += (x[i * xstride] - wm_x) * (wi / W);
          wm_y += (y[i * ystride] - wm_y) * (wi / W);
        }
    }

  W = 0;                        /* reset the total weight */

  for (i = 0; i < n; i++)
    {
      const double wi = w[i * wstride];

      if (wi > 0)
        {
          const double dx = x[i * xstride] - wm_x;
          const double dy = y[i * ystride] - wm_y;

          W += wi;
          wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
          wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
        }
    }

  /* In terms of y = b x */

  {
    double d2 = 0;
    double b = (wm_x * wm_y + wm_dxdy) / (wm_x * wm_x + wm_dx2);

    *c1 = b;

    *cov_11 = 1 / (W * (wm_x * wm_x + wm_dx2));

    /* Compute chi^2 = \sum w_i (y_i - b * x_i)^2 */

    for (i = 0; i < n; i++)
      {
        const double wi = w[i * wstride];

        if (wi > 0)
          {
            const double dx = x[i * xstride] - wm_x;
            const double dy = y[i * ystride] - wm_y;
            const double d = (wm_y - b * wm_x) + (dy - b * dx);
            d2 += wi * d * d;
          }
      }

    *chisq = d2;
  }

  return GSL_SUCCESS;
}

int
gsl_fit_mul_est (const double x,
                 const double c1, const double c11,
                 double *y, double *y_err)
{
  *y = c1 * x;
  *y_err = sqrt (c11) * fabs (x);
  return GSL_SUCCESS;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once
#include "iarglist.h"
#include "optioninfo.h"
#include "configaccess.h"
#include "format.h"
#include "streamreader.h"
#include <gsl/gsl>
#include <cmdlime/errors.h>
#include <cmdlime/customnames.h>
#include <vector>
#include <sstream>
#include <functional>
#include <memory>

namespace cmdlime::detail{

template <typename T>
class ArgList : public IArgList{
public:
    ArgList(std::string name,
            std::string type,
            std::function<std::vector<T>&()> argListGetter)
        : info_(std::move(name), {}, std::move(type))
        , argListGetter_(std::move(argListGetter))
    {
    }

    void setDefaultValue(const std::vector<T>& value)
    {
        hasValue_ = true;
        defaultValue_ = value;
    }

    OptionInfo& info() override
    {
        return info_;
    }

    const OptionInfo& info() const override
    {
        return info_;
    }

private:
    bool read(const std::string& data) override
    {
        if (!isDefaultValueOverwritten_){
            argListGetter_().clear();
            isDefaultValueOverwritten_ = true;
        }
        auto stream = std::stringstream{data};
        argListGetter_().emplace_back();
        if (!readFromStream(stream, argListGetter_().back()))
            return false;
        hasValue_ = true;
        return true;
    }

    bool hasValue() const override
    {
        return hasValue_;
    }

    bool isOptional() const override
    {
        return defaultValue_.has_value();
    }

    std::string defaultValue() const override
    {
        if (!defaultValue_.has_value())
            return {};
        auto stream = std::stringstream{};
        stream << "{";
        auto firstVal = true;
        for (auto& val : defaultValue_.value()){
            if (firstVal)
                stream << val;
            else
                stream << ", " << val;
            firstVal = false;
        }
        stream << "}";
        return stream.str();
    }

private:
    OptionInfo info_;
    std::function<std::vector<T>&()> argListGetter_;
    bool hasValue_ = false;
    std::optional<std::vector<T>> defaultValue_;
    bool isDefaultValueOverwritten_ = false;
};

template <>
inline bool ArgList<std::string>::read(const std::string& data)
{
    argListGetter_().push_back(data);
    hasValue_ = true;
    return true;
}

template<typename T, typename TConfig>
class ArgListCreator{
    using NameProvider = typename Format<ConfigAccess<TConfig>::format()>::nameProvider;

public:
    ArgListCreator(TConfig& cfg,
                   const std::string& varName,
                   const std::string& type,
                   std::function<std::vector<T>&()> argListGetter)
        : cfg_(cfg)
    {
        Expects(!varName.empty());
        Expects(!type.empty());
        argList_ = std::make_unique<ArgList<T>>(NameProvider::fullName(varName),
                                                NameProvider::valueName(type),
                                                std::move(argListGetter));
    }

    ArgListCreator<T, TConfig>& operator<<(const std::string& info)
    {
        argList_->info().addDescription(info);
        return *this;
    }

    ArgListCreator<T, TConfig>& operator<<(const Name& customName)
    {
        argList_->info().resetName(customName.value());
        return *this;
    }

    ArgListCreator<T, TConfig>& operator<<(const ValueName& valueName)
    {
        argList_->info().resetValueName(valueName.value());
        return *this;
    }

    ArgListCreator<T, TConfig>& operator()(std::vector<T> defaultValue = {})
    {
        defaultValue_ = std::move(defaultValue);
        argList_->setDefaultValue(defaultValue_);
        return *this;
    }

    operator std::vector<T>()
    {
        ConfigAccess<TConfig>{cfg_}.setArgList(std::move(argList_));
        return defaultValue_;
    }

private:
    std::unique_ptr<ArgList<T>> argList_;
    std::vector<T> defaultValue_;
    TConfig& cfg_;
};

template <typename T, typename TConfig>
ArgListCreator<T, TConfig> makeArgListCreator(TConfig& cfg,
                                              const std::string& varName,
                                              const std::string& type,
                                              std::function<std::vector<T>&()> argListGetter)
{
    return ArgListCreator<T, TConfig>{cfg, varName, type, std::move(argListGetter)};
}

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_multiroots.h>

BC_INLINE3(GSL_MULTIROOT_FN_EVAL,gsl_multiroot_function*,gsl_vector*,gsl_vector*,int)
BC_INLINE3(GSL_MULTIROOT_FN_EVAL_F,gsl_multiroot_function_fdf*,gsl_vector*,gsl_vector*,int)
BC_INLINE3(GSL_MULTIROOT_FN_EVAL_DF,gsl_multiroot_function_fdf*,gsl_vector*,gsl_matrix*,int)
BC_INLINE4(GSL_MULTIROOT_FN_EVAL_F_DF,gsl_multiroot_function_fdf*,gsl_vector*,gsl_vector*,gsl_matrix*,int)


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include <string.h>

struct quadratic_params
{
    double a, b, c;
};

double
quadratic (double x, void *params)
{
    struct quadratic_params *p
            = (struct quadratic_params *) params;

    double a = p->a;
    double b = p->b;
    double c = p->c;

    return (a * x + b) * x + c;
}

double
quadratic_deriv (double x, void *params)
{
    struct quadratic_params *p
            = (struct quadratic_params *) params;

    double a = p->a;
    double b = p->b;

    return 2.0 * a * x + b;
}

void
quadratic_fdf (double x, void *params,
               double *y, double *dy)
{
    struct quadratic_params *p
            = (struct quadratic_params *) params;

    double a = p->a;
    double b = p->b;
    double c = p->c;

    *y = (a * x + b) * x + c;
    *dy = 2.0 * a * x + b;
}

int
main (int argc, char *argv[])
{
    int status;
    int iter = 0, max_iter = 100;
    const gsl_root_fdfsolver_type *T;
    gsl_root_fdfsolver *s;
    double x0, x = 5.0, r_expected = sqrt (5.0);
    gsl_function_fdf FDF;
    struct quadratic_params params = {1.0, -2.0, 1.0};

    FDF.f = &quadratic;
    FDF.df = &quadratic_deriv;
    FDF.fdf = &quadratic_fdf;
    FDF.params = &params;

    if (argc < 2) {
        errno = EINVAL;
        perror("");
        return errno;
    }
    char *method = argv[1];
    if (0 == strcmp("newton", method)) {
        T = gsl_root_fdfsolver_newton;
    } else if (0 == strcmp("secant", method)) {
        T = gsl_root_fdfsolver_secant;
    } else if (0 == strcmp("steffenson", method)) {
        T = gsl_root_fdfsolver_steffenson;
    } else {
        errno = EINVAL;
        perror("");
        return errno;
    }
    s = gsl_root_fdfsolver_alloc (T);
    gsl_root_fdfsolver_set (s, &FDF, x);

    printf ("using %s method\n",
            gsl_root_fdfsolver_name (s));

    printf ("%-5s %10s %10s %10s\n",
            "iter", "root", "err", "err(est)");
    do
    {
        iter++;
        status = gsl_root_fdfsolver_iterate (s);
        x0 = x;
        x = gsl_root_fdfsolver_root (s);
        status = gsl_root_test_delta (x, x0, 0, 1e-3);

        if (status == GSL_SUCCESS)
            printf ("Converged:\n");

        printf ("%5d %10.7f %+10.7f %10.7f\n",
                iter, x, x - r_expected, x - x0);
    }
    while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fdfsolver_free (s);
    return status;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


text:
/* specfunc/airy_der.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_airy.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval_mode.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC aide.f, bide.f, aid.f, bid.f, r9admp.f */

/*
 series for aif on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   5.22e-18
                                         log weighted error  17.28
                               significant figures required  16.01
                                    decimal places required  17.73
*/
static double aif_data[8] = {
   0.10527461226531408809,
   0.01183613628152997844,
   0.00012328104173225664,
   0.00000062261225638140,
   0.00000000185298887844,
   0.00000000000363328873,
   0.00000000000000504622,
   0.00000000000000000522
};
static cheb_series aif_cs = {
  aif_data,
  7,
  -1, 1,
  7
};

/*
 series for aig on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   3.14e-19
                                         log weighted error  18.50
                               significant figures required  17.44
                                    decimal places required  18.98
*/
static double aig_data[9] = {
   0.021233878150918666852,
   0.086315930335214406752,
   0.001797594720383231358,
   0.000014265499875550693,
   0.000000059437995283683,
   0.000000000152403366479,
   0.000000000000264587660,
   0.000000000000000331562,
   0.000000000000000000314
};
static cheb_series aig_cs = {
  aig_data,
  8,
  -1, 1,
  8
};

/*
 series for aip2 on the interval  0.00000e+00 to  1.25000e-01
                                        with weighted error   2.15e-17
                                         log weighted error  16.67
                               significant figures required  14.27
                                    decimal places required  17.26
*/
static double aip2_data[15] = {
    0.0065457691989713757,
    0.0023833724120774592,
   -0.0000430700770220586,
    0.0000015629125858629,
   -0.0000000815417186163,
    0.0000000054103738057,
   -0.0000000004284130883,
    0.0000000000389497963,
   -0.0000000000039623161,
    0.0000000000004428184,
   -0.0000000000000536297,
    0.0000000000000069650,
   -0.0000000000000009620,
    0.0000000000000001403,
   -0.0000000000000000215
};
static cheb_series aip2_cs = {
  aip2_data,
  14,
  -1, 1,
  9
};

/*
 series for aip1 on the interval  1.25000e-01 to  1.00000e+00
                                        with weighted error   2.60e-17
                                         log weighted error  16.58
                               significant figures required  14.91
                                    decimal places required  17.28
*/
static double aip1_data[25] = {
    0.0358865097808301538,
    0.0114668575627764899,
   -0.0007592073583861400,
    0.0000869517610893841,
   -0.0000128237294298592,
    0.0000022062695681038,
   -0.0000004222295185921,
    0.0000000874686415726,
   -0.0000000192773588418,
    0.0000000044668460054,
   -0.0000000010790108052,
    0.0000000002700029447,
   -0.0000000000696480108,
    0.0000000000184489907,
   -0.0000000000050027817,
    0.0000000000013852243,
   -0.0000000000003908218,
    0.0000000000001121536,
   -0.0000000000000326862,
    0.0000000000000096619,
   -0.0000000000000028935,
    0.0000000000000008770,
   -0.0000000000000002688,
    0.0000000000000000832,
   -0.0000000000000000260
};
static cheb_series aip1_cs = {
  aip1_data,
  24,
  -1, 1,
  14
};

/*
 series for bif on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   9.05e-18
                                         log weighted error  17.04
                               significant figures required  15.83
                                    decimal places required  17.49
*/
static double bif_data[8] = {
   0.1153536790828570243,
   0.0205007894049192875,
   0.0002135290278902876,
   0.0000010783960614677,
   0.0000000032094708833,
   0.0000000000062930407,
   0.0000000000000087403,
   0.0000000000000000090
};
static cheb_series bif_cs = {
  bif_data,
  7,
  -1, 1,
  7
};

/*
 series for big on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   5.44e-19
                                         log weighted error  18.26
                               significant figures required  17.46
                                    decimal places required  18.74
*/
static double big_data[9] = {
   -0.097196440416443537390,
    0.149503576843167066571,
    0.003113525387121326042,
    0.000024708570579821297,
    0.000000102949627731379,
    0.000000000263970373987,
    0.000000000000458279271,
    0.000000000000000574283,
    0.000000000000000000544
};
static cheb_series big_cs = {
  big_data,
  8,
  -1, 1,
  8
};

/*
 series for bif2 on the interval  1.00000e+00 to  8.00000e+00
                                        with weighted error   3.82e-19
                                         log weighted error  18.42
                               significant figures required  17.68
                                    decimal places required  18.92
*/
static double bif2_data[10] = {
   0.323493987603522033521,
   0.086297871535563559139,
   0.002994025552655397426,
   0.000051430528364661637,
   0.000000525840250036811,
   0.000000003561751373958,
   0.000000000017146864007,
   0.000000000000061663520,
   0.000000000000000171911,
   0.000000000000000000382
};
static cheb_series bif2_cs = {
  bif2_data,
  9,
  -1, 1,
  9
};

/*
 series for big2 on the interval  1.00000e+00 to  8.00000e+00
                                        with weighted error   3.35e-17
                                         log weighted error  16.48
                               significant figures required  16.52
                                    decimal places required  16.98
*/
static double big2_data[10] = {
   1.6062999463621294578,
   0.7449088819876088652,
   0.0470138738610277380,
   0.0012284422062548239,
   0.0000173222412256624,
   0.0000001521901652368,
   0.0000000009113560249,
   0.0000000000039547918,
   0.0000000000000130017,
   0.0000000000000000335
};
static cheb_series big2_cs = {
  big2_data,
  9,
  -1, 1,
  9
};

/*
 series for bip2 on the interval  0.00000e+00 to  1.25000e-01
                                        with weighted error   2.07e-18
                                         log weighted error  17.69
                               significant figures required  16.51
                                    decimal places required  18.42
*/
static double bip2_data[29] = {
    -0.13269705443526630495,
    -0.00568443626045977481,
    -0.00015643601119611610,
    -0.00001136737203679562,
    -0.00000143464350991284,
    -0.00000018098531185164,
     0.00000000926177343611,
     0.00000001710005490721,
     0.00000000476698163504,
    -0.00000000035195022023,
    -0.00000000058890614316,
    -0.00000000006678499608,
     0.00000000006395565102,
     0.00000000001554529427,
    -0.00000000000792397000,
    -0.00000000000258326243,
     0.00000000000121655048,
     0.00000000000038707207,
    -0.00000000000022487045,
    -0.00000000000004953477,
     0.00000000000004563782,
     0.00000000000000332998,
    -0.00000000000000921750,
     0.00000000000000094157,
     0.00000000000000167154,
    -0.00000000000000055134,
    -0.00000000000000022369,
     0.00000000000000017487,
     0.00000000000000000207
};
static cheb_series bip2_cs = {
  bip2_data,
  28,
  -1, 1,
  14
};

/*
 series for bip1 on the interval  1.25000e-01 to  3.53553e-01
                                        with weighted error   1.86e-17
                                         log weighted error  16.73
                               significant figures required  15.67
                                    decimal places required  17.42
*/
static double bip1_data[24] = {
   -0.1729187351079553719,
   -0.0149358492984694364,
   -0.0005471104951678566,
    0.0001537966292958408,
    0.0000154353476192179,
   -0.0000065434113851906,
    0.0000003728082407879,
    0.0000002072078388189,
   -0.0000000658173336470,
    0.0000000074926746354,
    0.0000000011101336884,
   -0.0000000007265140553,
    0.0000000001782723560,
   -0.0000000000217346352,
   -0.0000000000020302035,
    0.0000000000019311827,
   -0.0000000000006044953,
    0.0000000000001209450,
   -0.0000000000000125109,
   -0.0000000000000019917,
    0.0000000000000015154,
   -0.0000000000000004977,
    0.0000000000000001155,
   -0.0000000000000000186
};
static cheb_series bip1_cs = {
  bip1_data,
  23,
  -1, 1,
  13
};

/*
 series for an22 on the interval -1.00000e+00 to -1.25000e-01
                                        with weighted error   3.30e-17
                                         log weighted error  16.48
                               significant figures required  14.95
                                    decimal places required  17.24
*/
static double an22_data[33] = {
    0.0537418629629794329,
   -0.0126661435859883193,
   -0.0011924334106593007,
   -0.0002032327627275655,
   -0.0000446468963075164,
   -0.0000113359036053123,
   -0.0000031641352378546,
   -0.0000009446708886149,
   -0.0000002966562236472,
   -0.0000000969118892024,
   -0.0000000326822538653,
   -0.0000000113144618964,
   -0.0000000040042691002,
   -0.0000000014440333684,
   -0.0000000005292853746,
   -0.0000000001967763374,
   -0.0000000000740800096,
   -0.0000000000282016314,
   -0.0000000000108440066,
   -0.0000000000042074801,
   -0.0000000000016459150,
   -0.0000000000006486827,
   -0.0000000000002574095,
   -0.0000000000001027889,
   -0.0000000000000412846,
   -0.0000000000000166711,
   -0.0000000000000067657,
   -0.0000000000000027585,
   -0.0000000000000011296,
   -0.0000000000000004645,
   -0.0000000000000001917,
   -0.0000000000000000794,
   -0.0000000000000000330
};
static cheb_series an22_cs = {
  an22_data,
  32,
  -1, 1,
  18
};

/*
 series for an21 on the interval -1.25000e-01 to -1.56250e-02
                                        with weighted error   3.43e-17
                                         log weighted error  16.47
                               significant figures required  14.48
                                    decimal places required  17.16
*/
static double an21_data[24] = {
    0.0198313155263169394,
   -0.0029376249067087533,
   -0.0001136260695958196,
   -0.0000100554451087156,
   -0.0000013048787116563,
   -0.0000002123881993151,
   -0.0000000402270833384,
   -0.0000000084996745953,
   -0.0000000019514839426,
   -0.0000000004783865344,
   -0.0000000001236733992,
   -0.0000000000334137486,
   -0.0000000000093702824,
   -0.0000000000027130128,
   -0.0000000000008075954,
   -0.0000000000002463214,
   -0.0000000000000767656,
   -0.0000000000000243883,
   -0.0000000000000078831,
   -0.0000000000000025882,
   -0.0000000000000008619,
   -0.0000000000000002908,
   -0.0000000000000000993,
   -0.0000000000000000343
};
static cheb_series an21_cs = {
  an21_data,
  23,
  -1, 1,
  12
};

/*
 series for an20 on the interval -1.56250e-02 to  0.00000e+00
                                        with weighted error   4.41e-17
                                         log weighted error  16.36
                               significant figures required  14.16
                                    decimal places required  16.96
*/
static double an20_data[16] = {
    0.0126732217145738027,
   -0.0005212847072615621,
   -0.0000052672111140370,
   -0.0000001628202185026,
   -0.0000000090991442687,
   -0.0000000007438647126,
   -0.0000000000795494752,
   -0.0000000000104050944,
   -0.0000000000015932426,
   -0.0000000000002770648,
   -0.0000000000000535343,
   -0.0000000000000113062,
   -0.0000000000000025772,
   -0.0000000000000006278,
   -0.0000000000000001621,
   -0.0000000000000000441
};
static cheb_series an20_cs = {
  an20_data,
  15,
  -1, 1,
  8
};

/*
 series for aph2 on the interval -1.00000e+00 to -1.25000e-01
                                        with weighted error   2.94e-17
                                         log weighted error  16.53
                               significant figures required  15.58
                                    decimal places required  17.28
*/
static double aph2_data[32] = {
   -0.2057088719781465107,
    0.0422196961357771922,
    0.0020482560511207275,
    0.0002607800735165006,
    0.0000474824268004729,
    0.0000105102756431612,
    0.0000026353534014668,
    0.0000007208824863499,
    0.0000002103236664473,
    0.0000000644975634555,
    0.0000000205802377264,
    0.0000000067836273921,
    0.0000000022974015284,
    0.0000000007961306765,
    0.0000000002813860610,
    0.0000000001011749057,
    0.0000000000369306738,
    0.0000000000136615066,
    0.0000000000051142751,
    0.0000000000019351689,
    0.0000000000007393607,
    0.0000000000002849792,
    0.0000000000001107281,
    0.0000000000000433412,
    0.0000000000000170801,
    0.0000000000000067733,
    0.0000000000000027017,
    0.0000000000000010835,
    0.0000000000000004367,
    0.0000000000000001769,
    0.0000000000000000719,
    0.0000000000000000294
};
static cheb_series aph2_cs = {
  aph2_data,
  31,
  -1, 1,
  16
};

/*
 series for aph1 on the interval -1.25000e-01 to -1.56250e-02
                                        with weighted error   6.38e-17
                                         log weighted error  16.20
                               significant figures required  14.91
                                    decimal places required  16.87
*/
static double aph1_data[22] = {
  -0.1024172908077571694,
   0.0071697275146591248,
   0.0001209959363122329,
   0.0000073361512841220,
   0.0000007535382954272,
   0.0000001041478171741,
   0.0000000174358728519,
   0.0000000033399795033,
   0.0000000007073075174,
   0.0000000001619187515,
   0.0000000000394539982,
   0.0000000000101192282,
   0.0000000000027092778,
   0.0000000000007523806,
   0.0000000000002156369,
   0.0000000000000635283,
   0.0000000000000191757,
   0.0000000000000059143,
   0.0000000000000018597,
   0.0000000000000005950,
   0.0000000000000001934,
   0.0000000000000000638
};
static cheb_series aph1_cs = {
  aph1_data,
  21,
  -1, 1,
  10
};

/*
 series for aph0 on the interval -1.56250e-02 to  0.00000e+00
                                        with weighted error   2.29e-17
                                         log weighted error  16.64
                               significant figures required  15.27
                                    decimal places required  17.23
*/
static double aph0_data[15] = {
 -0.0855849241130933257,
  0.0011214378867065261,
  0.0000042721029353664,
  0.0000000817607381483,
  0.0000000033907645000,
  0.0000000002253264423,
  0.0000000000206284209,
  0.0000000000023858763,
  0.0000000000003301618,
  0.0000000000000527010,
  0.0000000000000094555,
  0.0000000000000018709,
  0.0000000000000004024,
  0.0000000000000000930,
  0.0000000000000000229
};
static cheb_series aph0_cs = {
  aph0_data,
  14,
  -1, 1,
  7
};

static
int
airy_deriv_mod_phase(const double x, gsl_mode_t mode,
                     gsl_sf_result * ampl, gsl_sf_result * phi)
{
  const double pi34 = 2.356194490192344928847;
  gsl_sf_result result_a;
  gsl_sf_result result_p;
  double a, p;
  double sqx;

  if(x <= -4.0) {
    double z = 128.0/(x*x*x) + 1.0;
    cheb_eval_mode_e(&an20_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph0_cs, z, mode, &result_p);
  }
  else if(x <= -2.0) {
    double z = (128.0/(x*x*x) + 9.0) / 7.0;
    cheb_eval_mode_e(&an21_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph1_cs, z, mode, &result_p);
  }
  else if(x <= -1.0) {
    double z = (16.0/(x*x*x) + 9.0) / 7.0;
    cheb_eval_mode_e(&an22_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph2_cs, z, mode, &result_p);
  }
  else {
    ampl->val = 0.0;
    ampl->err = 0.0;
    phi->val  = 0.0;
    phi->err  = 0.0;
    GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
  }

  a =  0.3125 + result_a.val;
  p = -0.625  + result_p.val;

  sqx = sqrt(-x);

  ampl->val = sqrt(a * sqx);
  ampl->err = fabs(ampl->val) * (GSL_DBL_EPSILON + fabs(result_a.err/result_a.val));
  phi->val  = pi34 - x * sqx * p;
  phi->err = fabs(phi->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));

  return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_airy_Ai_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double c    = cos(p.val);
    result->val  = a.val * c;
    result->err  = fabs(result->val * p.err) + fabs(c * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x <= 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&aif_cs, x3, mode, &result_c0);
    cheb_eval_mode_e(&aig_cs, x3, mode, &result_c1);

    result->val  = (x2*(0.125 + result_c0.val) - result_c1.val) - 0.25;
    result->err  = fabs(x2*result_c0.val) + result_c1.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
      /* scale only if x is positive */
      double s = exp(2.0*x*sqrt(x)/3.0);
      result->val *= s;
      result->err *= s;
    }

    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double sqrtx = sqrt(x);
    const double z = (16.0/(x*sqrtx) - 9.0)/7.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&aip1_cs, z, mode, &result_c0);
    result->val  = -(0.28125 + result_c0.val) * s;
    result->err  = result_c0.err * s;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sqrtx = sqrt(x);
    const double z = 16.0/(x*sqrtx) - 1.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&aip2_cs, z, mode, &result_c0);
    result->val  = -(0.28125 + result_c0.val) * s;
    result->err  = result_c0.err * s;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_airy_Ai_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double c    = cos(p.val);
    result->val  = a.val * c;
    result->err  = fabs(result->val * p.err) + fabs(c * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&aif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&aig_cs, x3, mode, &result_c2);
    result->val  = (x*x*(0.125 + result_c1.val) - result_c2.val) - 0.25;
    result->err  = fabs(x*x*result_c1.err) + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x*x*x < 9.0/4.0 * GSL_LOG_DBL_MIN*GSL_LOG_DBL_MIN) {
    gsl_sf_result result_aps;
    const double arg = -2.0*x*sqrt(x)/3.0;
    const int stat_a = gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result_aps);
    const int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
                                                result_aps.val, result_aps.err,
                                                result);
    return GSL_ERROR_SELECT_2(stat_e, stat_a);
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

int
gsl_sf_airy_Bi_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  const double atr =  8.7506905708484345;   /* 16./(sqrt(8)-1) */
  const double btr = -2.0938363213560543;   /* -(sqrt(8)+1)/(sqrt(8)-1) */

  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double s     = sin(p.val);
    result->val  = a.val * s;
    result->err  = fabs(result->val * p.err) + fabs(s * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
    result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x2 * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
      /* scale only if x is positive */
      const double s = exp(-2.0*x*sqrt(x)/3.0);
      result->val *= s;
      result->err *= s;
    }

    return GSL_SUCCESS;
  }
  else if(x < 2.0) {
    const double z = (2.0*x*x*x - 9.0) / 7.0;
    const double s = exp(-2.0*x*sqrt(x)/3.0);
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
    result->val  = s * (x*x * (0.25 + result_c0.val) + 0.5 + result_c1.val);
    result->err  = s * (x*x * result_c0.err + result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 4.0) {
    const double sqrtx = sqrt(x);
    const double z = atr/(x*sqrtx) + btr;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&bip1_cs, z, mode, &result_c0);
    result->val  = s * (0.625 + result_c0.val);
    result->err  = s * result_c0.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sqrtx = sqrt(x);
    const double z = 16.0/(x*sqrtx) - 1.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&bip2_cs, z, mode, &result_c0);
    result->val  = s * (0.625 + result_c0.val);
    result->err  = s * result_c0.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_airy_Bi_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double s    = sin(p.val);
    result->val  = a.val * s;
    result->err  = fabs(result->val * p.err) + fabs(s * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
    result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x2 * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 2.0) {
    const double z = (2.0*x*x*x - 9.0) / 7.0;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c1);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c2);
    result->val  = x*x * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x*x * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_ROOT3_DBL_MAX*GSL_ROOT3_DBL_MAX) {
    gsl_sf_result result_bps;
    const double arg = 2.0*(x*sqrt(x)/3.0);
    int stat_b = gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result_bps);
    int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
                                          result_bps.val, result_bps.err,
                                          result);
    return GSL_ERROR_SELECT_2(stat_e, stat_b);
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_airy_Ai_deriv_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result));
}

double gsl_sf_airy_Ai_deriv(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_deriv_e(x, mode, &result));
}

double gsl_sf_airy_Bi_deriv_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result));
}

double gsl_sf_airy_Bi_deriv(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_deriv_e(x, mode, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 *
 * @file core_dgbelr.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2011-05-15
 * @generated d Tue Jan  7 11:44:50 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

#define A(_m, _n)  (double *)plasma_geteltaddr(A, ((_m)-1), ((_n)-1), eltsize)
#define V(_m)      &(V[(_m)-1])
#define TAU(_m)    &(TAU[(_m)-1])

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dgbelr is a kernel that will operate on a region (triangle) of data
 *  bounded by st and ed. This kernel eliminate a column by an element-wise
 *  annihiliation, and for each annihiliation, it apply a left update, followed
 *  by an right update on the diagonal 2x2 element, where it create a new nnz,
 *  eliminate it and update V and TAU, then it continue until finishing the the
 *  whole column. When this is done, it take advantage that data are on cache
 *  and will apply the right on the remaining part of this region that has not
 *  been updated by the right yet.
 *  Note that the column to be eliminated is located at st-1.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *         @arg PlasmaLower:
 *         @arg PlasmaUpper:
 *
 * @param[in] N
 *          The order of the matrix A.
 *
 * @param[in, out] A
 *          A pointer to the descriptor of the matrix A.
 *
 * @param[out] V
 *          double array, dimension (N).
 *          The scalar elementary reflectors are written in this
 *          array. So it is used as a workspace for V at each step
 *          of the bulge chasing algorithm.
 *
 * @param[out] TAU
 *          double array, dimension (N).
 *          The scalar factors of the elementary reflectors are written
 *          in thisarray. So it is used as a workspace for TAU at each step
 *          of the bulge chasing algorithm.
 *
 * @param[in] st
 *          A pointer to the start index where this kernel will operate.
 *
 * @param[in] ed
 *          A pointer to the end index where this kernel will operate.
 *
 * @param[in] eltsize
 *          PLASMA internal value which refer to the size of the precision.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
/***************************************************************************//**
 *  TYPE 1-BDL Householder
 *  add -1 because of C
 ******************************************************************************/
int
CORE_dgbelr(PLASMA_enum uplo, int N,
            PLASMA_desc *A,
            double *V,
            double *TAU,
            int st,
            int ed,
            int eltsize)
{
    int    NB, J1, J2;
    int    len1, len2, t1ed, t2st;
    int    i;
    static double zzero = 0.0;
    PLASMA_desc vA=*A;

    /* Check input arguments */
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (ed <= st) {
        coreblas_error(6, "Illegal value of st and ed (internal)");
        return -6;
    }

    /* Quick return */
    if (N == 0)
        return PLASMA_SUCCESS;

    NB = A->mb;
    if( uplo == PlasmaLower ){
        /* ========================
         *       LOWER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* generate Householder to annihilate a(i+k-1,i) within the band*/
            *V(i)          = *A(i, (st-1));
            *A(i, (st-1))  = zzero;
            LAPACKE_dlarfg_work( 2, A((i-1),(st-1)), V(i), 1, TAU(i));

            /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
            J1    = st;
            J2    = i-2;
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1;
            len2  = J2-t2st+1;
            if(len1>0)CORE_dlarfx2(PlasmaLeft, len1 , *V(i), (*TAU(i)), A(i-1, J1  ), ELTLDD(vA, (i-1)),  A(i,  J1 ), ELTLDD(vA, i) );
            if(len2>0)CORE_dlarfx2(PlasmaLeft, len2 , *V(i), (*TAU(i)), A(i-1, t2st), ELTLDD(vA, (i-1)),  A(i, t2st), ELTLDD(vA, i) );
            CORE_dlarfx2ce(PlasmaLower, V(i), TAU(i), A(i-1,i-1), A(i,i-1), A(i,i));
        }
        /* APPLY RIGHT ON THE REMAINING ELEMENT OF KERNEL 1 */
        for (i = ed; i >= st+1 ; i--){
            J1    = i+1;
            J2    = min(ed,N);
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1;
            len2  = J2-t2st+1;
            if(len1>0)CORE_dlarfx2(PlasmaRight, len1, (*V(i)), (*TAU(i)), A(J1,i-1),   ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
            if(len2>0)CORE_dlarfx2(PlasmaRight, len2, (*V(i)), (*TAU(i)), A(t2st,i-1), ELTLDD(vA, t2st), A(t2st, i), ELTLDD(vA, t2st) );
        }
    } else {
        /* ========================
         *       UPPER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* generate Householder to annihilate a(i+k-1,i) within the band*/
            *V(i)          = *A((st-1),  i);
            *A((st-1),  i) = zzero;
            LAPACKE_dlarfg_work( 2, A((st-1), (i-1)), V(i), 1, TAU(i));

            /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
            J1    = st;
            J2    = i-2;
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1;
            len2  = J2-t2st+1;
            if(len1>0)CORE_dlarfx2(PlasmaRight, len1, (*V(i)), (*TAU(i)), A(J1,i-1),   ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
            if(len2>0)CORE_dlarfx2(PlasmaRight, len2, (*V(i)), (*TAU(i)), A(t2st,i-1), ELTLDD(vA, t2st), A(t2st, i), ELTLDD(vA, t2st) );
            CORE_dlarfx2ce(PlasmaUpper, V(i), TAU(i), A((i-1),(i-1)), A((i-1), i), A(i,i));
        }
        /* APPLY LEFT ON THE REMAINING ELEMENT OF KERNEL 1*/
        for (i = ed; i >= st+1 ; i--){
            J1    = i+1;
            J2    = min(ed,N);
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1;
            len2  = J2-t2st+1;
            if(len1>0)CORE_dlarfx2(PlasmaLeft, len1 , *V(i), (*TAU(i)), A(i-1, J1  ), ELTLDD(vA, (i-1)),  A(i,  J1 ), ELTLDD(vA, i) );
            if(len2>0)CORE_dlarfx2(PlasmaLeft, len2 , *V(i), (*TAU(i)), A(i-1, t2st), ELTLDD(vA, (i-1)),  A(i, t2st), ELTLDD(vA, i) );
        }
    }  /* end of else for the upper case*/

    return PLASMA_SUCCESS;

}


text:
#ifndef Photometry_h
#define Photometry_h
/**
 * @file
 * $Revision: 1.5 $
 * $Date: 2008/07/09 19:40:52 $
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are
 *   public domain. See individual third-party library and package descriptions
 *   for intellectual property information, user agreements, and related
 *   information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or
 *   implied, is made by the USGS as to the accuracy and functioning of such
 *   software and related material nor shall the fact of distribution
 *   constitute any such warranty, and no responsibility is assumed by the
 *   USGS in connection therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
 *   in a browser or see the Privacy & Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <string>
#include <vector>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_roots.h>

namespace Isis {
  class Pvl;
  class PhotoModel;
  class AtmosModel;
  class NormModel;
  /**
   * @author ????-??-?? Unknown
   *
   * @internal
   *  @history 2007-08-02 Steven Lambright - Fixed memory leak
   *  @history 2008-03-07 Janet Barrett - Added SetPhotomWl method to allow
   *                      the application to set the p_normWavelength variable for
   *                      use by MoonAlbedo normalization.
   *  @history 2008-06-18 Steven Koechle - Fixed Documentation Errors
   *  @history 2008-07-09 Steven Lambright - Fixed unit test
   *  @history 2011-08-19 Sharmila Prasad - Implemented brentminimizer using GSL
   *  @history 2011-09-15 Sharmila Prasad - Implemented brent's root solver using GSL
   */
  class Photometry {
    public:
      Photometry(Pvl &pvl);
      Photometry() {};
      virtual ~Photometry();

      //! Calculate the surface brightness
      void Compute(double pha, double inc, double ema, double dn,
                   double &albedo, double &mult, double &base);
      void Compute(double pha, double inc, double ema, double deminc,
                   double demema, double dn, double &albedo,
                   double &mult, double &base);

      //! Set the wavelength
      virtual void SetPhotomWl(double wl);

      //! Double precision version of bracketing algorithm ported from Python.
      //! Solution bracketing for 1-D minimization routine.
      static void minbracket(double &xa, double &xb, double &xc, double &fa,
          double &fb, double &fc, double Func(double par, void *params),
          void *params);

      //! Brent's method 1-D minimization routine using GSL's r8Brent minimization Algorithm
      static int brentminimizer(double x_lower, double x_upper, gsl_function *Func,
          double & x_minimum, double tolerance);

      //! GSL's the Brent-Dekker method (Brent's method) combines an interpolation strategy
      //! with the bisection algorithm to estimate the root of the quadratic function.
      static int brentsolver(double x_lo, double x_hi, gsl_function *Func, double tolerance, double &root);

      PhotoModel *GetPhotoModel() const {
        return p_phtPmodel;
      }

      AtmosModel *GetAtmosModel() const {
        return p_phtAmodel;
      }

      NormModel *GetNormModel() const {
        return p_phtNmodel;
      }

    protected:
      AtmosModel *p_phtAmodel;
      PhotoModel *p_phtPmodel;
      NormModel *p_phtNmodel;
  };
};

#endif


text:
#ifndef TTT_BOARD_H
#define TTT_BOARD_H

#include <SDL_rect.h>

#include <gsl/util>

#include <array>
#include <optional>

namespace ttt
{
	class Game;

	enum class CellState
	{
		empty, x, o
	};

	class Board final
	{
	public:
		SDL_FPoint pos;
		SDL_FPoint cellSize;
		static constexpr int size{3};

		CellState preview{CellState::empty};
		SDL_Point previewCell{};

		Board(const SDL_FPoint &pos, const SDL_FPoint &cellSize, const Game &game) noexcept;

		void update() noexcept;
		void render() const noexcept;

		[[nodiscard]] CellState getCellState(const SDL_Point &pos) const noexcept;
		void setCellState(const SDL_Point &pos, CellState state) noexcept;
		void clear() noexcept;

		[[nodiscard]] std::optional<SDL_Point> getCell(const SDL_Point &point) const noexcept;
		[[nodiscard]] bool isEmpty() const noexcept;

	private:
		const Game &game;
		std::array<CellState, 9> grid{};

		static constexpr float lineWidth{16.0f};
		float previewOffset{0.0f};
		bool shrinkPreview{false};

		const SDL_Rect lineSrc{0, 0, 8, 576};
		const SDL_Rect xSrc{8, 0, 48, 48};
		const SDL_Rect oSrc{8, 48, 48, 48};
	};

	inline CellState Board::getCellState(const SDL_Point &pos) const noexcept
	{
		return gsl::at(grid, pos.x + pos.y * size);
	}

	inline void Board::setCellState(const SDL_Point &pos, CellState state) noexcept
	{
		gsl::at(grid, pos.x + pos.y * size) = state;
	}
}

#endif

max_stars_repo_name: itsArtem/TicTacToe

max_issues_repo_name: itsArtem/TicTacToe

max_forks_repo_name: itsArtem/TicTacToe


text:
#ifndef TTT_UTIL_FONTCACHE_H
#define TTT_UTIL_FONTCACHE_H

#include <SDL_ttf.h>

#include <gsl/pointers>
#include <gsl/util>

#include <vector>
#include <memory>
#include <string_view>

namespace ttt::gfx
{
	// A class for caching fonts for reuse.
	class FontCache final
	{
	public:
		using SizeType = gsl::index;

		TTF_Font  *load(std::string_view path, int size, int index = 0);
		void unload(SizeType index) noexcept;
		void erase(SizeType index) noexcept;
		void clear() noexcept;

		[[nodiscard]] TTF_Font *at(SizeType index) const noexcept;
		[[nodiscard]] TTF_Font *operator [](SizeType index) const noexcept;

		[[nodiscard]] SizeType getSize() const noexcept;
		[[nodiscard]] SizeType getMaxSize() const noexcept;

	private:
		std::vector<std::unique_ptr<TTF_Font, decltype(&TTF_CloseFont)>> fonts;
	};

	inline void FontCache::unload(SizeType index) noexcept { gsl::at(fonts, index).reset(); }
	inline void FontCache::erase(SizeType index) noexcept { fonts.erase(fonts.cbegin() + index); }
	inline void FontCache::clear() noexcept { fonts.clear(); }

	inline TTF_Font *FontCache::at(SizeType index) const noexcept { return gsl::at(fonts, index).get(); }
	inline TTF_Font *FontCache::operator [](SizeType index) const noexcept { return gsl::at(fonts, index).get(); }

	inline FontCache::SizeType FontCache::getSize() const noexcept { return fonts.size(); }
	inline FontCache::SizeType FontCache::getMaxSize() const noexcept { return fonts.max_size(); }
}

#endif


text:
/**
 *
 * @file core_claset.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_claset - Sets the elements of the matrix A on the diagonal
 *  to beta and on the off-diagonals to alpha
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which elements of the matrix are to be set
 *          = PlasmaUpper: Upper part of A is set;
 *          = PlasmaLower: Lower part of A is set;
 *          = PlasmaUpperLower: ALL elements of A are set.
 *
 * @param[in] M
 *          The number of rows of the matrix A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the matrix A.  N >= 0.
 *
 * @param[in] alpha
 *         The constant to which the off-diagonal elements are to be set.
 *
 * @param[in] beta
 *         The constant to which the diagonal elements are to be set.
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile A.
 *         On exit, A has been set accordingly.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claset = PCORE_claset
#define CORE_claset PCORE_claset
#endif
void CORE_claset(PLASMA_enum uplo, int M, int N,
                 PLASMA_Complex32_t alpha, PLASMA_Complex32_t beta,
                 PLASMA_Complex32_t *A, int LDA)
{
    LAPACKE_claset_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, alpha, beta, A, LDA);
}