{-
   Copyright 2016, Dominic Orchard, Andrew Rice, Mistral Contrastin, Matthew Danish

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
-}

{-# LANGUAGE ForeignFunctionInterface   #-}
{-# LANGUAGE TupleSections              #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

module Camfort.Specification.Units.InferenceBackendFlint where

import           Control.Monad
import           Data.List (partition)
-- import           Debug.Trace (trace, traceM, traceShowM)
import           Foreign
import           Foreign.C.Types
import           Numeric.LinearAlgebra (atIndex, rows, cols)
import qualified Numeric.LinearAlgebra as H
import           System.IO.Unsafe (unsafePerformIO)

-- | Units of measure extension to Fortran: Flint backend components
--
-- Some notes on the Flint library to aid comprehension of the original C and
-- this interface:
--
--   * They use a @typedef TYPE TYPE_t[1]@ convention to do call-by-reference
--     without an explicit pointer. It appears to be an unmentioned convention
--     borrowed from related library & depedency GMP, explained in a GMP doc
--     page: https://gmplib.org/manual/Parameter-Conventions . Any time one of
--     these is a function parameter, it is correct to use 'Ptr a' in Haskell.
--   * Flint extensively uses two typedefs @ulong@ and @slong@, which are "long
--     integers" in unsigned and signed representations respectively. However,
--     the story is more complicated in cross-platform contexts, because 64-bit
--     Linux's @long@s are 64 bits (8 bytes), while 64-bit Windows kept them at
--     32 bits (4 bytes). That type is 'CLong' in Haskell, and it doesn't match
--     up with Flint's @slong@, so we roll our own newtypes instead. (See the
--     definition for further explanation.)

-- | @typedef slong fmpz@
--
-- GHC's generalized newtype deriving handles deriving all the instances we
-- require for us.
newtype FMPZ = FMPZ { FMPZ -> SLong
unFMPZ :: SLong }
  deriving (Ptr b -> Int -> IO FMPZ
Ptr b -> Int -> FMPZ -> IO ()
Ptr FMPZ -> IO FMPZ
Ptr FMPZ -> Int -> IO FMPZ
Ptr FMPZ -> Int -> FMPZ -> IO ()
Ptr FMPZ -> FMPZ -> IO ()
FMPZ -> Int
(FMPZ -> Int)
-> (FMPZ -> Int)
-> (Ptr FMPZ -> Int -> IO FMPZ)
-> (Ptr FMPZ -> Int -> FMPZ -> IO ())
-> (forall b. Ptr b -> Int -> IO FMPZ)
-> (forall b. Ptr b -> Int -> FMPZ -> IO ())
-> (Ptr FMPZ -> IO FMPZ)
-> (Ptr FMPZ -> FMPZ -> IO ())
-> Storable FMPZ
forall b. Ptr b -> Int -> IO FMPZ
forall b. Ptr b -> Int -> FMPZ -> IO ()
forall a.
(a -> Int)
-> (a -> Int)
-> (Ptr a -> Int -> IO a)
-> (Ptr a -> Int -> a -> IO ())
-> (forall b. Ptr b -> Int -> IO a)
-> (forall b. Ptr b -> Int -> a -> IO ())
-> (Ptr a -> IO a)
-> (Ptr a -> a -> IO ())
-> Storable a
poke :: Ptr FMPZ -> FMPZ -> IO ()
$cpoke :: Ptr FMPZ -> FMPZ -> IO ()
peek :: Ptr FMPZ -> IO FMPZ
$cpeek :: Ptr FMPZ -> IO FMPZ
pokeByteOff :: Ptr b -> Int -> FMPZ -> IO ()
$cpokeByteOff :: forall b. Ptr b -> Int -> FMPZ -> IO ()
peekByteOff :: Ptr b -> Int -> IO FMPZ
$cpeekByteOff :: forall b. Ptr b -> Int -> IO FMPZ
pokeElemOff :: Ptr FMPZ -> Int -> FMPZ -> IO ()
$cpokeElemOff :: Ptr FMPZ -> Int -> FMPZ -> IO ()
peekElemOff :: Ptr FMPZ -> Int -> IO FMPZ
$cpeekElemOff :: Ptr FMPZ -> Int -> IO FMPZ
alignment :: FMPZ -> Int
$calignment :: FMPZ -> Int
sizeOf :: FMPZ -> Int
$csizeOf :: FMPZ -> Int
Storable, FMPZ -> FMPZ -> Bool
(FMPZ -> FMPZ -> Bool) -> (FMPZ -> FMPZ -> Bool) -> Eq FMPZ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: FMPZ -> FMPZ -> Bool
$c/= :: FMPZ -> FMPZ -> Bool
== :: FMPZ -> FMPZ -> Bool
$c== :: FMPZ -> FMPZ -> Bool
Eq, Eq FMPZ
Eq FMPZ
-> (FMPZ -> FMPZ -> Ordering)
-> (FMPZ -> FMPZ -> Bool)
-> (FMPZ -> FMPZ -> Bool)
-> (FMPZ -> FMPZ -> Bool)
-> (FMPZ -> FMPZ -> Bool)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> Ord FMPZ
FMPZ -> FMPZ -> Bool
FMPZ -> FMPZ -> Ordering
FMPZ -> FMPZ -> FMPZ
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: FMPZ -> FMPZ -> FMPZ
$cmin :: FMPZ -> FMPZ -> FMPZ
max :: FMPZ -> FMPZ -> FMPZ
$cmax :: FMPZ -> FMPZ -> FMPZ
>= :: FMPZ -> FMPZ -> Bool
$c>= :: FMPZ -> FMPZ -> Bool
> :: FMPZ -> FMPZ -> Bool
$c> :: FMPZ -> FMPZ -> Bool
<= :: FMPZ -> FMPZ -> Bool
$c<= :: FMPZ -> FMPZ -> Bool
< :: FMPZ -> FMPZ -> Bool
$c< :: FMPZ -> FMPZ -> Bool
compare :: FMPZ -> FMPZ -> Ordering
$ccompare :: FMPZ -> FMPZ -> Ordering
$cp1Ord :: Eq FMPZ
Ord, Integer -> FMPZ
FMPZ -> FMPZ
FMPZ -> FMPZ -> FMPZ
(FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ)
-> (FMPZ -> FMPZ)
-> (FMPZ -> FMPZ)
-> (Integer -> FMPZ)
-> Num FMPZ
forall a.
(a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (Integer -> a)
-> Num a
fromInteger :: Integer -> FMPZ
$cfromInteger :: Integer -> FMPZ
signum :: FMPZ -> FMPZ
$csignum :: FMPZ -> FMPZ
abs :: FMPZ -> FMPZ
$cabs :: FMPZ -> FMPZ
negate :: FMPZ -> FMPZ
$cnegate :: FMPZ -> FMPZ
* :: FMPZ -> FMPZ -> FMPZ
$c* :: FMPZ -> FMPZ -> FMPZ
- :: FMPZ -> FMPZ -> FMPZ
$c- :: FMPZ -> FMPZ -> FMPZ
+ :: FMPZ -> FMPZ -> FMPZ
$c+ :: FMPZ -> FMPZ -> FMPZ
Num, Num FMPZ
Ord FMPZ
Num FMPZ -> Ord FMPZ -> (FMPZ -> Rational) -> Real FMPZ
FMPZ -> Rational
forall a. Num a -> Ord a -> (a -> Rational) -> Real a
toRational :: FMPZ -> Rational
$ctoRational :: FMPZ -> Rational
$cp2Real :: Ord FMPZ
$cp1Real :: Num FMPZ
Real, Int -> FMPZ
FMPZ -> Int
FMPZ -> [FMPZ]
FMPZ -> FMPZ
FMPZ -> FMPZ -> [FMPZ]
FMPZ -> FMPZ -> FMPZ -> [FMPZ]
(FMPZ -> FMPZ)
-> (FMPZ -> FMPZ)
-> (Int -> FMPZ)
-> (FMPZ -> Int)
-> (FMPZ -> [FMPZ])
-> (FMPZ -> FMPZ -> [FMPZ])
-> (FMPZ -> FMPZ -> [FMPZ])
-> (FMPZ -> FMPZ -> FMPZ -> [FMPZ])
-> Enum FMPZ
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: FMPZ -> FMPZ -> FMPZ -> [FMPZ]
$cenumFromThenTo :: FMPZ -> FMPZ -> FMPZ -> [FMPZ]
enumFromTo :: FMPZ -> FMPZ -> [FMPZ]
$cenumFromTo :: FMPZ -> FMPZ -> [FMPZ]
enumFromThen :: FMPZ -> FMPZ -> [FMPZ]
$cenumFromThen :: FMPZ -> FMPZ -> [FMPZ]
enumFrom :: FMPZ -> [FMPZ]
$cenumFrom :: FMPZ -> [FMPZ]
fromEnum :: FMPZ -> Int
$cfromEnum :: FMPZ -> Int
toEnum :: Int -> FMPZ
$ctoEnum :: Int -> FMPZ
pred :: FMPZ -> FMPZ
$cpred :: FMPZ -> FMPZ
succ :: FMPZ -> FMPZ
$csucc :: FMPZ -> FMPZ
Enum, Enum FMPZ
Real FMPZ
Real FMPZ
-> Enum FMPZ
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> FMPZ)
-> (FMPZ -> FMPZ -> (FMPZ, FMPZ))
-> (FMPZ -> FMPZ -> (FMPZ, FMPZ))
-> (FMPZ -> Integer)
-> Integral FMPZ
FMPZ -> Integer
FMPZ -> FMPZ -> (FMPZ, FMPZ)
FMPZ -> FMPZ -> FMPZ
forall a.
Real a
-> Enum a
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> (a, a))
-> (a -> a -> (a, a))
-> (a -> Integer)
-> Integral a
toInteger :: FMPZ -> Integer
$ctoInteger :: FMPZ -> Integer
divMod :: FMPZ -> FMPZ -> (FMPZ, FMPZ)
$cdivMod :: FMPZ -> FMPZ -> (FMPZ, FMPZ)
quotRem :: FMPZ -> FMPZ -> (FMPZ, FMPZ)
$cquotRem :: FMPZ -> FMPZ -> (FMPZ, FMPZ)
mod :: FMPZ -> FMPZ -> FMPZ
$cmod :: FMPZ -> FMPZ -> FMPZ
div :: FMPZ -> FMPZ -> FMPZ
$cdiv :: FMPZ -> FMPZ -> FMPZ
rem :: FMPZ -> FMPZ -> FMPZ
$crem :: FMPZ -> FMPZ -> FMPZ
quot :: FMPZ -> FMPZ -> FMPZ
$cquot :: FMPZ -> FMPZ -> FMPZ
$cp2Integral :: Enum FMPZ
$cp1Integral :: Real FMPZ
Integral)

