0votos

Mejor factorización triple en Haskell

por josejuan hace 5 años

Siempre que hay que empaquetar algo, puede usarse programación lineal. En este caso es mixta. Se transforman las ecuaciones pasando a logaritmos, de esta forma transformamos el producto de potencias (la factorización) en suma de logaritmos y así poder aplicar LP. "Parece" funcionar un poco mejor, pero sigue sin ser eficiente.

Calcular la factorización de n en un producto de tres números (a,b,c) tales que 1 ≤ a ≤ b ≤ c y c/a es mínima.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import Data.List 
import System.Environment 
import Control.Monad 
import Data.LinearProgram 
import Data.LinearProgram.GLPK 
import qualified Data.Map as M 
import Data.Numbers.Primes 
 
{- 
 
    Sea la descomposición en factores primos 
 
        n = A^a B^b C^c 
 
    (con A, B, ... primos y a, b, ... naturales) 
 
    Tomando logaritmos tenemos que 
 
        ln n = a ln A + b ln B + ... 
 
    Se pide obtener una descomposición 
 
        n = A^(a1+a2+a3) B^(b1+b2+b3) ... 
 
    Sean los grupos 1, 2 y 3 correspondientes con `a`, `b` y `c` 
    del enunciado, entonces es 
 
        G1 = A^a1 B^b1 ... 
        G2 = A^a2 B^b2 ... 
        G3 = A^a3 B^b3 ... 
 
    Tales que 
 
        1 <= G1 <= G2 <= G3 
 
    que corresponde con la restricción del enunciado 
 
        1 <= a <= b <= c 
 
    Usando la representación en logaritmos, las ecuaciones 
    quedan 
 
        g1 = a1 ln A + b1 ln B + ... 
        g2 = a2 ln A + b2 ln B + ... 
        g3 = a3 ln A + b3 ln B + ... 
 
    y las restricciones 
 
         0 <= g1 
        g1 <= g2 
        g2 <= g3 
 
    El objetivo original de optimización es 
 
        minimizar `c` / `a` 
 
    Pero tomando logaritmos queda 
 
        mín: g3 - g1 
 
    es decir 
 
        mín: (a3 - a1) ln A + (b3 - b1) ln B + ... 
 
-} 
 
-- Las incógnitas Aij vienen dadas por cada grupo i=1,2,3 y cada j=primo 
a :: Int -> Int -> (Int, Int) 
a grupo primo = (grupo, primo) 
 
matchProblem :: -- (Floating c, Num (M.Map (Int, Int) c), Group c) => 
                [(Int, Int)] -> LP (Int, Int) Double 
matchProblem fs = 
  execLPM $ 
  do 
 
    setDirection Min 
 
    let geng g = linCombination [(log $ fromIntegral primo, a g primo) | (primo, _) <- fs] 
        [g1, g2, g3] = map geng [1..3] 
        n = fromIntegral $ sum [ primo^potencia | (primo, potencia) <- fs ] 
 
    forM_ fs $ \(primo, potencia) -> do 
                                        (var (a 1 primo) ^+^ var (a 2 primo) ^+^ var (a 3 primo)) `equalTo` (fromIntegral potencia) 
                                        forM_ [1..3] $ \g -> do 
                                                                setVarKind (a g primo) IntVar 
                                                                (var $ a g primo) `geqTo` 0 
 
    (       g1) `geqTo` 0 
    (g2 ^-^ g1) `geqTo` 0 
    (g3 ^-^ g2) `geqTo` 0 
 
    setObjective $ g3 ^-^ g1 
 
solveProblem fs = do 
  (_, Just (_, m)) <- glpSolveVars (mipDefaults { msgLev = MsgErr }) $ matchProblem fs 
  return 
    $ map (map (\((_,f),p) -> (f, floor p)) . filter (\(_, p) -> p > 0)) . 
      groupBy (\((a,_),_) ((b,_),_) -> a == b) . 
      sortBy (\((a,_),_) ((b,_),_) -> a `compare` b) . 
      M.toList $ m 
 
solveProblem' fs = solveProblem fs >>= return . map (product . map (uncurry (^))) 
 
potenciasPrimas = map (\ps -> (head ps, length ps)) . group . primeFactors 
 
solveProblem'' = solveProblem' . potenciasPrimas 
 
 
-- Para probar podría hacerse 
main = do 
  xs <- getArgs 
  rs <- solveProblem' $ readFactores xs 
  print rs 
 
readFactores [] = [] 
readFactores (x:y:xs) = (read x, read y):readFactores xs 
 
 
{- 
 
solveet]$ time -f "%E, %M" ./trifacILP 2 5 3 5 5 9 7 3 11 2 13 6 17 3 
[246138750,246359750,246503400] 
0:54.09, 8984 
 
-- n = 14947641017139643312500000 
 
-} 

Comenta la solución

Tienes que identificarte para poder publicar tu comentario.