/***************************************************************************//**
**  \mainpage Mesh-based Monte Carlo (MMC) - a 3D photon simulator
**
**  \author Qianqian Fang <q.fang at neu.edu>
**  \copyright Qianqian Fang, 2010-2024
**
**  \section sref Reference:
**  \li \c (\b Fang2010) Qianqian Fang, <a href="http://www.opticsinfobase.org/abstract.cfm?uri=boe-1-1-165">
**          "Mesh-based Monte Carlo Method Using Fast Ray-Tracing
**          in Plucker Coordinates,"</a> Biomed. Opt. Express, 1(1) 165-175 (2010).
**  \li \c (\b Fang2012) Qianqian Fang and David R. Kaeli,
**           <a href="https://www.osapublishing.org/boe/abstract.cfm?uri=boe-3-12-3223">
**          "Accelerating mesh-based Monte Carlo method on modern CPU architectures,"</a>
**          Biomed. Opt. Express 3(12), 3223-3230 (2012)
**  \li \c (\b Yao2016) Ruoyang Yao, Xavier Intes, and Qianqian Fang,
**          <a href="https://www.osapublishing.org/boe/abstract.cfm?uri=boe-7-1-171">
**          "Generalized mesh-based Monte Carlo for wide-field illumination and detection
**           via mesh retessellation,"</a> Biomed. Optics Express, 7(1), 171-184 (2016)
**  \li \c (\b Fang2019) Qianqian Fang and Shijie Yan,
**          <a href="http://dx.doi.org/10.1117/1.JBO.24.11.115002">
**          "Graphics processing unit-accelerated mesh-based Monte Carlo photon transport
**           simulations,"</a> J. of Biomedical Optics, 24(11), 115002 (2019)
**  \li \c (\b Yuan2021) Yaoshen Yuan, Shijie Yan, and Qianqian Fang,
**          <a href="https://www.osapublishing.org/boe/fulltext.cfm?uri=boe-12-1-147">
**          "Light transport modeling in highly complex tissues using the implicit
**           mesh-based Monte Carlo algorithm,"</a> Biomed. Optics Express, 12(1) 147-161 (2021)
**
**  \section slicense License
**          GPL v3, see LICENSE.txt for details
*******************************************************************************/

/**
   Copyright (C) 1995, 1997, 2001 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, August 1995.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library 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
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.
*/

/***************************************************************************//**
\file    drand48_r_libgw32c.c

\brief   POSIX 48bit multi-threaded RNG for win32 port via LibGW32C
*******************************************************************************/

#include "mmc_rand_drand48.h"

struct drand48_data __libc_drand48_data;

int __drand48_iterate (unsigned short int xsubi[3], struct drand48_data* buffer ) {
    uint64 X;
    uint64 result;

    /* Initialize buffer, if not yet done. */

    // built in bullshit removed
    // if (__builtin_expect (!buffer->__init, 0))
    if (buffer->__init == 0) {
        buffer->__a = 0x5deece66dull;
        buffer->__c = 0xb;
        buffer->__init = 1;
    }

    /* Do the real work. We choose a data type which contains at least
    48 bits. Because we compute the modulus it does not care how
    many bits really are computed. */

    X = (uint64) xsubi[2] << 32 | (uint32) xsubi[1] << 16 | xsubi[0];

    result = X * buffer->__a + buffer->__c;

    xsubi[0] = result & 0xffff;
    xsubi[1] = (result >> 16) & 0xffff;
    xsubi[2] = (result >> 32) & 0xffff;

    return 0;
}

int __erand48_r (unsigned short int xsubi[3], struct drand48_data* buffer, double* result) {
    union ieee754_double temp;

    /* Compute next state.  */
    if (__drand48_iterate (xsubi, buffer) < 0) {
        return -1;
    }

    /* Construct a positive double with the 48 random bits distributed over
       its fractional part so the resulting FP number is [0.0,1.0).  */

    temp.ieee.negative = 0;
    temp.ieee.exponent = IEEE754_DOUBLE_BIAS;
    temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12);
    temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4);

    /* Please note the lower 4 bits of mantissa1 are always 0.  */
    *result = temp.d - 1.0;

    return 0;
}
//weak_alias (__erand48_r, erand48_r)

int drand48_r (struct drand48_data* buffer, double* result) {
    return __erand48_r (buffer->__x, buffer, result);
}

int seed48_r (unsigned short int seed16v[3], struct drand48_data* buffer) {
    /* Save old value at a private place to be used as return value.  */
    memcpy (buffer->__old_x, buffer->__x, sizeof (buffer->__x));

    /* Install new state.  */
    buffer->__x[2] = seed16v[2];
    buffer->__x[1] = seed16v[1];
    buffer->__x[0] = seed16v[0];
    buffer->__a = 0x5deece66dull;
    buffer->__c = 0xb;
    buffer->__init = 1;

    return 0;
}

double erand48 (unsigned short int xsubi[3]) {
    double result;
    (void) __erand48_r (xsubi, &__libc_drand48_data, &result);
    return result;
}

//weak_alias (__seed48_r, seed48_r)
/*
int main()
{
  struct drand48_data st;
  double t;
  drand48_r(&st,&t);
  printf("%f \n", t );
  return 0;
}
*/