-- | Flint's long signed integer type @slong@ (= GMP's @mp_limb_signed_t@).
--
-- As described in their Portability doc page
-- https://flintlib.org/doc/portability.html , this replaces @long@ (long signed
-- integer). Importantly, it is *always* 64-bits, regardless of platform. @long@
-- on Windows is usually 32-bits (whether on a 32-bit or 64-bit install), and
-- you're meant to use @long long@ instead.
--
-- We tie the typedef to Haskell's 'Int64', since that should be the appropriate
-- size for any regular platform. Better would be to do some CPP or hsc2hs magic
-- to check the size of an @slong@ and use the appropriate Haskell signed
-- integer type.
--
-- GHC's generalized newtype deriving handles deriving all the instances we
-- require for us.
newtype SLong = SLong { SLong -> Int64
unSLong :: Int64 }
  deriving (Ptr b -> Int -> IO SLong
Ptr b -> Int -> SLong -> IO ()
Ptr SLong -> IO SLong
Ptr SLong -> Int -> IO SLong
Ptr SLong -> Int -> SLong -> IO ()
Ptr SLong -> SLong -> IO ()
SLong -> Int
(SLong -> Int)
-> (SLong -> Int)
-> (Ptr SLong -> Int -> IO SLong)
-> (Ptr SLong -> Int -> SLong -> IO ())
-> (forall b. Ptr b -> Int -> IO SLong)
-> (forall b. Ptr b -> Int -> SLong -> IO ())
-> (Ptr SLong -> IO SLong)
-> (Ptr SLong -> SLong -> IO ())
-> Storable SLong
forall b. Ptr b -> Int -> IO SLong
forall b. Ptr b -> Int -> SLong -> IO ()
forall a.
(a -> Int)
-> (a -> Int)
-> (Ptr a -> Int -> IO a)
-> (Ptr a -> Int -> a -> IO ())
-> (forall b. Ptr b -> Int -> IO a)
-> (forall b. Ptr b -> Int -> a -> IO ())
-> (Ptr a -> IO a)
-> (Ptr a -> a -> IO ())
-> Storable a
poke :: Ptr SLong -> SLong -> IO ()
$cpoke :: Ptr SLong -> SLong -> IO ()
peek :: Ptr SLong -> IO SLong
$cpeek :: Ptr SLong -> IO SLong
pokeByteOff :: Ptr b -> Int -> SLong -> IO ()
$cpokeByteOff :: forall b. Ptr b -> Int -> SLong -> IO ()
peekByteOff :: Ptr b -> Int -> IO SLong
$cpeekByteOff :: forall b. Ptr b -> Int -> IO SLong
pokeElemOff :: Ptr SLong -> Int -> SLong -> IO ()
$cpokeElemOff :: Ptr SLong -> Int -> SLong -> IO ()
peekElemOff :: Ptr SLong -> Int -> IO SLong
$cpeekElemOff :: Ptr SLong -> Int -> IO SLong
alignment :: SLong -> Int
$calignment :: SLong -> Int
sizeOf :: SLong -> Int
$csizeOf :: SLong -> Int
Storable, SLong -> SLong -> Bool
(SLong -> SLong -> Bool) -> (SLong -> SLong -> Bool) -> Eq SLong
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: SLong -> SLong -> Bool
$c/= :: SLong -> SLong -> Bool
== :: SLong -> SLong -> Bool
$c== :: SLong -> SLong -> Bool
Eq, Eq SLong
Eq SLong
-> (SLong -> SLong -> Ordering)
-> (SLong -> SLong -> Bool)
-> (SLong -> SLong -> Bool)
-> (SLong -> SLong -> Bool)
-> (SLong -> SLong -> Bool)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> Ord SLong
SLong -> SLong -> Bool
SLong -> SLong -> Ordering
SLong -> SLong -> SLong
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: SLong -> SLong -> SLong
$cmin :: SLong -> SLong -> SLong
max :: SLong -> SLong -> SLong
$cmax :: SLong -> SLong -> SLong
>= :: SLong -> SLong -> Bool
$c>= :: SLong -> SLong -> Bool
> :: SLong -> SLong -> Bool
$c> :: SLong -> SLong -> Bool
<= :: SLong -> SLong -> Bool
$c<= :: SLong -> SLong -> Bool
< :: SLong -> SLong -> Bool
$c< :: SLong -> SLong -> Bool
compare :: SLong -> SLong -> Ordering
$ccompare :: SLong -> SLong -> Ordering
$cp1Ord :: Eq SLong
Ord, Integer -> SLong
SLong -> SLong
SLong -> SLong -> SLong
(SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong)
-> (SLong -> SLong)
-> (SLong -> SLong)
-> (Integer -> SLong)
-> Num SLong
forall a.
(a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (Integer -> a)
-> Num a
fromInteger :: Integer -> SLong
$cfromInteger :: Integer -> SLong
signum :: SLong -> SLong
$csignum :: SLong -> SLong
abs :: SLong -> SLong
$cabs :: SLong -> SLong
negate :: SLong -> SLong
$cnegate :: SLong -> SLong
* :: SLong -> SLong -> SLong
$c* :: SLong -> SLong -> SLong
- :: SLong -> SLong -> SLong
$c- :: SLong -> SLong -> SLong
+ :: SLong -> SLong -> SLong
$c+ :: SLong -> SLong -> SLong
Num, Num SLong
Ord SLong
Num SLong -> Ord SLong -> (SLong -> Rational) -> Real SLong
SLong -> Rational
forall a. Num a -> Ord a -> (a -> Rational) -> Real a
toRational :: SLong -> Rational
$ctoRational :: SLong -> Rational
$cp2Real :: Ord SLong
$cp1Real :: Num SLong
Real, Int -> SLong
SLong -> Int
SLong -> [SLong]
SLong -> SLong
SLong -> SLong -> [SLong]
SLong -> SLong -> SLong -> [SLong]
(SLong -> SLong)
-> (SLong -> SLong)
-> (Int -> SLong)
-> (SLong -> Int)
-> (SLong -> [SLong])
-> (SLong -> SLong -> [SLong])
-> (SLong -> SLong -> [SLong])
-> (SLong -> SLong -> SLong -> [SLong])
-> Enum SLong
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: SLong -> SLong -> SLong -> [SLong]
$cenumFromThenTo :: SLong -> SLong -> SLong -> [SLong]
enumFromTo :: SLong -> SLong -> [SLong]
$cenumFromTo :: SLong -> SLong -> [SLong]
enumFromThen :: SLong -> SLong -> [SLong]
$cenumFromThen :: SLong -> SLong -> [SLong]
enumFrom :: SLong -> [SLong]
$cenumFrom :: SLong -> [SLong]
fromEnum :: SLong -> Int
$cfromEnum :: SLong -> Int
toEnum :: Int -> SLong
$ctoEnum :: Int -> SLong
pred :: SLong -> SLong
$cpred :: SLong -> SLong
succ :: SLong -> SLong
$csucc :: SLong -> SLong
Enum, Enum SLong
Real SLong
Real SLong
-> Enum SLong
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> SLong)
-> (SLong -> SLong -> (SLong, SLong))
-> (SLong -> SLong -> (SLong, SLong))
-> (SLong -> Integer)
-> Integral SLong
SLong -> Integer
SLong -> SLong -> (SLong, SLong)
SLong -> SLong -> SLong
forall a.
Real a
-> Enum a
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> (a, a))
-> (a -> a -> (a, a))
-> (a -> Integer)
-> Integral a
toInteger :: SLong -> Integer
$ctoInteger :: SLong -> Integer
divMod :: SLong -> SLong -> (SLong, SLong)
$cdivMod :: SLong -> SLong -> (SLong, SLong)
quotRem :: SLong -> SLong -> (SLong, SLong)
$cquotRem :: SLong -> SLong -> (SLong, SLong)
mod :: SLong -> SLong -> SLong
$cmod :: SLong -> SLong -> SLong
div :: SLong -> SLong -> SLong
$cdiv :: SLong -> SLong -> SLong
rem :: SLong -> SLong -> SLong
$crem :: SLong -> SLong -> SLong
quot :: SLong -> SLong -> SLong
$cquot :: SLong -> SLong -> SLong
$cp2Integral :: Enum SLong
$cp1Integral :: Real SLong
Integral)

{-
    typedef struct
    {
        fmpz * entries;
        slong r;
        slong c;
        fmpz ** rows;
    } fmpz_mat_struct;
-}
data FMPZMat

instance Storable FMPZMat where
  sizeOf :: FMPZMat -> Int
sizeOf FMPZMat
_ =
        Ptr Any -> Int
forall a. Storable a => a -> Int
sizeOf Ptr Any
forall a. Ptr a
nullPtr
      Int -> Int -> Int
forall a. Num a => a -> a -> a
+ SLong -> Int
forall a. Storable a => a -> Int
sizeOf (SLong
forall a. HasCallStack => a
undefined :: SLong)
      Int -> Int -> Int
forall a. Num a => a -> a -> a
+ SLong -> Int
forall a. Storable a => a -> Int
sizeOf (SLong
forall a. HasCallStack => a
undefined :: SLong)
      Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Ptr Any -> Int
forall a. Storable a => a -> Int
sizeOf Ptr Any
forall a. Ptr a
nullPtr
  alignment :: FMPZMat -> Int
alignment FMPZMat
_ = SLong -> Int
forall a. Storable a => a -> Int
alignment (SLong
forall a. HasCallStack => a
undefined :: SLong) -- TODO: ??
  peek :: Ptr FMPZMat -> IO FMPZMat
peek Ptr FMPZMat
_ = IO FMPZMat
forall a. HasCallStack => a
undefined
  poke :: Ptr FMPZMat -> FMPZMat -> IO ()
poke Ptr FMPZMat
_ = FMPZMat -> IO ()
forall a. HasCallStack => a
undefined

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_init" fmpz_mat_init :: Ptr FMPZMat -> SLong -> SLong -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_set" fmpz_mat_set :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_entry" fmpz_mat_entry :: Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
foreign import ccall unsafe "flint/fmpz.h fmpz_set_si" fmpz_set_si :: Ptr FMPZ -> SLong -> IO ()
foreign import ccall unsafe "flint/fmpz.h fmpz_get_si" fmpz_get_si :: Ptr FMPZ -> IO SLong
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_clear" fmpz_mat_clear :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_print_pretty" fmpz_mat_print_pretty :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_mul" fmpz_mat_mul :: Ptr FMPZMat -> Ptr FMPZMat -> Ptr FMPZMat -> IO ()

-- fmpz_mat_window_init(fmpz_mat_t window, const fmpz_mat_t mat, slong r1, slong c1, slong r2, slong c2)
--
-- Initializes the matrix window to be an r2 - r1 by c2 - c1 submatrix
-- of mat whose (0,0) entry is the (r1, c1) entry of mat. The memory
-- for the elements of window is shared with mat.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_init" fmpz_mat_window_init :: Ptr FMPZMat -> Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> IO ()

-- Frees the window (leaving underlying matrix alone).
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_clear" fmpz_mat_window_clear :: Ptr FMPZMat -> IO ()

-- r <- fmp_mat_rref B den A
--
-- Uses fraction-free Gauss-Jordan elimination to set (B, den) to the
-- reduced row echelon form of A and returns the rank of A. Aliasing
-- of A and B is allowed. r is rank of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rref" fmpz_mat_rref :: Ptr FMPZMat -> Ptr FMPZ -> Ptr FMPZMat -> IO SLong

-- r <- fmp_mat_inv B den A
--
-- Returns 1 if @A@ is nonsingular and 0 if @A@ is singular.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_inv" fmpz_mat_inv :: Ptr FMPZMat -> Ptr FMPZ -> Ptr FMPZMat -> IO CInt

-- fmpz_mat_hnf H A
--
-- H is the Hermite Normal Form of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_hnf" fmpz_mat_hnf :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rank" fmpz_mat_rank :: Ptr FMPZMat -> IO SLong

rref :: H.Matrix Double -> (H.Matrix Double, Int, Int)
rref :: Matrix Double -> (Matrix Double, Int, Int)
rref Matrix Double
m = IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int)
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      (Ptr FMPZ -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZ -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr FMPZ -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZ
den -> do
        let numRows :: SLong
numRows = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
        let numCols :: SLong
numCols = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
        Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM SLong
numRows SLong
numCols
        Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM SLong
numRows SLong
numCols
        [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
          [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
            Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
inputM SLong
i SLong
j
            Ptr FMPZ -> SLong -> IO ()
fmpz_set_si Ptr FMPZ
e (Double -> SLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
i, SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
j)))
        SLong
r <- Ptr FMPZMat -> Ptr FMPZ -> Ptr FMPZMat -> IO SLong
fmpz_mat_rref Ptr FMPZMat
outputM Ptr FMPZ
den Ptr FMPZMat
inputM
        FMPZ
d <- Ptr FMPZ -> IO FMPZ
forall a. Storable a => Ptr a -> IO a
peek Ptr FMPZ
den

        -- DEBUG:
        -- fmpz_mat_print_pretty outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        [[Double]]
lists <- [SLong] -> (SLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO [Double]) -> IO [[Double]])
-> (SLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
          [SLong] -> (SLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO Double) -> IO [Double])
-> (SLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
            Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
outputM SLong
i SLong
j
            SLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (SLong -> Double) -> IO SLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr FMPZ -> IO SLong
fmpz_get_si Ptr FMPZ
e
        let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
        (Matrix Double, Int, Int) -> IO (Matrix Double, Int, Int)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double
m', FMPZ -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral FMPZ
d, SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
r)

hnf :: H.Matrix Double -> H.Matrix Double
hnf :: Matrix Double -> Matrix Double
hnf Matrix Double
m = IO (Matrix Double) -> Matrix Double
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double) -> Matrix Double)
-> IO (Matrix Double) -> Matrix Double
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double))
-> (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double))
-> (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      let numRows :: SLong
numRows = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
      let numCols :: SLong
numCols = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
      Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM SLong
numRows SLong
numCols
      Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM SLong
numRows SLong
numCols
      [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
        [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
          Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
inputM SLong
i SLong
j
          Ptr FMPZ -> SLong -> IO ()
fmpz_set_si Ptr FMPZ
e (Double -> SLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
i, SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
j)))
      Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM Ptr FMPZMat
inputM
      SLong
r <- Ptr FMPZMat -> IO SLong
fmpz_mat_rank Ptr FMPZMat
outputM

      -- DEBUG:
      -- fmpz_mat_print_pretty outputM
      -- traceM $ "rank = " ++ show r
      --

      [[Double]]
lists <- [SLong] -> (SLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
rSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO [Double]) -> IO [[Double]])
-> (SLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
        [SLong] -> (SLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO Double) -> IO [Double])
-> (SLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
          Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
outputM SLong
i SLong
j
          SLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (SLong -> Double) -> IO SLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr FMPZ -> IO SLong
fmpz_get_si Ptr FMPZ
e
      let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
      Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
      Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
      Matrix Double -> IO (Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix Double
m'

inv :: H.Matrix Double -> Maybe (H.Matrix Double, Int)
inv :: Matrix Double -> Maybe (Matrix Double, Int)
inv Matrix Double
m = IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int))
-> IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int)
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      (Ptr FMPZ -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZ -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr FMPZ -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZ
den -> do
        let numRows :: SLong
numRows = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
        let numCols :: SLong
numCols = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
        Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM SLong
numRows SLong
numCols
        Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM SLong
numRows SLong
numCols
        [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
          [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
            Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
inputM SLong
i SLong
j
            Ptr FMPZ -> SLong -> IO ()
fmpz_set_si Ptr FMPZ
e (Double -> SLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
i, SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
j)))
        CInt
r <- Ptr FMPZMat -> Ptr FMPZ -> Ptr FMPZMat -> IO CInt
fmpz_mat_inv Ptr FMPZMat
outputM Ptr FMPZ
den Ptr FMPZMat
inputM
        FMPZ
d <- Ptr FMPZ -> IO FMPZ
forall a. Storable a => Ptr a -> IO a
peek Ptr FMPZ
den

        -- DEBUG:
        Ptr FMPZMat -> IO ()
fmpz_mat_print_pretty Ptr FMPZMat
outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        [[Double]]
lists <- [SLong] -> (SLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO [Double]) -> IO [[Double]])
-> (SLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
          [SLong] -> (SLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO Double) -> IO [Double])
-> (SLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
            Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
outputM SLong
i SLong
j
            SLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (SLong -> Double) -> IO SLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr FMPZ -> IO SLong
fmpz_get_si Ptr FMPZ
e
        let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
        if CInt
r CInt -> CInt -> Bool
forall a. Eq a => a -> a -> Bool
== CInt
1 then
          Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int)))
-> Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ (Matrix Double, Int) -> Maybe (Matrix Double, Int)
forall a. a -> Maybe a
Just (Matrix Double
m', FMPZ -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral FMPZ
d)
        else
          Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Matrix Double, Int)
forall a. Maybe a
Nothing

withMatrix :: H.Matrix Double -> ((SLong, SLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix :: Matrix Double -> ((SLong, SLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix Matrix Double
m (SLong, SLong, Ptr FMPZMat) -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    let numRows :: SLong
numRows = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
    let numCols :: SLong
numCols = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> SLong) -> Int -> SLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
    Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM SLong
numRows SLong
numCols
    [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0 .. SLong
numRows SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
- SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i ->
      [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0 .. SLong
numCols SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
- SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
        Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
outputM SLong
i SLong
j
        Ptr FMPZ -> SLong -> IO ()
fmpz_set_si Ptr FMPZ
e (Double -> SLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
i, SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
j)) :: SLong)
    b
x <- (SLong, SLong, Ptr FMPZMat) -> IO b
f (SLong
numRows, SLong
numCols, Ptr FMPZMat
outputM)
    Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

withBlankMatrix :: SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix :: SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix SLong
numRows SLong
numCols Ptr FMPZMat -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    Ptr FMPZMat -> SLong -> SLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM (SLong -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
numRows) (SLong -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral SLong
numCols)
    b
x <- Ptr FMPZMat -> IO b
f Ptr FMPZMat
outputM
    Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

withWindow :: Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withWindow :: Ptr FMPZMat
-> SLong
-> SLong
-> SLong
-> SLong
-> (Ptr FMPZMat -> IO b)
-> IO b
withWindow Ptr FMPZMat
underM SLong
r1 SLong
c1 SLong
r2 SLong
c2 Ptr FMPZMat -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
window -> do
    Ptr FMPZMat
-> Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> IO ()
fmpz_mat_window_init Ptr FMPZMat
window Ptr FMPZMat
underM SLong
r1 SLong
c1 SLong
r2 SLong
c2
    b
x <- Ptr FMPZMat -> IO b
f Ptr FMPZMat
window
    Ptr FMPZMat -> IO ()
fmpz_mat_window_clear Ptr FMPZMat
window
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

flintToHMatrix :: SLong -> SLong -> Ptr FMPZMat -> IO (H.Matrix Double)
flintToHMatrix :: SLong -> SLong -> Ptr FMPZMat -> IO (Matrix Double)
flintToHMatrix SLong
numRows SLong
numCols Ptr FMPZMat
flintM = do
  [[Double]]
lists <- [SLong] -> (SLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numRowsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO [Double]) -> IO [[Double]])
-> (SLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
    [SLong] -> (SLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO Double) -> IO [Double])
-> (SLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ SLong
j -> do
      Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
flintM SLong
i SLong
j
      SLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (SLong -> Double) -> IO SLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr FMPZ -> IO SLong
fmpz_get_si Ptr FMPZ
e
  Matrix Double -> IO (Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double -> IO (Matrix Double))
-> Matrix Double -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists

pokeM :: Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM :: Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
flintM SLong
i SLong
j SLong
v = do
  Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
flintM SLong
i SLong
j
  Ptr FMPZ -> SLong -> IO ()
fmpz_set_si Ptr FMPZ
e SLong
v

peekM :: Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM :: Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
flintM SLong
i SLong
j = do
  Ptr FMPZ
e <- Ptr FMPZMat -> SLong -> SLong -> IO (Ptr FMPZ)
fmpz_mat_entry Ptr FMPZMat
flintM SLong
i SLong
j
  Ptr FMPZ -> IO SLong
fmpz_get_si Ptr FMPZ
e

copyMatrix :: Ptr FMPZMat -> Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> IO ()
copyMatrix :: Ptr FMPZMat
-> Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> IO ()
copyMatrix Ptr FMPZMat
m1 Ptr FMPZMat
m2 SLong
r1 SLong
c1 SLong
r2 SLong
c2 =
  [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
r1 .. SLong
r2SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i ->
    [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
c1 .. SLong
c2SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j ->
      Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
m2 SLong
i SLong
j IO SLong -> (SLong -> IO ()) -> IO ()
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
m1 SLong
i SLong
j


-- copyMatrix' m1 m2 r1 c1 r2 c2 =
--   withWindow m2 r1 c1 r2 c2 $ \ w2 ->
--     withWindow m1 r1 c1 r2 c2 $ \ w1 -> do
--       fmpz_mat_set w1 w2

--------------------------------------------------
-- 'normalising' Hermite Normal Form, for lack of a better name
--
-- The problem with HNF is that it will happily return a matrix where
-- the leading-coefficient diagonal contains integers > 1.
--
-- In some cases the leading-coefficient does not divide some of the
-- remaining numbers in the row.
--
-- This corresponds to the case where one of the solutions would be
-- fractional if we were dealing with rational matrices.
--
-- But we don't want fractional solutions. We need to bump the matrix
-- so that this situation does not arise.
--
-- This tends to happen due to implicit polymorphism.
--
-- We do this by adding another column that copies the column with the
-- problematic leading-coefficient.
--
-- We then set up a new row expressing the constraint that the old
-- column is equal to the new column, 1:1.
--
-- This prevents a 'fractional' solution.
--
-- Along the way we look for leading-coefficients > 1 that do divide
-- their entire row. Then we simply apply elementary row scaling to
-- the row so that the leading-coefficients is 1.
--
-- The result is a matrix where the leading-coefficients of all the
-- rows that matter are 1. It's not quite a 'reduced' form though
-- because it can be bigger than the original matrix.
--
-- As a result we also return a list of columns that were cloned this
-- way.
--
-- The resulting matrix can have non-zeroes above the
-- leading-coefficients in the same column, so I apply some elementary
-- row operations to fix that and bring things into RREF.
--
-- Running time is computation of Hermite Normal Form twice, plus
-- construction of a slightly larger matrix (possibly), plus scanning
-- through the matrix for leading-coefficients, and then again to
-- divide them out sometimes.

normHNF :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF :: Matrix Double -> (Matrix Double, [Int])
normHNF Matrix Double
m
  | (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m, Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m) (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
1, Int
1) = (Int -> Matrix Double
forall a. (Num a, Element a) => Int -> Matrix a
H.ident (if Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
H.atIndex Matrix Double
m (Int
0, Int
0) Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 then Int
1 else Int
0), [])
  | Bool
otherwise = (Matrix Double
m', [Int]
indices)
  where
    numCols :: Int
numCols = Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
    indexLookup :: Int -> Int
indexLookup Int
j | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
numCols = Int
j
                  | Bool
otherwise = [Int]
indices [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
numCols)
    indices :: [Int]
indices = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Int
indexLookup ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ ((Matrix Double, [Int]) -> [Int])
-> [(Matrix Double, [Int])] -> [Int]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Matrix Double, [Int]) -> [Int]
forall a b. (a, b) -> b
snd [(Matrix Double, [Int])]
rs1
    ([(Matrix Double, [Int])]
rs1, (Matrix Double
m',[]):[(Matrix Double, [Int])]
_) = ((Matrix Double, [Int]) -> Bool)
-> [(Matrix Double, [Int])]
-> ([(Matrix Double, [Int])], [(Matrix Double, [Int])])
forall a. (a -> Bool) -> [a] -> ([a], [a])
break ([Int] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([Int] -> Bool)
-> ((Matrix Double, [Int]) -> [Int])
-> (Matrix Double, [Int])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Matrix Double, [Int]) -> [Int]
forall a b. (a, b) -> b
snd) [(Matrix Double, [Int])]
results
    results :: [(Matrix Double, [Int])]
results = [(Matrix Double, [Int])] -> [(Matrix Double, [Int])]
forall a. [a] -> [a]
tail ([(Matrix Double, [Int])] -> [(Matrix Double, [Int])])
-> [(Matrix Double, [Int])] -> [(Matrix Double, [Int])]
forall a b. (a -> b) -> a -> b
$ ((Matrix Double, [Int]) -> (Matrix Double, [Int]))
-> (Matrix Double, [Int]) -> [(Matrix Double, [Int])]
forall a. (a -> a) -> a -> [a]
iterate (Matrix Double -> (Matrix Double, [Int])
normHNF' (Matrix Double -> (Matrix Double, [Int]))
-> ((Matrix Double, [Int]) -> Matrix Double)
-> (Matrix Double, [Int])
-> (Matrix Double, [Int])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Matrix Double, [Int]) -> Matrix Double
forall a b. (a, b) -> a
fst) (Matrix Double
m, [])

normHNF' :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF' :: Matrix Double -> (Matrix Double, [Int])
normHNF' Matrix Double
m = ([SLong] -> [Int])
-> (Matrix Double, [SLong]) -> (Matrix Double, [Int])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((SLong -> Int) -> [SLong] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map SLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral) ((Matrix Double, [SLong]) -> (Matrix Double, [Int]))
-> (IO (Matrix Double, [SLong]) -> (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
-> (Matrix Double, [Int])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix Double, [SLong]) -> (Matrix Double, [SLong])
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double, [SLong]) -> (Matrix Double, [Int]))
-> IO (Matrix Double, [SLong]) -> (Matrix Double, [Int])
forall a b. (a -> b) -> a -> b
$ Matrix Double
-> ((SLong, SLong, Ptr FMPZMat) -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall b.
Matrix Double -> ((SLong, SLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix Matrix Double
m (SLong, SLong, Ptr FMPZMat) -> IO (Matrix Double, [SLong])
normhnf

normhnf :: (SLong, SLong, Ptr FMPZMat) -> IO (H.Matrix Double, [SLong])
normhnf :: (SLong, SLong, Ptr FMPZMat) -> IO (Matrix Double, [SLong])
normhnf (SLong
numRows, SLong
numCols, Ptr FMPZMat
inputM) = do
  SLong
-> SLong
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall b. SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix SLong
numRows SLong
numCols ((Ptr FMPZMat -> IO (Matrix Double, [SLong]))
 -> IO (Matrix Double, [SLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM Ptr FMPZMat
inputM
    SLong
rank <- Ptr FMPZMat -> IO SLong
fmpz_mat_rank Ptr FMPZMat
outputM
    -- scale the rows so that it is a RREF, returning the indices
    -- where leading-coefficients had to be scaled to 1.
    [(SLong, SLong)]
indices <- Ptr FMPZMat -> SLong -> SLong -> IO [(SLong, SLong)]
elemrowscale Ptr FMPZMat
outputM SLong
rank SLong
numCols
    -- HNF allows non-zeroes above the leading-coefficients that
    -- were greater than 1, so also fix that to bring into RREF.
    Ptr FMPZMat -> SLong -> SLong -> [(SLong, SLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM SLong
rank SLong
numCols [(SLong, SLong)]
indices

    -- column indices of leading co-efficients > 1
    [((SLong, SLong), [SLong])]
lcoefs <- (((SLong, SLong), [SLong]) -> Bool)
-> [((SLong, SLong), [SLong])] -> [((SLong, SLong), [SLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter ((SLong -> SLong -> Bool
forall a. Ord a => a -> a -> Bool
> SLong
1) (SLong -> Bool)
-> (((SLong, SLong), [SLong]) -> SLong)
-> ((SLong, SLong), [SLong])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [SLong] -> SLong
forall a. [a] -> a
head ([SLong] -> SLong)
-> (((SLong, SLong), [SLong]) -> [SLong])
-> ((SLong, SLong), [SLong])
-> SLong
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((SLong, SLong), [SLong]) -> [SLong]
forall a b. (a, b) -> b
snd) ([((SLong, SLong), [SLong])] -> [((SLong, SLong), [SLong])])
-> IO [((SLong, SLong), [SLong])] -> IO [((SLong, SLong), [SLong])]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [SLong]
-> (SLong -> IO ((SLong, SLong), [SLong]))
-> IO [((SLong, SLong), [SLong])]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0 .. SLong
rankSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] (\ SLong
i -> do
      [(SLong, SLong)]
cs <- [SLong] -> [SLong] -> [(SLong, SLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [SLong
0..] ([SLong] -> [(SLong, SLong)]) -> IO [SLong] -> IO [(SLong, SLong)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` [IO SLong] -> IO [SLong]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [ Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j | SLong
j <- [SLong
0 .. SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ]
      let ([(SLong, SLong)]
_, (SLong
j, SLong
lcoef):[(SLong, SLong)]
rest) = ((SLong, SLong) -> Bool)
-> [(SLong, SLong)] -> ([(SLong, SLong)], [(SLong, SLong)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span ((SLong -> SLong -> Bool
forall a. Eq a => a -> a -> Bool
== SLong
0) (SLong -> Bool)
-> ((SLong, SLong) -> SLong) -> (SLong, SLong) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (SLong, SLong) -> SLong
forall a b. (a, b) -> b
snd) [(SLong, SLong)]
cs
      ((SLong, SLong), [SLong]) -> IO ((SLong, SLong), [SLong])
forall (m :: * -> *) a. Monad m => a -> m a
return ((SLong
i, SLong
j), SLong
lcoefSLong -> [SLong] -> [SLong]
forall a. a -> [a] -> [a]
:((SLong, SLong) -> SLong) -> [(SLong, SLong)] -> [SLong]
forall a b. (a -> b) -> [a] -> [b]
map (SLong, SLong) -> SLong
forall a b. (a, b) -> b
snd [(SLong, SLong)]
rest))

    -- split the identified rows into two categories: those that can
    -- be divided out by the leading co-efficient, and those that
    -- cannot.
    let ([((SLong, SLong), [SLong])]
multCands, [((SLong, SLong), [SLong])]
consCands) = (((SLong, SLong), [SLong]) -> Bool)
-> [((SLong, SLong), [SLong])]
-> ([((SLong, SLong), [SLong])], [((SLong, SLong), [SLong])])
forall a. (a -> Bool) -> [a] -> ([a], [a])
partition (\ ((SLong, SLong)
_, SLong
lcoef:[SLong]
rest) -> (SLong -> Bool) -> [SLong] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((SLong -> SLong -> Bool
forall a. Eq a => a -> a -> Bool
== SLong
0) (SLong -> Bool) -> (SLong -> SLong) -> SLong -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (SLong -> SLong -> SLong
forall a. Integral a => a -> a -> a
`rem` SLong
lcoef)) [SLong]
rest) [((SLong, SLong), [SLong])]
lcoefs

    -- apply elementary row scaling
    [((SLong, SLong), [SLong])]
-> (((SLong, SLong), [SLong]) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [((SLong, SLong), [SLong])]
multCands ((((SLong, SLong), [SLong]) -> IO ()) -> IO ())
-> (((SLong, SLong), [SLong]) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ ((SLong
i, SLong
j), SLong
lcoef:[SLong]
_) ->
      [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
j..SLong
numCols SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
- SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j' -> do
        SLong
x <- Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j'
        Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
outputM SLong
i SLong
j' (SLong -> IO ()) -> SLong -> IO ()
forall a b. (a -> b) -> a -> b
$ SLong
x SLong -> SLong -> SLong
forall a. Integral a => a -> a -> a
`div` SLong
lcoef

    -- identify columns that need additional constraints generated and
    -- their associated value d, which is the value of the non-leading
    -- co-efficient divided by the GCD of the row.
    let genColCons :: ((a, a), [b]) -> (a, b)
genColCons ((a
_, a
j), b
lcoef:[b]
rest) = (a
j, b
minNLcoef b -> b -> b
forall a. Integral a => a -> a -> a
`div` b -> b -> b
forall a. Integral a => a -> a -> a
gcd b
lcoef b
minNLcoef)
          where
            restABS :: [b]
restABS   = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
forall a. Num a => a -> a
abs [b]
rest
            minNLcoef :: b
minNLcoef = [b] -> b
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ((b -> Bool) -> [b] -> [b]
forall a. (a -> Bool) -> [a] -> [a]
filter (b -> b -> Bool
forall a. Eq a => a -> a -> Bool
/= b
0) [b]
restABS)
        genColCons ((a, a), [b])
_ = [Char] -> (a, b)
forall a. HasCallStack => [Char] -> a
error [Char]
"normhnf: genColCons: impossible: missing leading co-efficient"

    let consCols :: [(SLong, SLong)]
consCols = (((SLong, SLong), [SLong]) -> (SLong, SLong))
-> [((SLong, SLong), [SLong])] -> [(SLong, SLong)]
forall a b. (a -> b) -> [a] -> [b]
map ((SLong, SLong), [SLong]) -> (SLong, SLong)
forall b a a. Integral b => ((a, a), [b]) -> (a, b)
genColCons [((SLong, SLong), [SLong])]
consCands

    -- generate operations that poke 1.0 and -d into columns of the
    -- given row (matrix parameter supplied later); d is the value of
    -- the non-leading co-efficient that isn't divisible by the
    -- leading-coefficient, divided by the GCD of the row.
    let ops :: [Ptr FMPZMat -> SLong -> IO ()]
ops = [ \ Ptr FMPZMat
flintM SLong
i -> do Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
flintM SLong
i SLong
j SLong
1
                                 Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
flintM SLong
i (SLong
numCols SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
+ SLong
k) (-SLong
d)
              | ((SLong
j, SLong
d), SLong
k) <- [(SLong, SLong)] -> [SLong] -> [((SLong, SLong), SLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(SLong, SLong)]
consCols [SLong
0..] ]
    let numOps :: SLong
numOps = Int -> SLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Ptr FMPZMat -> SLong -> IO ()] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Ptr FMPZMat -> SLong -> IO ()]
ops)
    let numRows' :: SLong
numRows' = SLong
numOps SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
+ SLong
rank
    let numCols' :: SLong
numCols' = SLong
numOps SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
+ SLong
numCols

    SLong
-> SLong
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall b. SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix SLong
numRows' SLong
numCols' ((Ptr FMPZMat -> IO (Matrix Double, [SLong]))
 -> IO (Matrix Double, [SLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM' -> do
      -- create larger matrix containing outputM as submatrix
      Ptr FMPZMat
-> Ptr FMPZMat -> SLong -> SLong -> SLong -> SLong -> IO ()
copyMatrix Ptr FMPZMat
outputM' Ptr FMPZMat
outputM SLong
0 SLong
0 SLong
rank SLong
numCols

      -- apply the operations on the extra space to write the
      -- additional rows and columns providing the additional
      -- constraints
      [(SLong, Ptr FMPZMat -> SLong -> IO ())]
-> ((SLong, Ptr FMPZMat -> SLong -> IO ()) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ ([SLong]
-> [Ptr FMPZMat -> SLong -> IO ()]
-> [(SLong, Ptr FMPZMat -> SLong -> IO ())]
forall a b. [a] -> [b] -> [(a, b)]
zip [SLong
rank..] [Ptr FMPZMat -> SLong -> IO ()]
ops) (((SLong, Ptr FMPZMat -> SLong -> IO ()) -> IO ()) -> IO ())
-> ((SLong, Ptr FMPZMat -> SLong -> IO ()) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ (SLong
i, Ptr FMPZMat -> SLong -> IO ()
op) -> Ptr FMPZMat -> SLong -> IO ()
op Ptr FMPZMat
outputM' SLong
i

      -- re-run HNF
      SLong
-> SLong
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall b. SLong -> SLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix SLong
numRows' SLong
numCols' ((Ptr FMPZMat -> IO (Matrix Double, [SLong]))
 -> IO (Matrix Double, [SLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [SLong]))
-> IO (Matrix Double, [SLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM'' -> do
        Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM'' Ptr FMPZMat
outputM'
        SLong
rank' <- Ptr FMPZMat -> IO SLong
fmpz_mat_rank Ptr FMPZMat
outputM''
        -- scale the rows so that it is a RREF, returning the indices
        -- where leading-coefficients had to be scaled to 1.
        [(SLong, SLong)]
indices' <- Ptr FMPZMat -> SLong -> SLong -> IO [(SLong, SLong)]
elemrowscale Ptr FMPZMat
outputM'' SLong
rank' SLong
numCols'
        -- HNF allows non-zeroes above the leading-coefficients that
        -- were greater than 1, so also fix that to bring into RREF.
        Ptr FMPZMat -> SLong -> SLong -> [(SLong, SLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM'' SLong
rank' SLong
numCols' [(SLong, SLong)]
indices'
        -- convert back to HMatrix form
        Matrix Double
h <- SLong -> SLong -> Ptr FMPZMat -> IO (Matrix Double)
flintToHMatrix SLong
rank' SLong
numCols' Ptr FMPZMat
outputM''
        (Matrix Double, [SLong]) -> IO (Matrix Double, [SLong])
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double
h, ((SLong, SLong) -> SLong) -> [(SLong, SLong)] -> [SLong]
forall a b. (a -> b) -> [a] -> [b]
map (SLong, SLong) -> SLong
forall a b. (a, b) -> a
fst [(SLong, SLong)]
consCols)

-- find leading-coefficients that are greater than 1 and scale those
-- rows accordingly to reach RREF
--
-- precondition: matrix outputM is in HNF
elemrowscale :: Ptr FMPZMat -> SLong -> SLong -> IO [(SLong, SLong)]
elemrowscale :: Ptr FMPZMat -> SLong -> SLong -> IO [(SLong, SLong)]
elemrowscale Ptr FMPZMat
outputM SLong
rank SLong
numCols = do
  -- column indices of leading co-efficients > 1
  [((SLong, SLong), [SLong])]
lcoefs <- (((SLong, SLong), [SLong]) -> Bool)
-> [((SLong, SLong), [SLong])] -> [((SLong, SLong), [SLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter ((SLong -> SLong -> Bool
forall a. Ord a => a -> a -> Bool
> SLong
1) (SLong -> Bool)
-> (((SLong, SLong), [SLong]) -> SLong)
-> ((SLong, SLong), [SLong])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [SLong] -> SLong
forall a. [a] -> a
head ([SLong] -> SLong)
-> (((SLong, SLong), [SLong]) -> [SLong])
-> ((SLong, SLong), [SLong])
-> SLong
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((SLong, SLong), [SLong]) -> [SLong]
forall a b. (a, b) -> b
snd) ([((SLong, SLong), [SLong])] -> [((SLong, SLong), [SLong])])
-> IO [((SLong, SLong), [SLong])] -> IO [((SLong, SLong), [SLong])]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [SLong]
-> (SLong -> IO ((SLong, SLong), [SLong]))
-> IO [((SLong, SLong), [SLong])]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [SLong
0 .. SLong
rankSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] (\ SLong
i -> do
    [(SLong, SLong)]
cs <- [SLong] -> [SLong] -> [(SLong, SLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [SLong
0..] ([SLong] -> [(SLong, SLong)]) -> IO [SLong] -> IO [(SLong, SLong)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` [IO SLong] -> IO [SLong]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [ Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j | SLong
j <- [SLong
0 .. SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ]
    let ([(SLong, SLong)]
_, (SLong
j, SLong
lcoef):[(SLong, SLong)]
rest) = ((SLong, SLong) -> Bool)
-> [(SLong, SLong)] -> ([(SLong, SLong)], [(SLong, SLong)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span ((SLong -> SLong -> Bool
forall a. Eq a => a -> a -> Bool
== SLong
0) (SLong -> Bool)
-> ((SLong, SLong) -> SLong) -> (SLong, SLong) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (SLong, SLong) -> SLong
forall a b. (a, b) -> b
snd) [(SLong, SLong)]
cs
    ((SLong, SLong), [SLong]) -> IO ((SLong, SLong), [SLong])
forall (m :: * -> *) a. Monad m => a -> m a
return ((SLong
i, SLong
j), SLong
lcoefSLong -> [SLong] -> [SLong]
forall a. a -> [a] -> [a]
:((SLong, SLong) -> SLong) -> [(SLong, SLong)] -> [SLong]
forall a b. (a -> b) -> [a] -> [b]
map (SLong, SLong) -> SLong
forall a b. (a, b) -> b
snd [(SLong, SLong)]
rest))

  let multCands :: [((SLong, SLong), [SLong])]
multCands = (((SLong, SLong), [SLong]) -> Bool)
-> [((SLong, SLong), [SLong])] -> [((SLong, SLong), [SLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter (\ ((SLong, SLong)
_, SLong
lcoef:[SLong]
rest) -> (SLong -> Bool) -> [SLong] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((SLong -> SLong -> Bool
forall a. Eq a => a -> a -> Bool
== SLong
0) (SLong -> Bool) -> (SLong -> SLong) -> SLong -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (SLong -> SLong -> SLong
forall a. Integral a => a -> a -> a
`rem` SLong
lcoef)) [SLong]
rest) [((SLong, SLong), [SLong])]
lcoefs

  -- apply elementary row scaling
  [((SLong, SLong), [SLong])]
-> (((SLong, SLong), [SLong]) -> IO (SLong, SLong))
-> IO [(SLong, SLong)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [((SLong, SLong), [SLong])]
multCands ((((SLong, SLong), [SLong]) -> IO (SLong, SLong))
 -> IO [(SLong, SLong)])
-> (((SLong, SLong), [SLong]) -> IO (SLong, SLong))
-> IO [(SLong, SLong)]
forall a b. (a -> b) -> a -> b
$ \ ((SLong
i, SLong
j), SLong
lcoef:[SLong]
_) -> do
    [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
j..SLong
numCols SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
- SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j' -> do
      SLong
x <- Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j'
      Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
outputM SLong
i SLong
j' (SLong -> IO ()) -> SLong -> IO ()
forall a b. (a -> b) -> a -> b
$ SLong
x SLong -> SLong -> SLong
forall a. Integral a => a -> a -> a
`div` SLong
lcoef
    (SLong, SLong) -> IO (SLong, SLong)
forall (m :: * -> *) a. Monad m => a -> m a
return (SLong
i, SLong
j)

-- use indices to guide elementary row additions to reach RREF
--
-- precondition: matrix outputM is in HNF save for work done by
-- elemrowscale, and indices is a list of coordinates where we have
-- just scaled the leading-coefficient to 1 and now we must look to
-- see if there are any non-zeroes in the column above the
-- leading-coefficient, because that is allowed by HNF.
elemrowadds :: Ptr FMPZMat -> SLong -> SLong -> [(SLong, SLong)] -> IO ()
elemrowadds :: Ptr FMPZMat -> SLong -> SLong -> [(SLong, SLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM SLong
_ SLong
numCols [(SLong, SLong)]
indices = do
  -- look for non-zero members of the columns above the
  -- leading-coefficient and wipe them out.
  [(SLong, SLong)] -> ((SLong, SLong) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(SLong, SLong)]
indices (((SLong, SLong) -> IO ()) -> IO ())
-> ((SLong, SLong) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ (SLong
lcI, SLong
lcJ) -> do
    let j :: SLong
j = SLong
lcJ
    [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
0..SLong
lcISLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
i -> do
      -- (i, j) ranges over the elements of the column above leading
      -- co-efficient (lcI, lcJ).
      SLong
x <- Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j
      if SLong
x SLong -> SLong -> Bool
forall a. Eq a => a -> a -> Bool
== SLong
0 then
        () -> IO ()
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- nothing to do at row i
      else do
        -- (i, j) is non-zero and row i must be cancelled x times
        let sf :: SLong
sf = SLong
x -- scaling factor = non-zero magnitude
        [SLong] -> (SLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [SLong
lcJ..SLong
numColsSLong -> SLong -> SLong
forall a. Num a => a -> a -> a
-SLong
1] ((SLong -> IO ()) -> IO ()) -> (SLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ SLong
j' -> do
          -- (i,   j') ranges over the row where we discovered a non-zero
          -- (lcI, j') ranges over the row with the leading co-efficient
          SLong
x1 <- Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
i SLong
j'
          SLong
x2 <- Ptr FMPZMat -> SLong -> SLong -> IO SLong
peekM Ptr FMPZMat
outputM SLong
lcI SLong
j'
          -- add lower row scaled by sf to upper row
          Ptr FMPZMat -> SLong -> SLong -> SLong -> IO ()
pokeM Ptr FMPZMat
outputM SLong
i SLong
j' (SLong
x1 SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
- SLong
x2 SLong -> SLong -> SLong
forall a. Num a => a -> a -> a
* SLong
sf)

--------------------------------------------------------------------------------

-- testFlint :: IO ()
-- testFlint = do
--   traceM "***********************8 testFlint 8***********************"
--   alloca $ \ a -> do
--     alloca $ \ b -> do
--       alloca $ \ den -> do
--         let n = 10
--         fmpz_mat_init a n n
--         fmpz_mat_init b n n
--         forM_ [0..n-1] $ \ i -> do
--           forM_ [0..n-1] $ \ j -> do
--             e <- fmpz_mat_entry a i j
--             fmpz_set_si e (2 * i + j)

--         -- fmpz_mat_mul b a a
--         r <- fmpz_mat_rref b den a
--         fmpz_mat_print_pretty a
--         fmpz_mat_print_pretty b
--         d <- peek den
--         traceM $ "r = " ++ show r ++ " den = " ++ show d
--         fmpz_mat_hnf b a
--         fmpz_mat_print_pretty b
--         fmpz_mat_clear a
--         fmpz_mat_clear b

--------------------------------------------------------------------------------

-- m1 :: H.Matrix Double
-- m1 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -4.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  4.0, -3.0 ]

-- m2 :: H.Matrix Double
-- m2 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -6.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  6.0, -4.0 ]